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, 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=1}^L \sum_{k\ell} Z^{t}_{k\ell} U^{t}_{pk} U^{t}_{qk} U^{t}_{r\ell} U^{t}_{s\ell}\]

Here each \(Z^{(t)}\) is a real symmetric matrix, referred to as a “diagonal Coulomb matrix,” and each \(U^{t}\) is a unitary matrix, referred to as 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. Additionally, one can choose to optimize the coefficients stored in the tensor to achieve a “compressed” factorization. This option is enabled by setting the optimize parameter to True. The optimization attempts to minimize a least-squares objective function quantifying the error in the decomposition. It uses scipy.optimize.minimize, passing both the objective function and its gradient. 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.

  • method (str) – The optimization method. See the documentation of scipy.optimize.minimize for possible values.

  • callback – Callback function for the optimization. See the documentation of scipy.optimize.minimize for usage.

  • options (dict | None) – Options for the optimization. See the documentation of scipy.optimize.minimize for usage.

  • 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 trianglular indices, i.e., pairs \((i, j)\) where \(i \leq j\). Passing a list with lower triangular indices will raise an error. This parameter is only used if optimize is set to True.

  • 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.

Return type:

tuple[ndarray, ndarray]

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 rank of the decomposition and n is the number of orbitals.

Raises:

ValueError – diag_coulomb_indices contains lower triangular indices.

ffsim.linalg.double_factorized_t2(t2_amplitudes, *, tol=1e-08, max_vecs=None)[source]

Double-factorized decomposition of t2 amplitudes.

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

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

Here each \(Z^{(mk)}\) is a real-valued matrix, referred to as a “diagonal Coulomb matrix,” and each \(U^{(mk)}\) is a unitary matrix, referred to as 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.

Note: Currently, only real-valued t2 amplitudes are supported.

Parameters:

  • t2_amplitudes (ndarray) – The t2 amplitudes tensor.

  • 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 t2 amplitudes tensor. This argument overrides tol.

Return type:

tuple[ndarray, ndarray]

Returns:

  • The diagonal Coulomb matrices, as a Numpy array of shape (n_vecs, 2, norb, norb). The last two axes index the rows and columns of the matrices. The first axis indexes the eigenvectors of the decomposition and the second axis exists because each eigenvector gives rise to 2 terms in the decomposition.

  • The orbital rotations, as a Numpy array of shape (n_vecs, 2, norb, norb). The last two axes index the rows and columns of the orbital rotations. The first axis indexes the eigenvectors of the decomposition and the second axis exists because each eigenvector gives rise to 2 terms in the decomposition.

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

Double-factorized decomposition of alpha-beta t2 amplitudes.

Decompose alpha-beta t2 amplitudes into diagonal Coulomb matrices with orbital rotations. This function returns two arrays:

  • diagonal_coulomb_mats, with shape (n_vecs, 4, 3, norb, norb).

  • orbital_rotations, with shape (n_vecs, 4, 2, norb, norb).

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

The original t2 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_vecs = diag_coulomb_mats.shape[0]
    expanded_diag_coulomb_mats = np.zeros((n_vecs, 4, 2 * norb, 2 * norb))
    expanded_orbital_rotations = np.zeros(
        (n_vecs, 4, 2 * norb, 2 * norb), dtype=complex
    )
    for m, k in itertools.product(range(n_vecs), range(4)):
        (mat_aa, mat_ab, mat_bb) = diag_coulomb_mats[m, k]
        expanded_diag_coulomb_mats[m, k] = np.block(
            [[mat_aa, mat_ab], [mat_ab.T, mat_bb]]
        )
        orbital_rotation_a, orbital_rotation_b = orbital_rotations[m, k]
        expanded_orbital_rotations[m, k] = scipy.linalg.block_diag(
            orbital_rotation_a, orbital_rotation_b
        )
    return (
        2j
        * contract(
            "mkpq,mkap,mkip,mkbq,mkjq->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 t2 amplitudes are supported.

Parameters:

  • t2_amplitudes (ndarray) – The t2 amplitudes tensor.

  • 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 t2 amplitudes tensor. This argument overrides tol.

Return type:

tuple[ndarray, ndarray]

Returns:

  • The diagonal Coulomb matrices, as a Numpy array of shape (n_vecs, 4, 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 singular vectors of the decomposition and the second axis exists because each singular vector gives rise to 4 terms in the decomposition.

  • The orbital rotations, as a Numpy array of shape (n_vecs, 4, 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 singular vectors of the decomposition and the second axis exists because each singular vector gives rise to 4 terms in 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

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