ffsim.linalg¶

Linear algebra utilities.

class ffsim.linalg.GivensRotation(c: float, s: complex, i: int, j: int)[source]¶

Bases: NamedTuple

A Givens rotation.

A Givens rotation acts on the two-dimensional subspace spanned by the \(i\)-th and \(j\)-th basis vectors as

\[\begin{split}\begin{pmatrix} c & s \\ -s^* & c \\ \end{pmatrix}\end{split}\]

where \(c\) is a real number and \(s\) is a complex number.

c: float¶: Alias for field number 0

i: int¶: Alias for field number 2

j: int¶: Alias for field number 3

s: complex¶: Alias for field number 1

ffsim.linalg.apply_matrix_to_slices(target, mat, slices, *, out=None)[source]¶

Apply a matrix to slices of a target tensor.

Parameters:

target (ndarray) – The tensor containing the slices on which to apply the matrix.
mat (ndarray) – The matrix to apply to slices of the target tensor.
slices – The slices of the target tensor on which to apply the matrix.

Return type:

ndarray

Returns:

The resulting tensor.

ffsim.linalg.double_factorized(two_body_tensor, *, tol=1e-08, max_vecs=None, optimize=False, method='L-BFGS-B', callback=None, options=None, diag_coulomb_indices=None, return_optimize_result=False, cholesky=True)[source]¶

Double-factorized decomposition of a two-body tensor.

The double-factorized decomposition is a representation of a two-body tensor \(h_{pqrs}\) as

\[h_{pqrs} = \sum_{t=0}^{L - 1} \sum_{k\ell} Z^{(t)}_{k\ell} U^{(t)}_{pk} U^{(t)}_{qk} U^{(t)}_{r\ell} U^{(t)}_{s\ell},\]

where each \(Z^{(t)}\) is a real symmetric matrix representing a diagonal Coulomb operator and each \(U^{(t)}\) is a unitary matrix representing an orbital rotation.

The number of terms \(L\) in the decomposition depends on the allowed error threshold. A larger error threshold may yield a smaller number of terms. Furthermore, the max_vecs parameter specifies an optional upper bound on \(L\). The max_vecs parameter is always respected, so if it is too small, then the error of the decomposition may exceed the specified error threshold.

The default behavior of this routine is to perform a straightforward “exact” factorization of the two-body tensor based on a nested eigenvalue decomposition, and then truncate the terms based on the values of tol and max_vecs. If optimize is set to True, then the entries of the resulting tensors (the diagonal Coulomb matrices and orbital rotations) are further optimized with scipy.optimize.minimize to reduce the error in the factorization. The diagonal Coulomb matrices returned by the optimization can be optionally constrained to have only certain elements allowed to be nonzero. This is achieved by passing the diag_coulomb_indices parameter, which is a list of matrix entry indices (integer pairs) specifying where the diagonal Coulomb matrices are allowed to be nonzero. Since the diagonal Coulomb matrices are symmetric, only upper triangular indices should be given, i.e., pairs \((i, j)\) where \(i \leq j\).

References

Note: Currently, only real-valued two-body tensors are supported.

Parameters:

two_body_tensor (ndarray) – The two-body tensor to decompose.
tol (float) – Tolerance for error in the decomposition. The error is defined as the maximum absolute difference between an element of the original tensor and the corresponding element of the reconstructed tensor.
max_vecs (int | None) – An optional limit on the number of terms to keep in the decomposition of the two-body tensor. This argument overrides tol.
optimize (bool) – Whether to optimize the tensors returned by the decomposition to to minimize the error in the factorization.
method (str) – The optimization method. See the documentation of scipy.optimize.minimize for possible values. This argument is ignored if optimize is set to False.
callback – Callback function for the optimization. See the documentation of scipy.optimize.minimize for usage. This argument is ignored if optimize is set to False.
options (dict | None) – Options for the optimization. See the documentation of scipy.optimize.minimize for usage. This argument is ignored if optimize is set to False.
diag_coulomb_indices (list[tuple[int, int]] | None) – Allowed indices for nonzero values of the diagonal Coulomb matrices. Matrix entries corresponding to indices not in this list will be set to zero. This list should contain only upper triangular indices, i.e., pairs \((i, j)\) where \(i \leq j\). This argument is ignored if optimize is set to False.
return_optimize_result (bool) – Whether to also return the OptimizeResult returned by scipy.optimize.minimize. This argument is ignored if optimize is set to False.
cholesky (bool) – Whether to perform the factorization using a modified Cholesky decomposition. If False, a full eigenvalue decomposition is used instead, which can be much more expensive. This argument is ignored if optimize is set to True.

Returns:

The diagonal Coulomb matrices and the orbital rotations. Each list of matrices is collected into a Numpy array, so this method returns a tuple of two Numpy arrays, the first containing the diagonal Coulomb matrices and the second containing the orbital rotations. Each Numpy array will have shape (L, n, n) where L is the number of terms in the decomposition and n is the number of orbitals. If optimize and return_optimize_result are both set to True, the OptimizeResult returned by scipy.optimize.minimize is also returned.

Raises:

ValueError – diag_coulomb_indices contains lower triangular indices.

ffsim.linalg.double_factorized_t2(t2_amplitudes, *, tol=1e-08, max_terms=None, optimize=False, method='L-BFGS-B', callback=None, options=None, diag_coulomb_indices=None, regularization=0, multi_stage_start=None, multi_stage_step=None, return_optimize_result=False)[source]¶

Double-factorized decomposition of \(t_2\) amplitudes.

The double-factorized decomposition of a \(t_2\) amplitudes tensor \(t_{ijab}\) is

\[t_{ijab} = i \sum_{k=0}^{L - 1} \sum_{pq} Z^{(k)}_{pq} U^{(k)}_{ap} U^{(k)*}_{ip} U^{(k)}_{bq} U^{(k)*}_{jq},\]

where each \(Z^{(k)}\) is a real symmetric matrix representing a diagonal Coulomb operator, each \(U^{(k)}\) is a unitary matrix representing an orbital rotation, \(i\) and \(j\) run over occupied orbitals, and \(a\) and \(b\) run over virtual orbitals.

The number of terms \(L\) in the decomposition depends on the allowed error threshold. A larger error threshold may yield a smaller number of terms. Furthermore, the max_terms parameter specifies an optional upper bound on \(L\). The max_terms parameter is always respected, so if it is too small, then the error of the decomposition may exceed the specified error threshold.

The default behavior of this routine is to perform a straightforward “exact” factorization of the \(t_2\) amplitudes tensor based on a nested eigenvalue decomposition, and then truncate the terms based on the values of tol and max_terms. If optimize is set to True, then the entries of the resulting tensors (the diagonal Coulomb matrices and orbital rotations) are further optimized with scipy.optimize.minimize to reduce the error in the factorization. The diagonal Coulomb matrices returned by the optimization can be optionally constrained to have only certain elements allowed to be nonzero. This is achieved by passing the diag_coulomb_indices parameter, which is a list of matrix entry indices (integer pairs) specifying where the diagonal Coulomb matrices are allowed to be nonzero. Since the diagonal Coulomb matrices are symmetric, only upper triangular indices should be given, i.e., pairs \((i, j)\) where \(i \leq j\).

A “multi-stage” optimization can be requested by passing either multi_stage_start or multi_stage_step. In multi-stage optimization, the number of terms is iteratively reduced by multi_stage_step, starting from multi_stage_start, until there are only max_terms remaining. This procedure can yield a more accurate factorization at the cost of increased running time. If you only pass one of multi_stage_start or multi_stage_step, then the value of the other is chosen automatically: multi_stage_start defaults to the total number of terms returned by the exact factorization, and multi_stage_step defaults to 1.

The original \(t_2\) amplitudes tensor can be reconstructed, up to the error tolerance, using the following function:

def reconstruct_t2(
    diag_coulomb_mats: np.ndarray, orbital_rotations: np.ndarray, nocc: int
) -> np.ndarray:
    return (
        1j
        * np.einsum(
            "kpq,kap,kip,kbq,kjq->ijab",
            diag_coulomb_mats,
            orbital_rotations,
            orbital_rotations.conj(),
            orbital_rotations,
            orbital_rotations.conj(),
        )[:nocc, :nocc, nocc:, nocc:]
    )

Note: Currently, only real-valued \(t_2\) amplitudes are supported.

Parameters:

t2_amplitudes (ndarray) – The doubles amplitudes tensor of shape (nocc, nocc, nvrt, nvrt), where nocc is the number of occupied orbitals and nvrt is the number of virtual orbitals.
tol (float) – Tolerance for error in the decomposition. The error is defined as the maximum absolute difference between an element of the original tensor and the corresponding element of the reconstructed tensor.
max_terms (int | None) – An optional limit on the number of terms to keep in the decomposition of the \(t_2\) amplitudes tensor. This argument overrides tol.
optimize (bool) – Whether to optimize the tensors returned by the decomposition to to minimize the error in the factorization.
method (str) – The optimization method. See the documentation of scipy.optimize.minimize for possible values. This argument is ignored if optimize is set to False.
callback – Callback function for the optimization. See the documentation of scipy.optimize.minimize for usage. This argument is ignored if optimize is set to False.
options (dict | None) – Options for the optimization. See the documentation of scipy.optimize.minimize for usage. This argument is ignored if optimize is set to False.
diag_coulomb_indices (list[tuple[int, int]] | None) – Allowed indices for nonzero values of the diagonal Coulomb matrices. Matrix entries corresponding to indices not in this list will be set to zero. This list should contain only upper triangular indices, i.e., pairs \((i, j)\) where \(i \leq j\). This argument is ignored if optimize is set to False.
regularization (float) – Parameter for regularizing the norm of the diagonal Coulomb matrices. This specifies the coefficient to a term added to the loss function equal to the difference between the sums of the squares of the Frobenius norms of the diagonal Coulomb matrices from the optimized factorization and the exact factorization. This argument is ignored if optimize is set to False.
multi_stage_start (int | None) – The number of terms to start multi-stage optimization. This argument is ignored if optimize is set to False.
multi_stage_step (int | None) – The number of terms to eliminate in each stage of multi-stage optimization. This argument is ignored if optimize is set to False.
return_optimize_result (bool) – Whether to also return the OptimizeResult returned by scipy.optimize.minimize. This argument is ignored if optimize is set to False.

Return type:

tuple[ndarray, ndarray] | tuple[ndarray, ndarray, OptimizeResult]

Returns:

ffsim.linalg.double_factorized_t2_alpha_beta(t2_amplitudes, *, tol=1e-08, max_terms=None)[source]¶

Double-factorized decomposition of alpha-beta \(t_2\) amplitudes.

Decompose alpha-beta \(t_2\) amplitudes into diagonal Coulomb matrices with orbital rotations. This function returns two arrays:

diagonal_coulomb_mats, with shape (n_terms, 3, norb, norb).
orbital_rotations, with shape (n_terms, 2, norb, norb).

The value of n_terms depends on the error tolerance tol. A larger error tolerance might yield a smaller value for n_terms. You can also set an optional upper bound on n_terms using the max_terms argument.

The original \(t_2\) amplitudes tensor can be reconstructed, up to the error tolerance, using the following function:

def reconstruct_t2_alpha_beta(
    diag_coulomb_mats: np.ndarray,
    orbital_rotations: np.ndarray,
    norb: int,
    nocc_a: int,
    nocc_b: int,
) -> np.ndarray:
    n_terms = diag_coulomb_mats.shape[0]
    expanded_diag_coulomb_mats = np.zeros((n_terms, 2 * norb, 2 * norb))
    expanded_orbital_rotations = np.zeros(
        (n_terms, 2 * norb, 2 * norb), dtype=complex
    )
    for k in range(n_terms):
        (mat_aa, mat_ab, mat_bb) = diag_coulomb_mats[k]
        expanded_diag_coulomb_mats[k] = np.block(
            [[mat_aa, mat_ab], [mat_ab.T, mat_bb]]
        )
        orbital_rotation_a, orbital_rotation_b = orbital_rotations[k]
        expanded_orbital_rotations[k] = scipy.linalg.block_diag(
            orbital_rotation_a, orbital_rotation_b
        )
    return (
        2j
        * np.einsum(
            "kpq,kap,kip,kbq,kjq->ijab",
            expanded_diag_coulomb_mats,
            expanded_orbital_rotations,
            expanded_orbital_rotations.conj(),
            expanded_orbital_rotations,
            expanded_orbital_rotations.conj(),
        )[:nocc_a, norb : norb + nocc_b, nocc_a:norb, norb + nocc_b :]
    )

Note: Currently, only real-valued \(t_2\) amplitudes are supported.

Parameters:

t2_amplitudes (ndarray) – The doubles amplitudes tensor of shape (nocc_a, nocc_b, nvrt_a, nvrt_b), where nocc_a is the number of occupied spin-alpha orbitals, nvrt_a is the number of virtual spin-alpha orbitals, nocc_b is the number of occupied spin-beta orbitals, and nvrt_b is the number of virtual spin-beta orbitals.
tol (float) – Tolerance for error in the decomposition. The error is defined as the maximum absolute difference between an element of the original tensor and the corresponding element of the reconstructed tensor.
max_terms (int | None) – An optional limit on the number of terms to keep in the decomposition of the \(t_2\) amplitudes tensor. This argument overrides tol.

Return type:

tuple[ndarray, ndarray]

Returns:

The diagonal Coulomb matrices, as a Numpy array of shape (n_terms, 3, norb, norb). The last two axes index the rows and columns of the matrices, and the third from last axis, which has 3 dimensions, indexes the spin interaction type of the matrix: alpha-alpha, alpha-beta, and beta-beta (in that order). The first axis indexes the terms of the decomposition.
The orbital rotations, as a Numpy array of shape (n_terms, 2, norb, norb). The last two axes index the rows and columns of the orbital rotations, and the third from last axis, which has 2 dimensions, indexes the spin sector of the orbital rotation: first alpha, then beta. The first axis indexes the terms of the decomposition.

ffsim.linalg.expm_multiply_taylor(mat, vec, tol=1e-12)[source]¶

Compute expm(mat) @ vec using a Taylor series expansion.

Return type:: ndarray

ffsim.linalg.givens_decomposition(mat)[source]¶

Givens rotation decomposition of a unitary matrix.

The Givens rotation decomposition of an \(n \times n\) unitary matrix \(U\) is given by

\[U = D G_L^* G_{L-1}^* \cdots G_1^*\]

where \(D\) is a diagonal matrix and each \(G_k\) is a Givens rotation. Here, the star \(*\) denotes the element-wise complex conjugate. A Givens rotation acts on the two-dimensional subspace spanned by the \(i\)-th and \(j\)-th basis vectors as

\[\begin{split}\begin{pmatrix} c & s \\ -s^* & c \\ \end{pmatrix}\end{split}\]

where \(c\) is a real number and \(s\) is a complex number. Therefore, a Givens rotation is described by a 4-tuple \((c, s, i, j)\), where \(c\) and \(s\) are the numbers appearing in the rotation matrix, and \(i\) and \(j\) are the indices of the basis vectors of the subspace being rotated. This function always returns Givens rotations with the property that \(i\) and \(j\) differ by at most one, that is, either \(j = i + 1\) or \(j = i - 1\).

The number of Givens rotations \(L\) is at most \(\frac{n (n-1)}{2}\), but it may be less. If we think of Givens rotations acting on disjoint indices as operations that can be performed in parallel, then the entire sequence of rotations can always be performed using at most n layers of parallel operations. The decomposition algorithm is described in [1].

[1] William R. Clements et al. Optimal design for universal multiport interferometers.

Parameters:

mat (ndarray) – The unitary matrix to decompose into Givens rotations.

Return type:

tuple[list[GivensRotation], ndarray]

Returns:

A list containing the Givens rotations \(G_1, \ldots, G_L\). Each Givens rotation is represented as a 4-tuple \((c, s, i, j)\), where \(c\) and \(s\) are the numbers appearing in the rotation matrix, and \(i\) and \(j\) are the indices of the basis vectors of the subspace being rotated.
A Numpy array containing the diagonal elements of the matrix \(D\).

ffsim.linalg.is_antihermitian(mat, *, rtol=1e-05, atol=1e-08)[source]¶

Determine if a matrix is approximately anti-Hermitian.

Parameters:

mat (ndarray) – The matrix.
rtol (float) – Relative numerical tolerance.
atol (float) – Absolute numerical tolerance.

Return type:

bool

Returns:

Whether the matrix is anti-Hermitian within the given tolerance.

ffsim.linalg.is_hermitian(mat, *, rtol=1e-05, atol=1e-08)[source]¶

Determine if a matrix is approximately Hermitian.

Parameters:

mat (ndarray) – The matrix.
rtol (float) – Relative numerical tolerance.
atol (float) – Absolute numerical tolerance.

Return type:

bool

Returns:

Whether the matrix is Hermitian within the given tolerance.

ffsim.linalg.is_orthogonal(mat, *, rtol=1e-05, atol=1e-08)[source]¶

Determine if a matrix is approximately orthogonal.

Parameters:

mat (ndarray) – The matrix.
rtol (float) – Relative numerical tolerance.
atol (float) – Absolute numerical tolerance.

Return type:

Union[bool, bool]

Returns:

Whether the matrix is orthogonal within the given tolerance.

ffsim.linalg.is_real_symmetric(mat, *, rtol=1e-05, atol=1e-08)[source]¶

Determine if a matrix is real and approximately symmetric.

Parameters:

mat (ndarray) – The matrix.
rtol (float) – Relative numerical tolerance.
atol (float) – Absolute numerical tolerance.

Return type:

Union[bool, bool]

Returns:

Whether the matrix is real and symmetric within the given tolerance.

ffsim.linalg.is_special_orthogonal(mat, *, rtol=1e-05, atol=1e-08)[source]¶

Determine if a matrix is approximately special orthogonal.

Parameters:

mat (ndarray) – The matrix.
rtol (float) – Relative numerical tolerance.
atol (float) – Absolute numerical tolerance.

Return type:

Union[bool, bool]

Returns:

Whether the matrix is special orthogonal within the given tolerance.

ffsim.linalg.is_unitary(mat, *, rtol=1e-05, atol=1e-08)[source]¶

Determine if a matrix is approximately unitary.

Parameters:

mat (ndarray) – The matrix.
rtol (float) – Relative numerical tolerance.
atol (float) – Absolute numerical tolerance.

Return type:

Union[bool, bool]

Returns:

Whether the matrix is unitary within the given tolerance.

ffsim.linalg.lup(mat)[source]¶

Column-pivoted LU decomposition of a matrix.

The decomposition is: :rtype: tuple[ndarray, ndarray, ndarray]

\[A = L U P\]

where L is a lower triangular matrix with unit diagonal elements, U is upper triangular, and P is a permutation matrix.

ffsim.linalg.match_global_phase(a, b)[source]¶

Phase the given arrays so that their phases match at one entry.

Parameters:

a (ndarray) – A Numpy array.
b (ndarray) – Another Numpy array.

Return type:

tuple[ndarray, ndarray]

Returns:

A pair of arrays (a’, b’) that are equal if and only if a == b * exp(i phi) for some real number phi.

ffsim.linalg.modified_cholesky(mat, *, tol=1e-08, max_vecs=None)[source]¶

Modified Cholesky decomposition.

The modified Cholesky decomposition of a square matrix \(M\) has the form

\[M = \sum_{i=1}^N v_i v_i^\dagger\]

where each \(v_i\) is a vector. M must be positive definite. No checking is performed to verify whether M is positive definite. The number of terms \(N\) in the decomposition depends on the allowed error threshold. A larger error threshold may yield a smaller number of terms. Furthermore, the max_vecs parameter specifies an optional upper bound on \(N\). The max_vecs parameter is always respected, so if it is too small, then the error of the decomposition may exceed the specified error threshold.

References

arXiv:1711.02242

Parameters:

mat (ndarray) – The matrix to decompose.
tol (float) – Tolerance for error in the decomposition. The error is defined as the maximum absolute difference between an element of the original tensor and the corresponding element of the reconstructed tensor.
max_vecs (int | None) – The maximum number of vectors to include in the decomposition.

Return type:

ndarray

Returns:

The Cholesky vectors v_i assembled into a 2-dimensional Numpy array whose columns are the vectors.

ffsim.linalg.one_hot(shape, index, *, dtype=<class 'complex'>)[source]¶

Return an array of all zeros except for a one at a specified index.

Parameters:

shape (int | tuple[int, ...]) – The desired shape of the array.
index – The index at which to place a one.

Returns:

The one-hot vector.

ffsim.linalg.reduced_matrix(mat, vecs)[source]¶

Compute reduced matrix within a subspace spanned by some vectors.

Given a linear operator \(A\) and a list of vectors \(\{v_i\}\), return the matrix M where \(M_{ij} = v_i^\dagger A v_j\).

Return type:: ndarray