change double_factorized_hamiltonian to class static method

docs/tutorials/03-double-factorized.ipynb

@@ -63,14 +63,16 @@
     "mol_hamiltonian = mol_data.hamiltonian\n",
     "\n",
     "# Get the Hamiltonian in the double-factorized representation\n",
-    "df_hamiltonian = ffsim.double_factorized_hamiltonian(mol_hamiltonian)"
+    "df_hamiltonian = ffsim.DoubleFactorizedHamiltonian.from_molecular_hamiltonian(\n",
+    "    mol_hamiltonian\n",
+    ")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The function `double_factorized_hamiltonian` returns an object of type `DoubleFactorizedHamiltonian` which is just a dataclass that stores the updated one-body-tensor, diagonal Coulomb matrices, and orbital rotations. In the cell below, we print out the tensors describing the original and double-factorized representations."
+    "Here, `mol_hamiltonian` is an instance of `MolecularHamiltonian`, a dataclass that stores the one- and two-body tensors, and `df_hamiltonian` is an instance of `DoubleFactorizedHamiltonian`, a dataclass that stores the updated one-body-tensor, diagonal Coulomb matrices, and orbital rotations. In the cell below, we print out the tensors describing the original and double-factorized representations."
    ]
   },
   {
@@ -404,7 +406,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.11"
+   "version": "3.10.13"
   },
   "orig_nbformat": 4
  },

python/ffsim/__init__.py

@@ -44,11 +44,7 @@
     apply_orbital_rotation,
     apply_tunneling_interaction,
 )
-from ffsim.hamiltonians import (
-    DoubleFactorizedHamiltonian,
-    MolecularHamiltonian,
-    double_factorized_hamiltonian,
-)
+from ffsim.hamiltonians import DoubleFactorizedHamiltonian, MolecularHamiltonian
 from ffsim.molecular_data import MolecularData
 from ffsim.protocols import (
     SupportsApproximateEquality,

python/ffsim/hamiltonians/__init__.py

@@ -10,8 +10,5 @@
 
 """Classes for representing Hamiltonians."""
 
-from ffsim.hamiltonians.double_factorized_hamiltonian import (
-    DoubleFactorizedHamiltonian,
-    double_factorized_hamiltonian,
-)
+from ffsim.hamiltonians.double_factorized_hamiltonian import DoubleFactorizedHamiltonian
 from ffsim.hamiltonians.molecular_hamiltonian import MolecularHamiltonian

python/ffsim/hamiltonians/double_factorized_hamiltonian.py

@@ -118,6 +118,163 @@ def to_number_representation(self) -> "DoubleFactorizedHamiltonian":
             z_representation=False,
         )
 
+    @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",
+        options: dict | None = None,
+        diag_coulomb_mask: np.ndarray | None = None,
+        cholesky: bool = True,
+    ) -> DoubleFactorizedHamiltonian:
+        r"""Double-factorized decomposition of a molecular Hamiltonian.
+
+        The double-factorized decomposition acts on a Hamiltonian of the form
+
+        .. math::
+
+            H = \sum_{pq, \sigma} h_{pq} a^\dagger_{p, \sigma} a_{q, \sigma}
+                + \frac12 \sum_{pqrs, \sigma \tau} h_{pqrs}
+                a^\dagger_{p, \sigma} a^\dagger_{r, \tau} a_{s, \tau} a_{q, \sigma}
+                + \text{constant}.
+
+        The Hamiltonian is decomposed into the double-factorized form
+
+        .. math::
+
+            H = \sum_{pq, \sigma} \kappa_{pq} a^\dagger_{p, \sigma} a_{q, \sigma}
+            + \frac12 \sum_t \sum_{ij, \sigma\tau}
+            Z^{(t)}_{ij} n^{(t)}_{i, \sigma} n^{(t)}_{j, \tau}
+            + \text{constant}'.
+
+        where
+
+        .. math::
+
+            n^{(t)}_{i, \sigma} = \sum_{pq} U^{(t)}_{pi}
+            a^\dagger_{p, \sigma} a^\dagger_{q, \sigma} U^{(t)}_{qi}.
+
+        Here :math:`U^{(t)}_{ij}` and :math:`Z^{(t)}_{ij}` are tensors that are output
+        by the decomposition, and :math:`\kappa_{pq}` is an updated one-body tensor.
+        Each matrix :math:`U^{(t)}` is guaranteed to be unitary so that the
+        :math:`n^{(t)}_{i, \sigma}` are number operators in a rotated basis, and
+        each :math:`Z^{(t)}` is a real symmetric matrix.
+        The number of terms :math:`t` in the decomposition depends on the allowed
+        error threshold. A larger error threshold leads to a smaller number of terms.
+        Furthermore, the `max_rank` parameter specifies an optional upper bound
+        on :math:`t`.
+
+        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 low rank 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_mask` parameter, which is
+        an :math:`N \times N` matrix of boolean values where :math:`N` is the number
+        of orbitals. The nonzero elements of this matrix indicate where the diagonal
+        Coulomb matrices are allowed to be nonzero. Only the upper triangular part of
+        the matrix is used because the diagonal Coulomb matrices are symmetric.
+
+        **"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)}_{i, \sigma} = \frac{(1 - z^{(t)}_{i, \sigma})}{2}
+
+        where :math:`z^{(t)}_{i, \sigma}` 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_{pq, \sigma} \kappa'_{pq} a^\dagger_{p, \sigma} a_{q, \sigma}
+            + \frac18 \sum_t \sum_{ij, \sigma\tau}^*
+            Z^{(t)}_{ij} z^{(t)}_{i, \sigma} z^{(t)}_{j, \tau}
+            + \text{constant}''
+
+        where the asterisk denotes summation over indices :math:`ij, \sigma\tau`
+        where :math:`i \neq j` or :math:`\sigma \neq \tau`.
+
+        Note: Currently, only real-valued two-body tensors are supported.
+
+        Args:
+            one_body_tensor: The one-body tensor of the Hamiltonian.
+            two_body_tensor: The two-body tensor of the Hamiltonian.
+            z_representation: Whether to use the "Z" representation of the
+                low rank 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_mask: Diagonal Coulomb matrix mask to use in the optimization.
+                This is a matrix of boolean values where the nonzero elements indicate
+                where the diagonal Coulomb matrices returned by optimization are allowed
+                to be nonzero. This parameter is only used if `optimize` is set to
+                True, and only the upper triangular part of the matrix is used.
+            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.
+
+        References:
+            - `arXiv:1808.02625`_
+            - `arXiv:2104.08957`_
+
+        .. _arXiv:1808.02625: https://arxiv.org/abs/1808.02625
+        .. _arXiv:2104.08957: https://arxiv.org/abs/2104.08957
+        .. _scipy.optimize.minimize: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
+        """
+        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,
+            options=options,
+            diag_coulomb_mask=diag_coulomb_mask,
+            cholesky=cholesky,
+        )
+        df_hamiltonian = DoubleFactorizedHamiltonian(
+            one_body_tensor=one_body_tensor,
+            diag_coulomb_mats=diag_coulomb_mats,
+            orbital_rotations=orbital_rotations,
+        )
+
+        if z_representation:
+            df_hamiltonian = df_hamiltonian.to_z_representation()
+
+        return df_hamiltonian
+
 
 def _df_z_representation(
     diag_coulomb_mats: np.ndarray, orbital_rotations: np.ndarray
@@ -140,159 +297,3 @@ def _df_z_representation(
         diag_coulomb_mats
     )
     return one_body_correction, constant_correction
-
-
-def double_factorized_hamiltonian(
-    hamiltonian: MolecularHamiltonian,
-    *,
-    z_representation: bool = False,
-    tol: float = 1e-8,
-    max_vecs: int | None = None,
-    optimize: bool = False,
-    method: str = "L-BFGS-B",
-    options: dict | None = None,
-    diag_coulomb_mask: np.ndarray | None = None,
-    cholesky: bool = True,
-) -> DoubleFactorizedHamiltonian:
-    r"""Double-factorized decomposition of a molecular Hamiltonian.
-
-    The double-factorized decomposition acts on a Hamiltonian of the form
-
-    .. math::
-
-        H = \sum_{pq, \sigma} h_{pq} a^\dagger_{p, \sigma} a_{q, \sigma}
-            + \frac12 \sum_{pqrs, \sigma \tau} h_{pqrs}
-            a^\dagger_{p, \sigma} a^\dagger_{r, \tau} a_{s, \tau} a_{q, \sigma}
-            + \text{constant}.
-
-    The Hamiltonian is decomposed into the double-factorized form
-
-    .. math::
-
-        H = \sum_{pq, \sigma} \kappa_{pq} a^\dagger_{p, \sigma} a_{q, \sigma}
-        + \frac12 \sum_t \sum_{ij, \sigma\tau}
-        Z^{(t)}_{ij} n^{(t)}_{i, \sigma} n^{(t)}_{j, \tau}
-        + \text{constant}'.
-
-    where
-
-    .. math::
-
-        n^{(t)}_{i, \sigma} = \sum_{pq} U^{(t)}_{pi}
-        a^\dagger_{p, \sigma} a^\dagger_{q, \sigma} U^{(t)}_{qi}.
-
-    Here :math:`U^{(t)}_{ij}` and :math:`Z^{(t)}_{ij}` are tensors that are output by
-    the decomposition, and :math:`\kappa_{pq}` is an updated one-body tensor.
-    Each matrix :math:`U^{(t)}` is guaranteed to be unitary so that the
-    :math:`n^{(t)}_{i, \sigma}` are number operators in a rotated basis, and
-    each :math:`Z^{(t)}` is a real symmetric matrix.
-    The number of terms :math:`t` in the decomposition depends on the allowed
-    error threshold. A larger error threshold leads to a smaller number of terms.
-    Furthermore, the `max_rank` parameter specifies an optional upper bound
-    on :math:`t`.
-
-    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 low rank decomposition.
-    It uses `scipy.optimize.minimize`, passing both the objective function
-    and its gradient. The core tensors 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_mask` parameter, which is an :math:`N \times N` matrix of
-    boolean values where :math:`N` is the number of orbitals. The nonzero elements of
-    this matrix indicate where the core tensors are allowed to be nonzero. Only the
-    upper triangular part of the matrix is used because the core tensors are symmetric.
-
-    **"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)}_{i, \sigma} = \frac{(1 - z^{(t)}_{i, \sigma})}{2}
-
-    where :math:`z^{(t)}_{i, \sigma}` 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_{pq, \sigma} \kappa'_{pq} a^\dagger_{p, \sigma} a_{q, \sigma}
-        + \frac18 \sum_t \sum_{ij, \sigma\tau}^*
-        Z^{(t)}_{ij} z^{(t)}_{i, \sigma} z^{(t)}_{j, \tau}
-        + \text{constant}''
-
-    where the asterisk denotes summation over indices :math:`ij, \sigma\tau`
-    where :math:`i \neq j` or :math:`\sigma \neq \tau`.
-
-    Note: Currently, only real-valued two-body tensors are supported.
-
-    Args:
-        one_body_tensor: The one-body tensor of the Hamiltonian.
-        two_body_tensor: The two-body tensor of the Hamiltonian.
-        z_representation: Whether to use the "Z" representation of the
-            low rank 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_mask: Diagonal Coulomb matrix mask to use in the optimization.
-            This is a matrix of boolean values where the nonzero elements indicate where
-            the diagonal coulomb matrices returned by optimization are allowed to be
-            nonzero. This parameter is only used if `optimize` is set to `True`, and
-            only the upper triangular part of the matrix is used.
-        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.
-
-    References:
-        - `arXiv:1808.02625`_
-        - `arXiv:2104.08957`_
-
-    .. _arXiv:1808.02625: https://arxiv.org/abs/1808.02625
-    .. _arXiv:2104.08957: https://arxiv.org/abs/2104.08957
-    .. _scipy.optimize.minimize: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
-    """
-    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,
-        options=options,
-        diag_coulomb_mask=diag_coulomb_mask,
-        cholesky=cholesky,
-    )
-    df_hamiltonian = DoubleFactorizedHamiltonian(
-        one_body_tensor=one_body_tensor,
-        diag_coulomb_mats=diag_coulomb_mats,
-        orbital_rotations=orbital_rotations,
-    )
-
-    if z_representation:
-        df_hamiltonian = df_hamiltonian.to_z_representation()
-
-    return df_hamiltonian

tests/hamiltonians/double_factorized_hamiltonian_test.py

@@ -37,7 +37,7 @@ def test_double_factorized_hamiltonian(z_representation: bool):
     )
 
     # perform double factorization
-    df_hamiltonian = ffsim.double_factorized_hamiltonian(
+    df_hamiltonian = ffsim.DoubleFactorizedHamiltonian.from_molecular_hamiltonian(
         mol_hamiltonian,
         z_representation=z_representation,
     )
@@ -86,7 +86,7 @@ def test_z_representation_round_trip():
     one_body_tensor = ffsim.random.random_hermitian(norb, seed=2474)
     two_body_tensor = ffsim.random.random_two_body_tensor_real(norb, seed=7054)
 
-    df_hamiltonian = ffsim.double_factorized_hamiltonian(
+    df_hamiltonian = ffsim.DoubleFactorizedHamiltonian.from_molecular_hamiltonian(
         ffsim.MolecularHamiltonian(one_body_tensor, two_body_tensor)
     )
     df_hamiltonian_num = df_hamiltonian.to_z_representation().to_number_representation()