Add partially-randomized QPE resource estimator (prqpe)

.github/workflows/python-app.yml

@@ -89,6 +89,12 @@ jobs:
             config_file: ''
             config_setup: ''
 
+          - name: prqpe-unit-tests
+            directory: prqpe
+            script: pytest
+            config_file: ''
+            config_setup: ''
+
     name: Test ${{ matrix.test-suite.name }}
 
     steps:

analysis/algorithm.py

@@ -29,12 +29,13 @@
 
 import logging
 import math
+from dataclasses import dataclass
 
 from pyLIQTR.PhaseEstimation.pe import PhaseEstimation
 from pyLIQTR.qubitization.phase_estimation import QubitizedPhaseEstimation
 from qualtran.bloqs.phase_estimation import TextbookQPE
 
-from qhat.analysis.config_types import AlgorithmConfiguration, GeneralConfiguration
+from qhat.analysis.config_types import AlgorithmConfiguration, GeneralConfiguration, value
 
 logger = logging.getLogger(__name__)
 
@@ -140,6 +141,17 @@ def build_controlled_time_evolution(
 
 # -------------------------------------------------------------------------------------------------
 
+@dataclass(frozen=True)
+class PRQPEAlgorithm:
+    """Inert carrier of prqpe algorithm-level knobs. Builds no bloq; the prqpe estimator is analytic."""
+    method: str = "partially_randomized"
+    overlap: float = 1.0
+    xi: float | None = None
+    randomizer: str = "rte"
+    commuting_group_size: int | None = None
+
+# -------------------------------------------------------------------------------------------------
+
 def build_algorithm(
         config_algorithm: AlgorithmConfiguration,
         unitary,
@@ -158,6 +170,14 @@ def build_algorithm(
         return build_time_evolution(config_algorithm, unitary)
     elif config_algorithm.method.lower() in ("controlled time evolution",):
         return build_controlled_time_evolution(config_algorithm, unitary)
+    elif config_algorithm.method.lower() in ("qpe: partially randomized",):
+        return PRQPEAlgorithm(
+            method="partially_randomized",
+            overlap=value(config_algorithm.overlap, 1.0),
+            xi=config_algorithm.xi,
+            randomizer=value(config_algorithm.randomizer, "rte"),
+            commuting_group_size=config_algorithm.commuting_group_size,
+        )
     else:
         raise ValueError(f"Invalid algorithm method \"{config_algorithm.method}\".")
 

analysis/analysis.py

@@ -8,6 +8,7 @@
 
 from qhat.analysis.config_types import AnalysisConfiguration, GeneralConfiguration
 from qhat.analysis.file_io import save_matrix, load_state, save_state
+from qhat.analysis.prqpe_adapter import resource_estimation_prqpe
 
 logger = logging.getLogger(__name__)
 
@@ -49,12 +50,15 @@ def resource_estimation_pyliqtr(
 
 def estimate_resources(
         config_analysis: AnalysisConfiguration,
-        algorithm) -> dict:
+        algorithm,
+        hamiltonian=None) -> dict:
 
     if config_analysis.resource_estimator.lower() == "pyliqtr":
         return resource_estimation_pyliqtr(config_analysis, algorithm)
     elif config_analysis.resource_estimator.lower() == "cirq":
         return resource_estimation_cirq(config_analysis, algorithm)
+    elif config_analysis.resource_estimator.lower() == "prqpe":
+        return resource_estimation_prqpe(config_analysis, algorithm, hamiltonian=hamiltonian)
     else:
         raise ValueError(
                 f"Invalid resource estimator method \"{config_analysis.resource_estimator}\".")
@@ -590,7 +594,7 @@ def analyze_algorithm(
     # Dispatch to requested analyses
     if config_analysis.resource_estimator is not None:
         logger.info(f"Performing resource estimation using {config_analysis.resource_estimator}.")
-        results["resource_estimates"] = estimate_resources(config_analysis, algorithm)
+        results["resource_estimates"] = estimate_resources(config_analysis, algorithm, hamiltonian=hamiltonian)
 
     if config_analysis.matrix_output_file is not None:
         logger.info("Generating unitary matrix output.")

analysis/config_types.py

@@ -127,6 +127,9 @@ def encode_double_factorization(self, **kwargs):
         self._only_once()
         self.method = "double factorization"
         self.energy_error = kwargs["energy_error"]
+    def encode_none(self):
+        self._only_once()
+        self.method = "none"
     def _generate_TOML_table(self):
         table = tomlkit.table()
         table["method"] = self.method
@@ -146,12 +149,20 @@ def __init__(self):
         self.num_phase_qubits = None
         self.probability_of_failure = None
         self.energy_error = None
+        self.overlap = None
+        self.xi = None
+        self.randomizer = None
+        self.commuting_group_size = None
     def _generate_TOML_table(self):
         table = tomlkit.table()
         table["method"] = self.method
         self.save_if_present(table, "num_phase_qubits")
         self.save_if_present(table, "probability_of_failure")
         self.save_if_present(table, "energy_error")
+        self.save_if_present(table, "overlap")
+        self.save_if_present(table, "xi")
+        self.save_if_present(table, "randomizer")
+        self.save_if_present(table, "commuting_group_size")
         return table
 
 # -------------------------------------------------------------------------------------------------
@@ -166,6 +177,10 @@ def __init__(self):
         self.num_eigenvalues = 0
         self.eigendecomposition_matrices = 'approximate'
         self.which_eigenvalues = 'smallest'
+        self.prqpe_C_gs = None
+        self.prqpe_cgs_rule = None
+        self.prqpe_cgs_max_qubits = 14
+        self.prqpe_target_precision = 0.0016
 
     def _generate_TOML_table(self):
         table = tomlkit.table()
@@ -177,6 +192,11 @@ def _generate_TOML_table(self):
         self.save_if_present(table, "num_eigenvalues")
         self.save_if_present(table, "eigendecomposition_matrices")
         self.save_if_present(table, "which_eigenvalues")
+        self.save_if_present(table, "prqpe_C_gs")
+        if self.prqpe_cgs_rule is not None:
+            table["prqpe_cgs_rule"] = list(self.prqpe_cgs_rule)   # tomlkit needs a list, not a tuple
+        self.save_if_present(table, "prqpe_cgs_max_qubits")
+        self.save_if_present(table, "prqpe_target_precision")
         return table
 
 # -------------------------------------------------------------------------------------------------

analysis/examples/config_prqpe.py

@@ -0,0 +1,92 @@
+"""
+Partially-Randomized QPE Resource Estimation Example (analytic "prqpe" estimator)
+
+This configuration runs the analytic prqpe resource estimator end-to-end on a real
+molecule. Rather than constructing an explicit circuit and counting gates, prqpe
+computes Toffoli-based fault-tolerant costs (Toffoli counts, logical qubit counts,
+per-circuit Toffoli budgets, and number of circuits) directly from the Pauli
+decomposition of the Hamiltonian using closed-form formulas for partially-randomized
+QPE (Gunther et al., arXiv:2503.05647).
+
+Because the cost is analytic, the unitary encoding is a no-op (encode_none): the
+Hamiltonian is passed through unchanged to the estimator.
+
+This example uses the Be-H second-quantization tensors that ship alongside it
+(Be-H_1.30_sto-6g_as-003-003.tensors.npz, a 3-occupied/3-vacant active space).
+
+Usage (run from this examples directory so the relative .npz path resolves):
+    cd qhat/analysis/examples
+    python -m qhat.analysis.driver config_prqpe.py
+
+It writes a TOML results summary (named by a content hash) and a log file into the
+current working directory; look for the [resource_estimates] table.
+"""
+
+# =================================================================================================
+# GENERAL CONFIGURATION
+# =================================================================================================
+
+general.print_verbose()  # Options: print_default(), print_verbose(), print_debug()
+general.logfile = "prqpe_analysis.log"
+
+# =================================================================================================
+# HAMILTONIAN SPECIFICATION
+# =================================================================================================
+# Load a second-quantized Hamiltonian from file and map to qubits.
+# Supports: .npz (NumPy), .h5/.hdf5 (HDF5). The path is relative to the directory you
+# run the driver from (see the Usage note above).
+
+hamiltonian.load_second_quantization(
+    "Be-H_1.30_sto-6g_as-003-003.tensors.npz",
+    fermion_to_qubit_transform="JW"  # Jordan-Wigner
+)
+
+# To generate tensors for your own molecule, use qhat.hamiltonian_generator
+# (see qhat/hamiltonian_generator/README.md), or build them with openfermionpyscf as
+# shown in the companion script run_prqpe_molecule.py.
+#
+# Alternative: load a Pauli-string Hamiltonian directly
+# hamiltonian.load_pauli_strings("path/to/your_pauli_strings.json")
+
+# =================================================================================================
+# UNITARY ENCODING
+# =================================================================================================
+# The analytic prqpe estimator does not need an explicit circuit encoding, so the
+# unitary step is a no-op: the Hamiltonian is passed through unchanged.
+
+unitary.encode_none()
+
+# =================================================================================================
+# ALGORITHM SELECTION
+# =================================================================================================
+# Partially-randomized QPE.
+
+algorithm.method = "QPE: partially randomized"
+algorithm.overlap = 1.0  # Overlap of the initial (guiding) state with the target eigenstate
+
+# Optional partially-randomized QPE knobs:
+# algorithm.xi = None                  # depth-reduction factor (auto if None)
+# algorithm.randomizer = "rte"         # randomizer strategy ("rte" or "qdrift")
+# algorithm.commuting_group_size = None
+
+# =================================================================================================
+# ANALYSIS CONFIGURATION
+# =================================================================================================
+# Request the analytic prqpe resource estimate.
+
+analysis.resource_estimator = "prqpe"
+analysis.prqpe_target_precision = 0.0016  # Target phase-estimation precision (Hartree)
+
+# C_gs (ground-state Trotter constant):
+#   - Omit to auto-estimate on small systems (diagonalization-based, <= 14 qubits).
+#   - Set explicitly for large systems where auto-estimation is impractical:
+# analysis.prqpe_C_gs = 1.2e-3
+#   - Or extrapolate from a fitted (A, b): C_gs = A * lambda**b
+# analysis.prqpe_cgs_rule = (A, b)
+
+# Expected output (Toffoli-based costs, written to the [resource_estimates] TOML table):
+#   - toffoli_count: total Toffoli gate count
+#   - logical_qubits: number of logical qubits required
+#   - max_toffoli_per_circuit: maximum Toffoli count per circuit
+#   - num_circuits: number of circuits
+#   - C_gs / method / metadata

analysis/examples/run_prqpe_molecule.py

@@ -0,0 +1,127 @@
+#!/usr/bin/env python
+"""
+End-to-end PRQPE resource estimate for a molecule, the whole way through QHAT.
+
+This single, self-contained script:
+  1. builds a molecule's second-quantized Hamiltonian with openfermionpyscf,
+  2. writes the integral tensors in the .npz layout QHAT's loader expects,
+  3. writes a QHAT analysis config that selects the analytic "prqpe" estimator, and
+  4. runs the QHAT analysis driver on it -- the exact same entry point as
+     `python -m qhat.analysis.driver <config>` -- then prints the resource estimates.
+
+So it exercises the real, config-driven QHAT pipeline end-to-end (load second
+quantization -> Jordan-Wigner -> no-op "analytic" encoding -> partially-randomized
+QPE -> analytic prqpe resource estimate -> TOML summary).
+
+Run it with the project venv (qhat + openfermionpyscf installed):
+    python run_prqpe_molecule.py
+
+Edit MOLECULE below to estimate a different system. Small systems (<= 14 qubits)
+auto-estimate the ground-state Trotter constant C_gs by diagonalization; for larger
+systems set analysis.prqpe_C_gs in the config (see config_prqpe.py).
+"""
+
+import os
+import subprocess
+import sys
+import tempfile
+
+try:
+    import tomllib  # Python 3.11+ standard library
+except ModuleNotFoundError:  # pragma: no cover
+    import tomli as tomllib
+
+import numpy as np
+from openfermion import MolecularData
+from openfermionpyscf import run_pyscf
+
+# --- Molecule definition (edit me) -------------------------------------------------------------
+MOLECULE = {
+    "name": "H2_sto-3g",
+    "geometry": [("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.74))],  # Angstrom
+    "basis": "sto-3g",
+    "charge": 0,
+    "multiplicity": 1,  # spin multiplicity, 2S + 1
+}
+TARGET_PRECISION = 0.0016  # Hartree (chemical accuracy)
+
+
+def build_tensors_npz(path):
+    """Build the molecule and save (constant, one_body, two_body) as a QHAT-format .npz.
+
+    QHAT's load_second_quantization reads ``one_body`` and ``two_body`` (and an
+    optional scalar ``constant``) and reconstructs an openfermion InteractionOperator,
+    so we save exactly the tensors of ``get_molecular_hamiltonian()``.
+    """
+    molecule = MolecularData(
+        geometry=MOLECULE["geometry"],
+        basis=MOLECULE["basis"],
+        multiplicity=MOLECULE["multiplicity"],
+        charge=MOLECULE["charge"],
+    )
+    molecule = run_pyscf(molecule, run_scf=True)
+    ham = molecule.get_molecular_hamiltonian()  # openfermion InteractionOperator
+    np.savez_compressed(
+        path,
+        constant=ham.constant,
+        one_body=ham.one_body_tensor,
+        two_body=ham.two_body_tensor,
+    )
+    return ham.one_body_tensor.shape[0]  # number of spin-orbitals (= qubits under JW)
+
+
+CONFIG_TEMPLATE = '''\
+general.print_verbose()
+general.logfile = "prqpe_{name}.log"
+
+hamiltonian.load_second_quantization("{npz}", fermion_to_qubit_transform="JW")
+
+unitary.encode_none()
+
+algorithm.method = "QPE: partially randomized"
+algorithm.overlap = 1.0
+
+analysis.resource_estimator = "prqpe"
+analysis.prqpe_target_precision = {eps}
+'''
+
+
+def main():
+    name = MOLECULE["name"]
+    with tempfile.TemporaryDirectory() as workdir:
+        npz = os.path.join(workdir, name + ".tensors.npz")
+        n_qubits = build_tensors_npz(npz)
+        print(f"Built {name}: {n_qubits} spin-orbitals (qubits under Jordan-Wigner).")
+
+        config_path = os.path.join(workdir, "config.py")
+        with open(config_path, "w") as fh:
+            fh.write(CONFIG_TEMPLATE.format(name=name, npz=npz, eps=TARGET_PRECISION))
+
+        # Run the real QHAT driver, exactly like: python -m qhat.analysis.driver config.py
+        print("Running the QHAT analysis driver (prqpe estimator)...")
+        proc = subprocess.run(
+            [sys.executable, "-m", "qhat.analysis.driver", "config.py"],
+            cwd=workdir, capture_output=True, text=True,
+        )
+        if proc.returncode != 0:
+            sys.stderr.write(proc.stdout + "\n" + proc.stderr + "\n")
+            raise SystemExit(f"QHAT driver failed (exit {proc.returncode})")
+
+        tomls = [f for f in os.listdir(workdir) if f.endswith(".toml")]
+        if not tomls:
+            raise SystemExit("QHAT driver produced no TOML summary")
+        with open(os.path.join(workdir, tomls[0]), "rb") as fh:
+            summary = tomllib.load(fh)
+
+    est = summary["resource_estimates"]
+    print(f"\n=== PRQPE resource estimate for {name} ===")
+    print(f"  estimator              : {est['estimator']} ({est['method']})")
+    print(f"  toffoli_count          : {est['toffoli_count']:.6g}")
+    print(f"  logical_qubits         : {est['logical_qubits']}")
+    print(f"  max_toffoli_per_circuit: {est['max_toffoli_per_circuit']:.6g}")
+    print(f"  num_circuits           : {est['num_circuits']}")
+    print(f"  C_gs                   : {est['C_gs']:.6g}")
+
+
+if __name__ == "__main__":
+    main()