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, fully_diagonalize=())[source]#

Find the block diagonalization of a Hamiltonian order by order.

This uses a generalization of quasi-degenerate perturbation theory known as Lowdin perturbation theory, Schrieffer-Wolff transformation, or van Vleck transformation to the case of multiple blocks. Some blocks of the resulting Hamiltonian can be fully diagonalized, reproducing the usual Rayleigh-Schrodinger perturbation theory. Alternatively, the algorithm can perturbatively eliminate any subset of the off-diagonal elements of the Hamiltonian.

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.

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 into blocks. 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 have the same format as in the list case.
- 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,
  used directly.
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 sets of orthonormal eigenvectors to project the Hamiltonian onto and separate it into blocks. If None, the unperturbed Hamiltonian must be block diagonal. If some vectors are missing, the implicit method is used. Mutually exclusive with subspace_indices.
subspace_indices (tuple[int, ...] | ndarray | None) – An array indicating which state belongs to which subspace. If there are two blocks, 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.
fully_diagonalize (tuple[int] | dict[int, ndarray]) – Indices of the blocks that should be fully diagonalized. If the Hamiltonian only has one block, it is fully diagonalized by default. Alternatively can be a dictionary with the indices of the diagonal blocks as keys and as values appropriately shaped numpy boolean arrays marking which matrix elements should be eliminated by the diagonalization. Must be symmetric, and may not have any True values corresponding to matrix elements coupling degenerate eigenvalues.

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]

pymablock.operator_to_BlockSeries(operator, *, name=None, subspace_eigenvectors=None, subspace_indices=None, implicit=False, symbols=None, atol=1e-12, hermitian=False)[source]#

Convert an operator to BlockSeries format.

The normalization consists of the following steps:

Determine perturbation parameters from the input format,
Separate the operator into blocks based on eigenvectors or subspace indices,
Form a BlockSeries.

Parameters:

operator (list | dict | BlockSeries | MutableDenseMatrix) –
Full symbolic or numeric operator. It is normalized to a BlockSeries and separated into blocks corresponding to different subspaces.

Supported formats:
- A list,
  of the form [h_0, h_1, h_2, ...] where h_0 is the unperturbed operator and h_1, h_2, ... are the first order perturbations. The elements h_i may be Matrix, ndarray, spmatrix, that require separating into blocks. Otherwise, h_i may be a list of lists with the operator 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 operators. 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 operator.dimension_names. h_i have the same format as in the list case.
- 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,
  used directly.
name (str | None) – Name of the operator.
subspace_eigenvectors (tuple[Any, Any] | None) – A tuple with sets of orthonormal eigenvectors to project the operator onto and separate it into blocks. If some vectors are missing, the last block is defined implicitly. Mutually exclusive with subspace_indices.
subspace_indices (tuple[int, ...] | None) – An array indicating which state belongs to which subspace. If there are two blocks, the labels are 0 for the A (effective) subspace and 1 for the B (auxiliary) subspace. Mutually exclusive with subspace_eigenvectors.
implicit (bool) – Whether to wrap the operator of the last subspace into a linear operator.
symbols (list[Symbol] | None) – Symbols that label the perturbative parameters of a symbolic operator. The order of the symbols is mapped to the indices of the operator, see BlockSeries. If None, the perturbative parameters are taken from the unperturbed operator.
atol (float) – Absolute tolerance to consider matrices as exact zeros. This is used to validate that the unperturbed operator is block-diagonal.
hermitian (bool) – Whether the operator may be assumed to be Hermitian (for performance).

Returns:

operator – Transformed operator.

Return type:

BlockSeries

Solvers of Sylvester equation#

Algorithms for quasi-degenerate perturbation theory.

pymablock.block_diagonalization.solve_sylvester_KPM(h_0, subspace_eigenvectors, solver_options=None)[source]#

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.
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.
- auxiliary_vectors: np.ndarray
  Partial set of eigenvectors of the auxiliary subspace, used to speed up convergence of the KPM solver.

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, vecs_implicit=None, atol=1e-12)[source]#

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 (tuple[ndarray | MatrixBase, ...]) – List of arrays of eigenvalues in each subspace. If vecs_implicit are provided, the last subspace corresponds to the partial set of eigenvalues of the implicit subspace.
vecs_implicit (ndarray | None) – Partial eigenvectors of the implicit subspace. Only used for KPM. If provided, the last block is returned in the basis of the original Hilbert space rather than eigenbasis of the Hamiltonian.
atol (float) – Absolute tolerance to consider energy differences as exact zeros.

Returns:

solve_sylvester – Function that solves Sylvester’s equation.

Return type:

Callable

pymablock.block_diagonalization.solve_sylvester_direct(h_0, eigenvectors, **solver_options)[source]#

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 (list[ndarray]) – Eigenvectors of the effective subspaces 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)[source]#

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)[source]#

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)[source]#

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)[source]#

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)[source]#

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

Parameters:

hamiltonian (ndarray | spmatrix) – Hamiltonian of the system.
eps (float) – 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]]

Algorithms#

Tools for compiling optimized series computations.

pymablock.algorithm_parsing.series_computation(series, algorithm, scope=None, *, operator=None)[source]#

Compile a BlockSeries computation.

Given several series, functions to apply to their elements, and an algorithm, return the output series defined by the algorithm.

The algorithm parsing used by this function applies multiple optimizations used in the Pymablock main algorithm. While these could be generated by hand, the resulting code is complex and repetitive. The mini-language used to specify the algorithm allows to avoid this complexity and enabled multiple improvements of the current algorithm and simplifies development of new algorithms.

Specifically, code generation from the mini-language description of an algorithm:

Handles initialization of series and the definition of their evaluation functions.
Utilizes hermiticity and antihermiticity to reduce the number of evaluations.
Automatically handles the implicit mode when parts of the series are provided as linear operators.
Handles deletion of intermediate series terms that are only used once to reduce the memory usage.

Implementing a new algorithm is advanced usage, and familiarity with the codebase is highly recommended.

Parameters:

series (dict[str, BlockSeries]) – Dictionary with all input series, where the keys are the names of the series.
algorithm (Callable) – Algorithm to use for the block diagonalization. Should be passed as a callable whose contents follow the algorithm mini-language, see notes below.
scope (dict | None) – Extra variables to pass to pass to the algorithm. It is particularly relevant for passing custom functions or data.
operator (Callable | None) – (optional) function to use for matrix multiplication. Defaults to matmul.

Returns:

series (dict[str, BlockSeries]) – A dictionary with all the series used in the computation. The keys are the names of the series and the values are the corresponding BlockSeries.
linear_operator_series (dict[str, BlockSeries]) – The same series as above, but wrapped into linear operators. Only used in the implicit mode.

Return type:

tuple[dict[str, BlockSeries], dict[str, BlockSeries]]

Notes

The algorithm callable is not evaluated directly, but rather parsed to extract the computation that needs to be performed. It needs to follow the specification below.

Warning

This domain-specific language is experimental and may change in the future.

The function body contains multiple with statements that define the series and products of that algorithm. Throughout the definition the series and products are represented by their name using string literals.

A series definition allows the following statements:

start = ... to define the zeroth order of the series. Allowed values are "series_name", 0, 1.
hermitian or antihermitian to optionally mark the lower triangular blocks of a series to be evaluated using a conjugate transpose of the upper triangular blocks.
One or more expressions that define how to evaluate the series. If there are multiple expressions, they are summed together. The expression can contain the following:
- String literals to represent series.
- Attribute .adj access to represent the conjugate adjoint of a series.
- Integer literals.
- Unary and binary operations.
- Function calls. Using f("series") will call the function f with the series and the block index as arguments. Using f(expression) will call the function with the evaluated expression and block index as arguments.
if <condition>: differentiates evaluation based on the requested index. Allowed conditions are:
- diagonal: indices on the main diagonal.
- offdiagonal: indices not on the main diagonal.
- lower: indices in the lower triangle.

If a name contains an “@” symbol, it defines a Cauchy product of the terms in it. For example "A @ B @ C" is a Cauchy product of the series A, B, and C.

A product definition must contain one of the two following statements:

hermitian to mark the product as hermitian.
pass otherwise.

The final return statement in the function body defines a tuple of series that are the output of the algorithm and terms of which should not be deleted.

Example

The algorithm definition may look as follows (this example does not do anything useful):

def my_algorithm():
    with "B":
        start = 0
        hermitian
        if diagonal:
            "A" + f("B @ C")

    with "C":
        start = "A"
        if offdiagonal:
            "A" + "B" / 2
        "B @ C"

    with "B @ C":
        hermitian

    return "C"

Here "A" is an input, "B" and "C" are defined in the computation, and the function f must be provided using the scope.

For an extended example, see the main function in pymablock/algorithms.py.

Package reference

Contents

Package reference#

Block diagonalization#

Solvers of Sylvester equation#

Series#

Linear algebra#

Kernel polynomial method (KPM)#

Algorithms#