ffsim.linalg.double_factorized_t2

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)

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:

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.