Preprocessing examples
Classical preprocessing that shapes an electronic-structure Hamiltonian before it ever reaches a quantum circuit. These examples double-factorize the two-electron integrals and trace the result through to block-encoding cost.
See also
Preprocessing: molecule to qubit Hamiltonian — the chemistry bridges and double factorization these examples build on.
Double factorization of the two-electron integrals
Explicit (X-DF) and compressed (C-DF) double factorization of the two-electron integrals, following Cohn, Motta, and Parrish. Generates restricted Hartree-Fock molecular-orbital integrals for H2O/STO-3G with PySCF, then factorizes the ERI tensor both the exact and the least-squares-optimized way.
"""Explicit (X-DF) and compressed (C-DF) double factorization of the two-electron
integrals, following Cohn, Motta, and Parrish, PRX Quantum 2, 040352 (2021).
Generates restricted Hartree-Fock molecular-orbital integrals for H2O/STO-3G with
PySCF, then double-factorizes the ERI tensor two ways:
* X-DF: the exact factorization (rank-one cores) via pivoted Cholesky of
the ERI supermatrix (the default), truncated by residual pivot.
* C-DF: a least-squares-optimized factorization that reaches comparable
accuracy with far fewer leaves.
Heavy linear algebra uses the NVIDIA math libraries (cuSOLVER/cuBLAS via CuPy)
when a GPU is available; otherwise it falls back to NumPy/SciPy.
"""
from __future__ import annotations
import numpy as np
import cudaq_algorithms as algorithms
df = algorithms.double_factorization
def water_mo_eri():
from pyscf import ao2mo, gto, scf
mol = gto.M(atom="O 0 0 0; H 0 0 0.957; H 0 0.926 -0.24",
basis="sto-3g",
verbose=0)
mf = scf.RHF(mol).run()
n = mf.mo_coeff.shape[1]
return np.asarray(ao2mo.restore("s1", ao2mo.kernel(mol, mf.mo_coeff), n))
def main():
backend = "auto"
_, backend_name = df.resolve_backend(backend)
eri = water_mo_eri()
norm = np.linalg.norm(eri)
print("Double factorization of H2O/STO-3G two-electron integrals")
print(f" backend: {backend_name}")
print(f" orbitals: {eri.shape[0]} ||eri||_F: {norm:.6f}")
full = df.explicit_double_factorization(eri,
threshold=0.0,
backend=backend)
print("\nX-DF (explicit):")
print(f" full-rank leaves: {full.num_leaves} "
f"reconstruction error: {df.factorization_error(eri, full):.3e}")
for threshold in (1.0e-2, 1.0e-3, 1.0e-4):
truncated = df.explicit_double_factorization(eri,
threshold=threshold,
backend=backend)
rel = df.factorization_error(eri, truncated) / norm
print(
f" threshold {threshold:.0e}: leaves={truncated.num_leaves:2d} "
f"rel error={rel:.3e}")
print("\nC-DF (compressed) vs X-DF at equal leaf count:")
for num_leaves in (2, 4):
explicit = df.explicit_double_factorization(eri,
threshold=0.0,
max_num_leaves=num_leaves,
backend=backend)
compressed = df.compressed_double_factorization(eri,
num_leaves=num_leaves,
max_iterations=600,
backend=backend)
x_rel = df.factorization_error(eri, explicit) / norm
c_rel = df.factorization_error(eri, compressed) / norm
print(f" leaves={num_leaves:2d}: X-DF rel error={x_rel:.3e} "
f"C-DF rel error={c_rel:.3e}")
assert c_rel <= x_rel + 1.0e-6
if __name__ == "__main__":
main()
From double-factorized integrals to LCU/QSVT cost
Takes H2/STO-3G integrals, double-factorizes the two-electron tensor, and for
each leaf count reconstructs the truncated tensor and bridges it to a qubit
PauliLCU block encoding. The result is an end-to-end view of how DF
compression drives the downstream LCU and QSVT cost.
"""From double-factorized integrals to LCU/QSVT cost, end to end.
Takes H2/STO-3G molecular integrals (hardcoded, so the example needs no
PySCF), double-factorizes the two-electron tensor, and — for each leaf
count — reconstructs the truncated tensor, bridges it to a qubit
Hamiltonian (``chemistry.qubit_hamiltonian``: spin expansion +
Jordan-Wigner), and reports what the truncation buys on the quantum side:
the PauliLCU normalization ``alpha`` (which sets QSVT polynomial degree
via ``d ~ alpha * t`` for time evolution), the Pauli term count (SELECT
cost), and the exact ground-state shift it costs.
Uses ``fermion.jordan_wigner`` via ``chemistry.qubit_hamiltonian``.
"""
from __future__ import annotations
import numpy as np
import cudaq_algorithms as algorithms
from cudaq_algorithms import PauliLCU, chemistry
df = algorithms.double_factorization
# H2 / STO-3G at R = 0.7414 A: MO-basis core Hamiltonian and chemist-
# notation (pq|rs) two-electron integrals; FCI total energy -1.137270 Ha.
ONE_BODY = np.array([[-1.25246357, 0.0], [0.0, -0.47594871]])
ERI = np.zeros((2, 2, 2, 2))
ERI[0, 0, 0, 0] = 0.67449876
ERI[1, 1, 1, 1] = 0.69716349
ERI[0, 0, 1, 1] = ERI[1, 1, 0, 0] = 0.66347258
ERI[0, 1, 0, 1] = ERI[1, 0, 1, 0] = 0.18128881
ERI[0, 1, 1, 0] = ERI[1, 0, 0, 1] = 0.18128881
E_NUCLEAR = 0.71375697
def ground_energy(spin_op) -> float:
return float(np.min(np.linalg.eigvalsh(spin_op.to_matrix())))
def main() -> None:
exact_h = chemistry.qubit_hamiltonian(ONE_BODY,
ERI,
scalar_offset=E_NUCLEAR)
exact_energy = ground_energy(exact_h)
exact_lcu = PauliLCU(exact_h)
print(f"H2/STO-3G on {exact_h.qubit_count} qubits: "
f"exact ground state {exact_energy:+.6f} Ha, "
f"alpha = {exact_lcu.alpha:.4f}, "
f"{exact_lcu.num_terms} LCU terms")
full = df.explicit_double_factorization(ERI, threshold=0.0)
max_leaves = full.num_leaves
print(f"\nX-DF of the ERI tensor: {max_leaves} leaves at full rank")
header = (f"{'leaves':>6} {'tensor error':>13} {'alpha':>8} "
f"{'terms':>6} {'dE_ground (Ha)':>15}")
print(header)
print("-" * len(header))
for num_leaves in range(1, max_leaves + 1):
truncated = df.explicit_double_factorization(ERI,
max_num_leaves=num_leaves)
tensor_error = df.factorization_error(ERI, truncated)
h = chemistry.qubit_hamiltonian(ONE_BODY,
df.reconstruct_eri(truncated),
scalar_offset=E_NUCLEAR)
lcu = PauliLCU(h)
energy_shift = ground_energy(h) - exact_energy
print(f"{num_leaves:>6} {tensor_error:>13.3e} {lcu.alpha:>8.4f} "
f"{lcu.num_terms:>6} {energy_shift:>+15.3e}")
print("\nEach row's Hamiltonian is QSVT-ready: PauliLCU(h) is a block "
"encoding of h/alpha, so a smaller alpha means a lower-degree "
"polynomial (fewer walk steps) for the same simulated time.")
if __name__ == "__main__":
main()