add random_double_factorized_hamiltonian (#274)

python/ffsim/random/__init__.py

@@ -12,6 +12,7 @@
 
 from ffsim.random.random import (
     random_antihermitian,
+    random_double_factorized_hamiltonian,
     random_hermitian,
     random_molecular_hamiltonian,
     random_orthogonal,
@@ -30,6 +31,7 @@
 
 __all__ = [
     "random_antihermitian",
+    "random_double_factorized_hamiltonian",
     "random_hermitian",
     "random_molecular_hamiltonian",
     "random_orthogonal",

python/ffsim/random/random.py

@@ -275,7 +275,7 @@ def random_molecular_hamiltonian(
     """Sample a random molecular Hamiltonian.
 
     Args:
-        norb: The number of orbitals.
+        norb: The number of spatial orbitals.
         seed: A seed to initialize the pseudorandom number generator.
             Should be a valid input to ``np.random.default_rng``.
         dtype: The data type to use for the one- and two-body tensors. The constant
@@ -353,7 +353,7 @@ def random_ucj_op_spin_balanced(
     """Sample a random spin-balanced unitary cluster Jastrow (UCJ) operator.
 
     Args:
-        norb: The number of orbitals.
+        norb: The number of spatial orbitals.
         n_reps: The number of ansatz repetitions.
         with_final_orbital_rotation: Whether to include a final orbital rotation
             in the operator.
@@ -474,3 +474,37 @@ def random_ucj_op_spinless(
         orbital_rotations=orbital_rotations,
         final_orbital_rotation=final_orbital_rotation,
     )
+
+
+def random_double_factorized_hamiltonian(
+    norb: int, *, rank: int | None = None, z_representation: bool = False, seed=None
+) -> hamiltonians.DoubleFactorizedHamiltonian:
+    """Sample a random double-factorized Hamiltonian.
+
+    Args:
+        norb: The number of spatial orbitals.
+        rank: The desired number of terms in the two-body part of the Hamiltonian.
+            If not specified, it will be set to ``norb * (norb + 1) // 2``.
+        z_representation: Whether to return a Hamiltonian in the "Z" representation.
+        seed: A seed to initialize the pseudorandom number generator.
+            Should be a valid input to ``np.random.default_rng``.
+
+    Returns:
+        The sampled double-factorized Hamiltonian.
+    """
+    if rank is None:
+        rank = norb * (norb + 1) // 2
+    rng = np.random.default_rng(seed)
+    one_body_tensor = random_hermitian(norb, seed=rng)
+    diag_coulomb_mats = np.stack(
+        [random_real_symmetric_matrix(norb, seed=rng) for _ in range(rank)]
+    )
+    orbital_rotations = np.stack([random_unitary(norb, seed=rng) for _ in range(rank)])
+    constant = rng.standard_normal()
+    return hamiltonians.DoubleFactorizedHamiltonian(
+        one_body_tensor=one_body_tensor,
+        diag_coulomb_mats=diag_coulomb_mats,
+        orbital_rotations=orbital_rotations,
+        constant=constant,
+        z_representation=z_representation,
+    )

tests/python/trotter/double_factorized_test.py

@@ -22,11 +22,12 @@
 @pytest.mark.parametrize(
     "norb, nelec, time, n_steps, order, z_representation, target_fidelity",
     [
-        (3, (1, 1), 0.1, 10, 0, False, 0.99),
-        (3, (1, 1), 0.1, 3, 1, True, 0.99),
-        (4, (2, 1), 0.1, 3, 2, False, 0.99),
-        (4, (1, 2), 0.1, 3, 2, True, 0.99),
-        (4, (2, 2), 0.1, 8, 1, False, 0.99),
+        (3, (1, 1), 0.1, 10, 0, False, 0.999),
+        (3, (1, 1), 0.1, 2, 1, True, 0.999),
+        (4, (2, 1), 0.1, 1, 2, False, 0.999),
+        (4, (1, 2), 0.1, 3, 2, True, 0.999),
+        (4, (2, 2), 0.1, 4, 1, False, 0.999),
+        (5, (3, 2), 0.1, 5, 1, True, 0.999),
     ],
 )
 def test_random(
@@ -38,19 +39,15 @@ def test_random(
     z_representation: bool,
     target_fidelity: float,
 ):
+    """Test random Hamiltonian."""
     rng = np.random.default_rng(2488)
 
     # generate random Hamiltonian
     dim = ffsim.dim(norb, nelec)
-    one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng)
-    two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng, dtype=float)
-    mol_hamiltonian = ffsim.MolecularHamiltonian(one_body_tensor, two_body_tensor)
-    hamiltonian = ffsim.linear_operator(mol_hamiltonian, norb=norb, nelec=nelec)
-
-    # perform double factorization
-    df_hamiltonian = ffsim.DoubleFactorizedHamiltonian.from_molecular_hamiltonian(
-        mol_hamiltonian, z_representation=z_representation
+    hamiltonian = ffsim.random.random_double_factorized_hamiltonian(
+        norb, rank=norb, z_representation=z_representation, seed=rng
     )
+    linop = ffsim.linear_operator(hamiltonian, norb=norb, nelec=nelec)
 
     # generate initial state
     dim = ffsim.dim(norb, nelec)
@@ -59,18 +56,18 @@ def test_random(
 
     # compute exact state
     exact_state = scipy.sparse.linalg.expm_multiply(
-        -1j * time * hamiltonian,
+        -1j * time * linop,
         initial_state,
-        traceA=ffsim.trace(mol_hamiltonian, norb=norb, nelec=nelec),
+        traceA=1.0,
     )
 
     # make sure time is not too small
-    assert abs(np.vdot(exact_state, initial_state)) < 0.98
+    assert abs(np.vdot(exact_state, initial_state)) < 0.95
 
     # simulate
     final_state = ffsim.simulate_trotter_double_factorized(
         initial_state,
-        df_hamiltonian,
+        hamiltonian,
         time,
         norb=norb,
         nelec=nelec,