diff --git a/docs/tutorials/03-double-factorized.ipynb b/docs/tutorials/03-double-factorized.ipynb index 935060978..589df6c62 100644 --- a/docs/tutorials/03-double-factorized.ipynb +++ b/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 }, diff --git a/python/ffsim/__init__.py b/python/ffsim/__init__.py index c03d91001..07d965a20 100644 --- a/python/ffsim/__init__.py +++ b/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, diff --git a/python/ffsim/hamiltonians/__init__.py b/python/ffsim/hamiltonians/__init__.py index 32bef5b91..cd0a9b9c7 100644 --- a/python/ffsim/hamiltonians/__init__.py +++ b/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 diff --git a/python/ffsim/hamiltonians/double_factorized_hamiltonian.py b/python/ffsim/hamiltonians/double_factorized_hamiltonian.py index 0184da262..0a39c044a 100644 --- a/python/ffsim/hamiltonians/double_factorized_hamiltonian.py +++ b/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 diff --git a/tests/hamiltonians/double_factorized_hamiltonian_test.py b/tests/hamiltonians/double_factorized_hamiltonian_test.py index f2458a8e8..2c25167e5 100644 --- a/tests/hamiltonians/double_factorized_hamiltonian_test.py +++ b/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() diff --git a/tests/trotter/qdrift_test.py b/tests/trotter/qdrift_test.py index 9b6810f01..8f28fea05 100644 --- a/tests/trotter/qdrift_test.py +++ b/tests/trotter/qdrift_test.py @@ -267,7 +267,7 @@ def test_simulate_qdrift_double_factorized_h_chain( hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb=norb, nelec=nelec) # perform double factorization - df_hamiltonian = ffsim.double_factorized_hamiltonian( + df_hamiltonian = ffsim.DoubleFactorizedHamiltonian.from_molecular_hamiltonian( mol_hamiltonian, z_representation=z_representation, ) @@ -364,7 +364,7 @@ def test_simulate_qdrift_double_factorized_random( hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb=norb, nelec=nelec) # perform double factorization - df_hamiltonian = ffsim.double_factorized_hamiltonian( + df_hamiltonian = ffsim.DoubleFactorizedHamiltonian.from_molecular_hamiltonian( mol_hamiltonian, optimize=optimize, z_representation=z_representation, diff --git a/tests/trotter/trotter_test.py b/tests/trotter/trotter_test.py index 81d1a5a5d..6be8bc33c 100644 --- a/tests/trotter/trotter_test.py +++ b/tests/trotter/trotter_test.py @@ -48,7 +48,7 @@ def test_simulate_trotter_double_factorized_random( hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb=norb, nelec=nelec) # perform double factorization - df_hamiltonian = ffsim.double_factorized_hamiltonian( + df_hamiltonian = ffsim.DoubleFactorizedHamiltonian.from_molecular_hamiltonian( mol_hamiltonian, z_representation=z_representation )