ffsim.linalg.double_factorized_t2_alpha_beta

ffsim.linalg.double_factorized_t2_alpha_beta(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 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 (
        1j
        * 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.

  • 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 (tuple[list[tuple[int, int]] | None, list[tuple[int, int]] | None, list[tuple[int, int]] | None] | None) – Allowed indices for nonzero values of the diagonal Coulomb matrices. If specified, should be a tuple of 3 lists, for alpha-alpha, alpha-beta, and beta-beta matrices, in that order. For the alpha-alpha and beta-beta matrices, each index must be upper triangular, that is, of the form \((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:

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