Package reference

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_eigenvectors do not span the full space of the unperturbed Hamiltonian.

Parameters:

hamiltonian (list | dict | BlockSeries | MutableDenseMatrix | Expr) –
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.
- A single Expr containing the second-quantized Hamiltonian.
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. If neither subspace_eigenvectors nor subspace_indices are provided, the BlockSeries is defined with a single block.
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. If neither subspace_eigenvectors nor subspace_indices are provided, the BlockSeries is defined with a single block.
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_eigenvectors 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 input 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 | MutableDenseMatrix | Expr] | ndarray | MutableDenseMatrix | Expr) – 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. If there is only one block, the boolean array may be provided directly without a dictionary. Must be symmetric, and may not have any True values corresponding to matrix elements coupling degenerate eigenvalues. If the Hamiltonian has second-quantized operators, the array may be a sympy matrix instead, and its entries must correspond to sums of all possible operator powers to eliminate, e.g. a + Dagger(a) or a ** k + Dagger(a) ** k.

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. If neither subspace_eigenvectors nor subspace_indices are provided, the BlockSeries is defined with a single block.
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. If neither subspace_eigenvectors nor subspace_indices are provided, the BlockSeries is defined with a single block.
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

pymablock.series.zero = zero#

A sentinel value for missing terms in series.

Supports multiplication or addition by any other object (and returns itself). Any other usage should be handled by other code.

pymablock.series.one = one#

A sentinel value representing the identity operator.

Used for the zeroth order of a unitary transformation. Does not support any methods, and should be handled by other code.

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]

Number ordered form#

class pymablock.number_ordered_form.LadderOp(*args, **hints)[source]#

Bases: Operator

Ladder operator on a lattice.

Notes

Like a BosonOp, but creation and annihilation operators commute. This is useful for simulating Floquet systems. The operator implementation is minimal, and is meant to be used in combination with the NumberOperator class.

Implementation mostly copied from BosonOp.

class pymablock.number_ordered_form.NumberOperator(*args, **hints)[source]#

Bases: HermitianOperator

Number operator for bosonic, fermionic, and spin operators.

Notes

This class is used to simplify expressions with second-quantized operators. We do this ourselves, because sympy does not support this yet.

Construct a number operator for bosonic, ladder, fermionic, or spin mode.

Parameters:

operator – Operator that the number operator counts.
args – Length-1 list with the operator.
hints – Unused; required for compatibility with sympy.

doit(**hints)[source]#

Evaluate the operator.

For example,

>>> from sympy.physics.quantum.boson import BosonOp
>>> from pymablock.second_quantization import NumberOperator
>>> b = BosonOp('b')
>>> n = NumberOperator(b)
>>> n.doit()
Dagger(b)*b

Returns:: result – The evaluated operator.
Return type:: Expr

property name: Symbol#: Return the name of the operator.

class pymablock.number_ordered_form.NumberOrderedForm(operators, terms, *, validate=True, **hints)[source]#

Bases: Operator

Number ordered form of quantum operators.

A number ordered form represents quantum operators where: 1. All creation operators are on the left 2. All annihilation operators are on the right 3. Number operators (and other scalar expressions) are in the middle

This representation makes it easy to manipulate complex quantum expressions, because commuting a creation or annihilation operator through a function of a number operator simply replaces the corresponding number operator N with N ± 1.

See the second quantization documentation for a detailed description.

Parameters:

operators (Sequence[BosonOp | LadderOp | SigmaOpBase | FermionOp]) – List of quantum operators (annihilation operators only).
terms (dict[tuple[int | Integer, ...], Expr] | tuple[tuple[tuple[int | Integer, ...], Expr], ...] | Tuple) – Dictionary mapping operator power tuples to coefficient expressions. Negative powers represent creation operators, positive powers represent annihilation operators.
**hints (dict) – Additional hints passed to the parent class.
validate (bool)

applyfunc(func, *args, **kwargs)[source]#

Apply a SymPy function to the terms of this NumberOrderedForm.

Notes

NumberOrderedForm stores its coefficients with number operators replaced with integer placeholder symbols. This method does applies the function to those coefficients, but does not change the operators themselves.

Parameters:

func (Callable) – SymPy function to apply (e.g., sympy.simplify, sympy.factor)
*args – Additional positional arguments for the function
**kwargs – Additional keyword arguments for the function

Returns:

A new NumberOrderedForm with the function applied to its terms

Return type:

NumberOrderedForm

as_expr()[source]#

Convert the NumberOrderedForm to a standard SymPy expression.

Returns:: result – A standard SymPy expression equivalent to this NumberOrderedForm.
Return type:: Expr

Examples

Convert a NumberOrderedForm to a standard SymPy expression:

>>> from sympy.physics.quantum.boson import BosonOp
>>> from pymablock.number_ordered_form import (
...     NumberOrderedForm, NumberOperator
... )
>>> a = BosonOp('a')
>>> # Create NumberOrderedForm with creation and annihilation operators
>>> nof = NumberOrderedForm.from_expr(a.adjoint() * a + 2)
>>> nof
2 + N_a
>>> # Convert back to a standard SymPy expression
>>> expr = nof.as_expr()
>>> expr
2 + N_a
>>> # You can also use as_expr() with number operators
>>> nof2 = NumberOrderedForm.from_expr(NumberOperator(a) + 3)
>>> nof2.as_expr()
3 + N_a

doit(**hints)[source]#

Evaluate the NumberOrderedForm.

Parameters:: **hints – Additional hints passed to the parent class.
Returns:: result – The evaluated NumberOrderedForm.
Return type:: Expr

filter_terms(conditions, keep=False)[source]#

Filter the terms of this NumberOrderedForm based on given conditions.

Parameters:

conditions (tuple[tuple[Expr, ...], ...]) – Tuples of conditions to filter the terms. Each condition is a tuple containing the powers of operators, possibly symbolic, e.g. 3 + n.
keep (bool) – If True, keep the terms that satisfy any of the conditions. If False (default), keep the terms that do not satisfy any of the conditions.

Returns:

A new NumberOrderedForm with only the terms that do not satisfy any of the conditions.

Return type:

NumberOrderedForm

classmethod from_expr(expr, operators=None)[source]#

Create a NumberOrderedForm instance from a sympy expression.

Parameters:

expr (Expr) – Sympy expression with quantum operators.
operators (Sequence[BosonOp | LadderOp | SigmaOpBase | FermionOp] | None) – List of quantum operators to use. If None, operators will be extracted from the expression.

Returns:

A NumberOrderedForm instance representing the expression.

Return type:

NumberOrderedForm

Examples

Create a NumberOrderedForm from a sympy expression with bosonic operators:

>>> from sympy.physics.quantum import boson
>>> from pymablock.number_ordered_form import NumberOrderedForm
>>> a = BosonOp('a')
>>> expr = a.adjoint() * a + 1  # a^† * a + 1
>>> nof = NumberOrderedForm.from_expr(expr)
>>> nof
1 + N_a

Using a number operator:

>>> from pymablock.number_ordered_form import NumberOperator
>>> n_a = NumberOperator(a)
>>> expr = n_a + 2
>>> nof = NumberOrderedForm.from_expr(expr)
>>> nof
2 + N_a

property is_algebraic#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_complex#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_composite#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_even#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_extended_negative#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_extended_nonnegative#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_extended_nonpositive#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_extended_nonzero#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_extended_positive#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_extended_real#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_imaginary#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_integer#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_irrational#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_negative#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_noninteger#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_nonnegative#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_nonpositive#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_nonzero#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_odd#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

is_particle_conserving()[source]#

Check if this expression conserves particle numbers.

Returns:: True if the expression has no unpaired creation or annihilation operators, False otherwise.
Return type:: bool

property is_positive#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_prime#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_rational#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_real#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_transcendental#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property is_zero#: Returns True when the argument is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

property operators: list[BosonOp | LadderOp | SigmaOpBase | FermionOp]#: The list of included operators.

property terms: dict[tuple[int | Integer, ...], Expr] | tuple[tuple[tuple[int | Integer, ...], Expr], ...] | Tuple#

The dictionary of terms.

Notes

Internally, terms are stored as a tuple of (key, value) tuples for better performance. This property converts the internal representation to a dictionary for compatibility.

pymablock.number_ordered_form.find_operators(expr)[source]#

Find all quantum operators in a SymPy expression.

Parameters:: expr (Expr) – The expression to search for quantum operators.
Returns:: operators – A list of unique quantum operators found in the expression. Boson and spin operators are listed before fermion operators and both are sorted by their names.
Return type:: list[OperatorType]

Second quantization#

Second quantization tools.

See number_ordered_form_plan.md for the plan to implement NumberOrderedForm as an Operator subclass for better representation of number-ordered expressions.

pymablock.second_quantization.apply_mask_to_operator(operator, mask, keep=True)[source]#

Apply a mask to filter specific terms in a matrix operator.

This function selectively keeps terms in a symbolic matrix operator based on their powers of creation and annihilation operators.

See more details of how this works in the second quantization documentation.

Parameters:

operator (MatrixBase) – Matrix operator containing symbolic expressions with second quantized operators.
mask (ndarray) – A matrix with NumberOrderedForm elements that define selection criteria. Specifically, the elements of the operator[i, j] with powers matching any mask[i, j].terms are selected.
keep (bool) – If True (default), keep the terms that satisfy any of the conditions. If False discard the terms that satisfy any of the conditions. Used for inverting the mask.

Returns:

filtered – A new matrix with the same shape as the input, but containing only the selected terms.

Return type:

sympy.matrices.dense.MutableDenseMatrix

Examples

Let’s filter out terms in a Hamiltonian matrix based on their boson number operators:

>>> import sympy
>>> from sympy.physics.quantum import boson, Dagger
>>> from pymablock.number_ordered_form import NumberOrderedForm
>>> from pymablock.second_quantization import apply_mask_to_operator
>>>
>>> # Create bosonic operators
>>> a = boson.BosonOp('a')
>>> b = boson.BosonOp('b')
>>>
>>> # Create a matrix with different operator terms
>>> H = sympy.Matrix([[a * Dagger(a) + b * Dagger(b), a * Dagger(b)],
        [b * Dagger(a), a * Dagger(a) - b * Dagger(b)]])
>>> # Convert to NumberOrderedForm for easier handling
>>> H_nof = H.applyfunc(NumberOrderedForm.from_expr)
>>>
>>> # Create a mask that selects only terms with a specific power pattern
>>> # Select only terms with exactly one 'a' operator and one 'b' operator
>>> mask = sympy.Matrix([[sympy.S.Zero, NumberOrderedForm([a, b], {(1, -1): sympy.S.One})],
            [NumberOrderedForm([a, b], {(1, 1): sympy.S.One}), sympy.S.Zero]])
>>> H_filtered = apply_mask_to_operator(H_nof, mask, keep=True)
>>> # H_filtered now contains only the terms that match the mask:
>>> # [[0, Dagger(b)a], [0, 0]]

pymablock.second_quantization.solve_sylvester_2nd_quant(eigs)[source]#

Solve a Sylvester equation for 2nd quantized diagonal Hamiltonians.

Parameters:: eigs (tuple[tuple[Expr, ...], ...]) – Tuple of lists of expressions representing the diagonal Hamiltonian blocks.
Returns:: A function that takes a matrix of operators and a tuple of indices, and computes the element-wise solution to the Sylvester equation for those diagonal Hamiltonian blocks.
Return type:: Callable

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#

Number ordered form#

Second quantization#

Kernel polynomial method (KPM)#

Algorithms#