Source code for ffsim.hamiltonians.double_factorized_hamiltonian

# (C) Copyright IBM 2023.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

from __future__ import annotations

import dataclasses

import numpy as np
import scipy.linalg
from scipy.sparse.linalg import LinearOperator

from ffsim.contract.diag_coulomb import diag_coulomb_linop
from ffsim.contract.num_op_sum import num_op_sum_linop
from ffsim.hamiltonians.molecular_hamiltonian import MolecularHamiltonian
from ffsim.linalg import double_factorized
from ffsim.operators import FermionOperator
from ffsim.protocols import fermion_operator
from ffsim.states import dim


[docs] @dataclasses.dataclass(frozen=True) class DoubleFactorizedHamiltonian: r"""A Hamiltonian in the double-factorized representation. The double-factorized form of the molecular Hamiltonian is .. math:: H = \sum_{\sigma, pq} \kappa_{pq} a^\dagger_{\sigma, p} a_{\sigma, q} + \frac12 \sum_t \sum_{\sigma\tau, ij} Z^{(t)}_{ij} n^{(t)}_{\sigma, i} n^{(t)}_{\tau, j} + \text{constant}'. where .. math:: n^{(t)}_{\sigma, i} = \sum_{pq} U^{(t)}_{pi} a^\dagger_{\sigma, p} a_{\sigma, q} U^{(t)}_{qi}. Here each :math:`U^{(t)}` is a unitary matrix and each :math:`Z^{(t)}` is a real symmetric matrix. **"Z" representation** The "Z" representation of the double factorization is an alternative representation that sometimes yields simpler quantum circuits. Under the Jordan-Wigner transformation, the number operators take the form .. math:: n^{(t)}_{\sigma, i} = \frac{(1 - z^{(t)}_{\sigma, i})}{2} where :math:`z^{(t)}_{\sigma, i}` is the Pauli Z operator in the rotated basis. The "Z" representation is obtained by rewriting the two-body part in terms of these Pauli Z operators and updating the one-body term as appropriate: .. math:: H = \sum_{\sigma, pq} \kappa'_{pq} a^\dagger_{\sigma, p} a_{\sigma, q} + \frac18 \sum_t \sum_{\sigma\tau, ij}^* Z^{(t)}_{ij} z^{(t)}_{\sigma, i} z^{(t)}_{\tau, j} + \text{constant}'' where the asterisk denotes summation over indices :math:`\sigma\tau, ij` where :math:`\sigma \neq \tau` or :math:`i \neq j`. References: - `Low rank representations for quantum simulation of electronic structure`_ - `Quantum Filter Diagonalization with Double-Factorized Hamiltonians`_ Attributes: one_body_tensor (np.ndarray): The one-body tensor :math:`\kappa`. diag_coulomb_mats (np.ndarray): The diagonal Coulomb matrices. orbital_rotations (np.ndarray): The orbital rotations. constant (float): The constant. z_representation (bool): Whether the Hamiltonian is in the "Z" representation rather than the "number" representation. .. _Low rank representations for quantum simulation of electronic structure: https://arxiv.org/abs/1808.02625 .. _Quantum Filter Diagonalization with Double-Factorized Hamiltonians: https://arxiv.org/abs/2104.08957 .. _scipy.optimize.minimize: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html """ one_body_tensor: np.ndarray diag_coulomb_mats: np.ndarray orbital_rotations: np.ndarray constant: float = 0.0 z_representation: bool = False @property def norb(self) -> int: """The number of spatial orbitals.""" return self.one_body_tensor.shape[0]
[docs] def to_z_representation(self) -> "DoubleFactorizedHamiltonian": """Return the Hamiltonian in the "Z" representation.""" if self.z_representation: return self one_body_correction, constant_correction = _df_z_representation( self.diag_coulomb_mats, self.orbital_rotations ) return DoubleFactorizedHamiltonian( one_body_tensor=self.one_body_tensor + one_body_correction, diag_coulomb_mats=self.diag_coulomb_mats, orbital_rotations=self.orbital_rotations, constant=self.constant + constant_correction, z_representation=True, )
[docs] def to_number_representation(self) -> "DoubleFactorizedHamiltonian": """Return the Hamiltonian in the "number" representation.""" if not self.z_representation: return self one_body_correction, constant_correction = _df_z_representation( self.diag_coulomb_mats, self.orbital_rotations ) return DoubleFactorizedHamiltonian( one_body_tensor=self.one_body_tensor - one_body_correction, diag_coulomb_mats=self.diag_coulomb_mats, orbital_rotations=self.orbital_rotations, constant=self.constant - constant_correction, z_representation=False, )
[docs] @staticmethod def from_molecular_hamiltonian( hamiltonian: MolecularHamiltonian, *, z_representation: bool = False, tol: float = 1e-8, max_vecs: int | None = None, optimize: bool = False, method: str = "L-BFGS-B", callback=None, options: dict | None = None, diag_coulomb_indices: list[tuple[int, int]] | None = None, cholesky: bool = True, ) -> DoubleFactorizedHamiltonian: r"""Initialize a DoubleFactorizedHamiltonian from a MolecularHamiltonian. This function takes as input a :class:`MolecularHamiltonian`, which stores a one-body tensor, two-body tensor, and constant. It performs a double-factorized decomposition of the two-body tensor and computes a new one-body tensor and constant, and returns a :class:`DoubleFactorizedHamiltonian` storing the results. See :class:`DoubleFactorizedHamiltonian` for a description of the `z_representation` argument. See :func:`ffsim.linalg.double_factorized` for a description of the rest of the arguments. Args: hamiltonian: The Hamiltonian whose double-factorized representation to compute. z_representation: Whether to use the "Z" representation of the decomposition. tol: 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: An optional limit on the number of terms to keep in the decomposition of the two-body tensor. This argument overrides ``tol``. optimize: Whether to optimize the tensors returned by the decomposition. method: 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: Options for the optimization. See the documentation of `scipy.optimize.minimize`_ for usage. diag_coulomb_indices: 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 :math:`(i, j)` where :math:`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: 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 double-factorized Hamiltonian. """ one_body_tensor = hamiltonian.one_body_tensor - 0.5 * np.einsum( "prqr", hamiltonian.two_body_tensor ) diag_coulomb_mats, orbital_rotations = double_factorized( hamiltonian.two_body_tensor, tol=tol, max_vecs=max_vecs, optimize=optimize, method=method, callback=callback, options=options, diag_coulomb_indices=diag_coulomb_indices, cholesky=cholesky, ) df_hamiltonian = DoubleFactorizedHamiltonian( one_body_tensor=one_body_tensor, diag_coulomb_mats=diag_coulomb_mats, orbital_rotations=orbital_rotations, constant=hamiltonian.constant, ) if z_representation: df_hamiltonian = df_hamiltonian.to_z_representation() return df_hamiltonian
[docs] def to_molecular_hamiltonian(self): """Convert the DoubleFactorizedHamiltonian to a MolecularHamiltonian.""" df_hamiltonian = self.to_number_representation() two_body_tensor = np.einsum( "kij,kpi,kqi,krj,ksj->pqrs", df_hamiltonian.diag_coulomb_mats, df_hamiltonian.orbital_rotations, df_hamiltonian.orbital_rotations.conj(), df_hamiltonian.orbital_rotations, df_hamiltonian.orbital_rotations.conj(), ) one_body_tensor = df_hamiltonian.one_body_tensor + 0.5 * np.einsum( "prqr", two_body_tensor ) return MolecularHamiltonian( one_body_tensor=one_body_tensor, two_body_tensor=two_body_tensor, constant=df_hamiltonian.constant, )
def _linear_operator_(self, norb: int, nelec: tuple[int, int]) -> LinearOperator: """Return a SciPy LinearOperator representing the object.""" dim_ = dim(norb, nelec) eigs, vecs = scipy.linalg.eigh(self.one_body_tensor) num_linop = num_op_sum_linop(eigs, norb, nelec, orbital_rotation=vecs) diag_coulomb_linops = [ diag_coulomb_linop( diag_coulomb_mat, norb, nelec, orbital_rotation=orbital_rotation, z_representation=self.z_representation, ) for diag_coulomb_mat, orbital_rotation in zip( self.diag_coulomb_mats, self.orbital_rotations ) ] def matvec(vec: np.ndarray): vec = vec.astype(complex, copy=False) result = self.constant * vec result += num_linop @ vec for linop in diag_coulomb_linops: result += linop @ vec return result return LinearOperator( shape=(dim_, dim_), matvec=matvec, rmatvec=matvec, dtype=complex ) def _diag_(self, norb: int, nelec: tuple[int, int]) -> np.ndarray: """Return the diagonal entries of the Hamiltonian.""" return self.to_molecular_hamiltonian()._diag_(norb, nelec) def _fermion_operator_(self) -> FermionOperator: """Return a FermionOperator representing the object.""" return fermion_operator(self.to_molecular_hamiltonian())
def _df_z_representation( diag_coulomb_mats: np.ndarray, orbital_rotations: np.ndarray ) -> tuple[np.ndarray, float]: one_body_correction = 0.5 * ( np.einsum( "tij,tpi,tqi->pq", diag_coulomb_mats, orbital_rotations, orbital_rotations.conj(), ) + np.einsum( "tij,tpj,tqj->pq", diag_coulomb_mats, orbital_rotations, orbital_rotations.conj(), ) ) constant_correction = 0.25 * np.einsum("ijj->", diag_coulomb_mats) - 0.5 * np.sum( diag_coulomb_mats ) return one_body_correction, constant_correction