Package reference

Block diagonalization

Quasi-degenerate perturbation theory

pymablock.block_diagonalize(hamiltonian, *, solve_sylvester=None, subspace_eigenvectors=None, subspace_indices=None, direct_solver=True, solver_options=None, symbols=None, atol=1e-12)

Find the block diagonalization of a Hamiltonian order by order.

This uses quasi-degenerate perturbation theory known as Lowdin perturbation theory, Schrieffer-Wolff transformation, or van Vleck transformation.

This function does not yet perform the computation. Instead, it defines the computation as a BlockSeries object, which can be evaluated at any order.

This function accepts a Hamiltonian in several formats and it first brings it to the eigenbasis of the unperturbed Hamiltonian if the blocks, eigenvectors or indices of the eigenvalues are provided, see below.

For large numerical calculations with a sparse Hamiltonian and a low dimensional relevant subspace, the algorithm uses an implicit representation of the spectrum and does not require full diagonalization. This is enabled if subspace_vectors do not span the full space of the unperturbed Hamiltonian.

Parameters:

• hamiltonian (list | dict | BlockSeries | MutableDenseMatrix) –

  Full symbolic or numeric Hamiltonian to block diagonalize. The Hamiltonian is normalized to a BlockSeries by separating it into effective and auxiliary subspaces.

  Supported formats:

  • A list,

    of the form [h_0, h_1, h_2, ...] where h_0 is the unperturbed Hamiltonian and h_1, h_2, ... are the first order perturbations. The elements h_i may be Matrix, ndarray, spmatrix, that require separating the unperturbed Hamiltonian into effective and auxiliary subspaces. Otherwise, h_i may be a list of lists with the Hamiltonian blocks.

  • A dictionary,

    of the form {(0, 0): h_0, (1, 0): h_1, (0, 1): h_2}, or {1: h_0, x: h_1, y: h_2} for symbolic Hamiltonians. In the former case, the keys must be tuples of integers indicating the order of each perturbation. In the latter case, the keys must be monomials and the indices are ordered as in H.dimension_names. The values of the dictionary, h_i may be Matrix, ndarray, spmatrix, that require separating the unperturbed Hamiltonian into effective and auxiliary subspaces. Otherwise, h_i may be a list of lists with the Hamiltonian blocks.

  • A Matrix,

    unless a list of symbols is provided as perturbative parameters, all symbols will be treated as perturbative. The normalization to BlockSeries is done by Taylor expanding in symbols to the desired order.

  • A BlockSeries,

    returned unchanged.

• solve_sylvester (Callable | None) – A function that solves the Sylvester equation. If not provided, it is selected automatically based on the inputs.

• subspace_eigenvectors (tuple[ndarray | MutableDenseMatrix, ...] | None) – A tuple with orthonormal eigenvectors to project the Hamiltonian in and separate it into the A (effective) and B (auxiliary) blocks. The first element of the tuple has the eigenvectors of the A subspace, and the second element has the eigenvectors of the B subspace. If None, the unperturbed Hamiltonian must be block diagonal. For implicit, the (partial) auxiliary subspace may be missing or incomplete. Mutually exclusive with subspace_indices.

• subspace_indices (tuple[int, ...] | ndarray | None) – An array indicating which basis vector belongs to which subspace. The labels are 0 for the A (effective) subspace and 1 for the B (auxiliary) subspace. Only applicable if the unperturbed Hamiltonian is diagonal. Mutually exclusive with subspace_eigenvectors.

• solver_options (dict | None) – Dictionary containing the options to pass to the Sylvester solver. See docstrings of solve_sylvester_KPM and solve_sylvester_direct for details.

• direct_solver (bool) – Whether to use the direct solver that relies on MUMPS (default). Otherwise, the KPM solver is used. Only applicable if the implicit method is used (i.e. subspace_vectors is incomplete)

• symbols (Symbol | Sequence[Symbol] | None) – Symbols that label the perturbative parameters of a symbolic Hamiltonian. The order of the symbols is mapped to the indices of the Hamiltonian, see BlockSeries. If None, the perturbative parameters are taken from the unperturbed Hamiltonian.

• atol (float) – Absolute tolerance to consider matrices as exact zeros. This is used to validate that the unperturbed Hamiltonian is block-diagonal.

Returns:

• H_tilde (BlockSeries) – Block diagonalized Hamiltonian.

• U (BlockSeries) – Unitary matrix that block diagonalizes H such that U * H * U^H = H_tilde.

• U_adjoint (BlockSeries) – Adjoint of U.

Return type:

tuple[BlockSeries, BlockSeries, BlockSeries]

Solvers of Sylvester equation

Algorithms for quasi-degenerate perturbation theory.

pymablock.block_diagonalization.solve_sylvester_KPM(h_0, subspace_eigenvectors, solver_options=None)

Solve Sylvester energy division for the Kernel Polynomial Method (KPM).

General energy division for numerical problems through either full knowledge of the B-space or application of the KPM Green’s function.

Parameters:

• h_0 (ndarray | spmatrix) – Unperturbed Hamiltonian of the system.

• subspace_eigenvectors (tuple[ndarray, ...]) – Subspaces to project the unperturbed Hamiltonian and separate it into blocks. The first element of the tuple contains the effective subspace, and the second element contains the (partial) auxiliary subspace.

• solver_options (dict | None) –

  Dictionary containing any of the following options for KPM.

  • eps: float

    Tolerance for Hamiltonian rescaling.

  • atol: float

    Accepted precision of the Green’s function result in 2-norm.

  • max_moments: int

    Maximum number of expansion moments of the Green’s function.

Returns:

solve_sylvester – Function that applies divide by energies to the right hand side of Sylvester’s equation.

Return type:

Callable

pymablock.block_diagonalization.solve_sylvester_diagonal(eigs_A, eigs_B, vecs_B=None)

Define a function for solving a Sylvester’s equation for diagonal matrices.

Optionally, this function also applies the eigenvectors of the second matrix to the solution.

Parameters:

• eigs_A (ndarray | MatrixBase) – Eigenvalues of the effective (A) subspace of the unperturbed Hamiltonian.

• eigs_B (ndarray | MatrixBase) – Eigenvalues of auxiliary (B) subspace of the unperturbed Hamiltonian.

• vecs_B (ndarray | None) – Eigenvectors of the auxiliary (B) subspace of the unperturbed Hamiltonian.

Returns:

solve_sylvester – Function that solves Sylvester’s equation.

Return type:

Callable

pymablock.block_diagonalization.solve_sylvester_direct(h_0, eigenvectors, **solver_options)

Solve Sylvester equation using a direct sparse solver.

This function uses MUMPS, which is a parallel direct solver for sparse matrices. This solver is very efficient for large sparse matrices.

Parameters:

• h_0 (spmatrix) – Unperturbed Hamiltonian of the system.

• eigenvectors (ndarray) – Eigenvectors of the effective subspace of the unperturbed Hamiltonian.

• **solver_options (dict) – Keyword arguments to pass to the solver eps and atol, see pymablock.linalg.direct_greens_function.

Returns:

solve_sylvester – Function that solves the corresponding Sylvester equation.

Return type:

Callable[[numpy.ndarray], numpy.ndarray]

Series

class pymablock.series.BlockSeries(eval=None, data=None, shape=(), n_infinite=1, dimension_names=None, name=None)

Series with block structure and caching of evaluated items.

BlockSeries is used to represent Hamiltonians, unitary transformations, and other series needed in the block diagonalization algorithm.

Construct an infinite series that caches its items.

The series has finite and infinite dimensions.

Parameters:

• eval (Callable | None) – Function that takes an index and returns the corresponding item.

• data (dict[tuple[int, ...], Any] | None) – Dictionary {index: value} with initial data of the series.

• shape (tuple[int, ...]) – Shape of the finite dimensions.

• n_infinite (int) – Number of infinite dimensions.

• dimension_names (tuple[str | Symbol, ...] | None) – Names of the infinite dimensions. This is used for symbolic series.

• name (str | None) – Name of the series. Used for printing.

Series with block structure and their product operations.

pymablock.series.cauchy_dot_product(*series, operator=None, hermitian=False)

Multivariate Cauchy product of BlockSeries.

Notes: This treats a singleton one as the identity operator.

Parameters:

• series (BlockSeries) – Series to multiply using their block structure. Must be at least two.

• operator (Callable | None) – (optional) Function for multiplying elements of the series. Default is matrix multiplication matmul.

• hermitian (bool) – (optional) whether to assume that the result is Hermitian in order to reduce computational costs.

Returns:

A new series that is the Cauchy dot product of the given series.

Return type:

BlockSeries

Linear algebra

Linear algebra utilities.

pymablock.linalg.direct_greens_function(h, E, atol=1e-07, eps=0)

Compute the Green’s function of a Hamiltonian using MUMPS solver.

Parameters:

• h (spmatrix) – Hamiltonian matrix.

• E (float) – Energy at which to compute the Green’s function.

• atol (float) – Accepted precision of the desired result in infinity-norm.

• eps (float) – Tolerance for the MUMPS solver to identify null pivots. Passed through to MUMPS CNTL(3), see MUMPS user guide.

Returns:

greens_function – Function that solves \((E - H) sol = vec\).

Return type:

Callable[[np.ndarray], np.ndarray]

Kernel polynomial method (KPM)

Kernel Polynomial Method (KPM).

See arXiv:cond-mat/0504627 and arXiv:1909.09649.

pymablock.kpm.greens_function(hamiltonian, energy, vector, atol=1e-07, max_moments=1000000)

Return a solution of (energy - hamiltonian) @ x = vector.

Uses the Kernel polynomial method (KPM) with the Jackson kernel.

Parameters:

• hamiltonian (ndarray | spmatrix) – Hamiltonian with shape (N, N). It is required that the Hamiltonian is rescaled so that its spectrum lies in the interval [-1, 1].

• energy (float) – Rescaled energy at which to evaluate the Green’s function.

• vector (ndarray) – Vector with shape (N,).

• atol (float) – Accepted precision of the desired result in 2-norm.

• max_moments (int) – Maximum order of KPM expansion to compute.

Returns:

solution – Solution x of (E - H_0) * x = v.

Return type:

ndarray

pymablock.kpm.rescale(hamiltonian, eps=0.01, bounds=None, lower_bounds=None)

Rescale a Hamiltonian to the interval [-1 - eps/2, 1 + eps/2].

Adapted with modifications from kwant.kpm Copyright 2011-2016 Kwant developers, BSD simplified license https://gitlab.kwant-project.org/kwant/kwant/-/blob/v1.4.3/LICENSE.rst

Parameters:

• hamiltonian (ndarray | spmatrix) – Hamiltonian of the system.

• eps (float | None) – Extra tolerance to add to the spectral bounds as a fraction of the bandwidth.

• bounds (tuple[float, float] | None) – Estimated boundaries of the spectrum. If not provided the maximum and minimum eigenvalues are calculated.

• lower_bounds (tuple[float, float] | None) – Energy interval to definitely include within the [-1, 1] rescaled energies.

Returns:

• hamiltonian_rescaled – Rescaled Hamiltonian.

• (a, b) – Rescaling parameters such that h_rescaled = (h - b) / a.

Return type:

tuple[ndarray | spmatrix, tuple[float, float]]