diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index 036fad6..6422989 100644 --- a/.github/workflows/python-app.yml +++ b/.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: diff --git a/analysis/algorithm.py b/analysis/algorithm.py index 1f2338e..94c544f 100644 --- a/analysis/algorithm.py +++ b/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}\".") diff --git a/analysis/analysis.py b/analysis/analysis.py index e5b3ba2..c26a331 100644 --- a/analysis/analysis.py +++ b/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.") diff --git a/analysis/config_types.py b/analysis/config_types.py index ab6dcaf..ab72cd0 100644 --- a/analysis/config_types.py +++ b/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 # ------------------------------------------------------------------------------------------------- diff --git a/analysis/examples/config_prqpe.py b/analysis/examples/config_prqpe.py new file mode 100644 index 0000000..0cbbc91 --- /dev/null +++ b/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 diff --git a/analysis/examples/run_prqpe_molecule.py b/analysis/examples/run_prqpe_molecule.py new file mode 100644 index 0000000..eb3b7a6 --- /dev/null +++ b/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 ` -- 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() diff --git a/analysis/prqpe_adapter.py b/analysis/prqpe_adapter.py new file mode 100644 index 0000000..c8cefa3 --- /dev/null +++ b/analysis/prqpe_adapter.py @@ -0,0 +1,131 @@ +"""QHAT -> prqpe adapter for analytic phase-estimation resource estimation. + +This module marshals a QHAT ``Hamiltonian`` into the data the analytic ``prqpe`` +estimator consumes (descending-sorted Pauli weights ``|h_l|``, the system qubit +count, and the ground-state Trotter constant ``C_gs``) and shapes the resulting +:class:`~qhat.prqpe.api.ResourceEstimate` into a flat, TOML-serializable dict. + +The estimator itself implements the partially-randomized quantum phase +estimation cost model of Gunther et al., "Phase estimation with partially +randomized time evolution" (arXiv:2503.05647); this file is only the +QHAT-side glue and does not modify the prqpe core. +""" + +from __future__ import annotations + +import logging + +import numpy as np +from openfermion import QubitOperator + +from qhat.prqpe import estimate_resources +from qhat.prqpe.representation import build_representation, sorted_weights +from qhat.prqpe.cgs import estimate_cgs, cgs_rule + +logger = logging.getLogger(__name__) + +# ------------------------------------------------------------------------------------------------- + + +def _qubit_operator_from_qhat(hamiltonian) -> QubitOperator: + """Rebuild an openfermion ``QubitOperator`` from a QHAT Hamiltonian. + + The identity term (empty Pauli tuple ``()``) carries no resource cost and is + skipped here; ``build_representation`` also drops it downstream. + """ + qop = QubitOperator() + for pauli_tuple, coef in hamiltonian.get_all_pauli_strings(return_as="tuples").items(): + if pauli_tuple == (): # identity term -> classical constant, no cost + continue + qop += QubitOperator(pauli_tuple, coef) + return qop + + +def _pick(obj, name, default): + """``getattr`` with a fallback that also treats ``None`` as "use default".""" + v = getattr(obj, name, None) + return default if v is None else v + + +def _resolve_C_gs(pauli_ham, n_qubits, lam, config_analysis) -> float: + """Resolve the ground-state Trotter constant ``C_gs`` (tiered strategy).""" + explicit = getattr(config_analysis, "prqpe_C_gs", None) + if explicit is not None: + return float(explicit) + + rule = getattr(config_analysis, "prqpe_cgs_rule", None) + if rule is not None: + A, b = rule + return float(cgs_rule(lam, A, b)) + + max_qubits = getattr(config_analysis, "prqpe_cgs_max_qubits", 14) + if n_qubits <= max_qubits: + logger.info("Estimating C_gs by diagonalization (%d qubits)." % n_qubits) + return float(estimate_cgs(pauli_ham)) + + raise ValueError( + f"C_gs required for {n_qubits} qubits (> prqpe_cgs_max_qubits={max_qubits}): " + "the system is too large to diagonalize. Set analysis.prqpe_C_gs to an " + "explicit value or analysis.prqpe_cgs_rule=(A, b) for the C_gs = A * lambda**b " + "extrapolation." + ) + + +def _to_toml_primitive(v): + """Coerce a metadata value to a TOML-serializable python primitive.""" + if isinstance(v, np.floating): + return float(v) + if isinstance(v, np.integer): + return int(v) + if isinstance(v, bool): + return v + if isinstance(v, int): + return v + if isinstance(v, float): + return v + if isinstance(v, str): + return v + if isinstance(v, (np.complexfloating, complex)): + return str(v) + if isinstance(v, np.bool_): + return bool(v) + try: + return float(v) + except (TypeError, ValueError): + return str(v) + + +# ------------------------------------------------------------------------------------------------- + + +def resource_estimation_prqpe(config_analysis, algorithm, hamiltonian=None) -> dict: + """Estimate phase-estimation resources for ``hamiltonian`` with prqpe.""" + if hamiltonian is None: + raise ValueError("PRQPE resource estimation requires the Hamiltonian.") + + pauli_ham = build_representation(_qubit_operator_from_qhat(hamiltonian)) + weights = sorted_weights(pauli_ham) + lam = float(weights.sum()) + n_qubits = hamiltonian.num_qubits() + C_gs = _resolve_C_gs(pauli_ham, n_qubits, lam, config_analysis) + + est = estimate_resources( + weights, n_qubits, C_gs, + target_precision=_pick(config_analysis, "prqpe_target_precision", 0.0016), + method=_pick(algorithm, "method", "partially_randomized"), + overlap=_pick(algorithm, "overlap", 1.0), + xi=getattr(algorithm, "xi", None), + randomizer=_pick(algorithm, "randomizer", "rte"), + commuting_group_size=getattr(algorithm, "commuting_group_size", None), + ) + + return { + "estimator": "prqpe", + "method": est.method, + "toffoli_count": float(est.toffoli_count), + "logical_qubits": int(est.logical_qubits), + "max_toffoli_per_circuit": float(est.max_toffoli_per_circuit), + "num_circuits": int(est.num_circuits), + "C_gs": float(C_gs), + "metadata": {k: _to_toml_primitive(v) for k, v in est.metadata.items()}, + } diff --git a/analysis/tests/test_prqpe_adapter.py b/analysis/tests/test_prqpe_adapter.py new file mode 100644 index 0000000..e2e1440 --- /dev/null +++ b/analysis/tests/test_prqpe_adapter.py @@ -0,0 +1,215 @@ +""" +Light unit tests for the QHAT -> prqpe adapter (qhat.analysis.prqpe_adapter). + +These tests deliberately import ONLY the adapter, the configuration types, +openfermion, and the prqpe core. They do NOT import hamiltonian.py or +analysis.py, so they run without the heavy qhat stack (pyLIQTR, qualtran, +pyscf, ...). A small duck-typed stub stands in for a real QHAT Hamiltonian. + +Coverage: +- _qubit_operator_from_qhat drops the identity term and keeps the rest. +- resource_estimation_prqpe returns a well-formed, positive, TOML-serializable + result dict. +- _resolve_C_gs tiering: explicit value > cgs_rule > small-system auto > + ValueError when the system is too large and nothing is configured. +- A missing Hamiltonian raises ValueError. +""" + +import pytest +import tomlkit +from types import SimpleNamespace + +from openfermion import QubitOperator + +from qhat.analysis.prqpe_adapter import ( + _qubit_operator_from_qhat, + _resolve_C_gs, + resource_estimation_prqpe, +) +from qhat.prqpe.representation import build_representation, sorted_weights +from qhat.prqpe.cgs import cgs_rule + + +# ================================================================================== +# Shared stub + fixtures +# ================================================================================== + +class _StubHam: + """Duck-typed stand-in for a QHAT Hamiltonian (avoids the heavy stack).""" + + def __init__(self, terms, nq): + self._t, self._nq = terms, nq + + def get_all_pauli_strings(self, return_as="tuples"): + assert return_as == "tuples" + return self._t + + def num_qubits(self): + return self._nq + + +# A small Hermitian two-qubit example with an explicit identity term. +SAMPLE_TERMS = { + (): 0.7, + ((0, 'Z'),): 1.0, + ((1, 'Z'),): -0.5, + ((0, 'X'), (1, 'X')): 0.3, +} + + +@pytest.fixture +def sample_ham(): + """A small two-qubit Hermitian stub Hamiltonian.""" + return _StubHam(dict(SAMPLE_TERMS), 2) + + +@pytest.fixture +def config_analysis(): + """A minimal analysis config selecting prqpe with small-system auto C_gs.""" + return SimpleNamespace( + resource_estimator="prqpe", + prqpe_C_gs=None, + prqpe_cgs_rule=None, + prqpe_cgs_max_qubits=14, + prqpe_target_precision=0.0016, + ) + + +@pytest.fixture +def algorithm(): + """A minimal prqpe algorithm carrier (duck-typed).""" + return SimpleNamespace( + method="partially_randomized", + overlap=1.0, + xi=None, + randomizer="rte", + commuting_group_size=None, + ) + + +# ================================================================================== +# Test: _qubit_operator_from_qhat +# ================================================================================== + +class TestQubitOperatorFromQhat: + """Test the QHAT-Hamiltonian -> openfermion QubitOperator conversion.""" + + def test_drops_identity_keeps_rest(self, sample_ham): + """The identity (empty-tuple) term is dropped; all others survive.""" + qop = _qubit_operator_from_qhat(sample_ham) + + assert isinstance(qop, QubitOperator) + # Identity term must be absent (it carries no resource cost). + assert () not in qop.terms + # The three non-identity terms must be present with their coefficients. + assert len(qop.terms) == 3 + assert qop.terms[((0, 'Z'),)] == 1.0 + assert qop.terms[((1, 'Z'),)] == -0.5 + assert qop.terms[((0, 'X'), (1, 'X'))] == 0.3 + + def test_no_identity_term_is_noop(self): + """A Hamiltonian without an identity term converts all terms.""" + ham = _StubHam({((0, 'X'),): 0.5, ((1, 'Y'),): 0.25}, 2) + qop = _qubit_operator_from_qhat(ham) + + assert () not in qop.terms + assert len(qop.terms) == 2 + + +# ================================================================================== +# Test: resource_estimation_prqpe +# ================================================================================== + +class TestResourceEstimationPrqpe: + """Test the top-level adapter entry point.""" + + def test_returns_well_formed_dict(self, config_analysis, algorithm, sample_ham): + """A successful estimate is a flat dict with positive, typed fields.""" + out = resource_estimation_prqpe(config_analysis, algorithm, hamiltonian=sample_ham) + + assert out["estimator"] == "prqpe" + + assert isinstance(out["toffoli_count"], float) + assert out["toffoli_count"] > 0 + + assert isinstance(out["logical_qubits"], int) + assert out["logical_qubits"] > 0 + + assert isinstance(out["C_gs"], float) + assert out["C_gs"] > 0 + + assert isinstance(out["metadata"], dict) + + def test_result_is_toml_serializable(self, config_analysis, algorithm, sample_ham): + """The result dict round-trips through tomlkit (wrapped in a table).""" + out = resource_estimation_prqpe(config_analysis, algorithm, hamiltonian=sample_ham) + + document = tomlkit.document() + document["resource_estimates"] = out + rendered = tomlkit.dumps(document) + reloaded = tomlkit.loads(rendered) + + assert reloaded["resource_estimates"]["estimator"] == "prqpe" + assert reloaded["resource_estimates"]["logical_qubits"] == out["logical_qubits"] + assert reloaded["resource_estimates"]["toffoli_count"] == pytest.approx( + out["toffoli_count"] + ) + + def test_missing_hamiltonian_raises(self, config_analysis, algorithm): + """A missing Hamiltonian is a usage error.""" + with pytest.raises(ValueError, match="requires the Hamiltonian"): + resource_estimation_prqpe(config_analysis, algorithm, hamiltonian=None) + + +# ================================================================================== +# Test: _resolve_C_gs tiering +# ================================================================================== + +class TestResolveCgs: + """Test the tiered ground-state Trotter constant resolution.""" + + def test_explicit_value_wins(self): + """An explicit prqpe_C_gs takes precedence over everything else.""" + config = SimpleNamespace( + prqpe_C_gs=3.14, + prqpe_cgs_rule=(0.5, 2.0), # present, but must be ignored + prqpe_cgs_max_qubits=14, + ) + # pauli_ham is unused on this path, so None is fine. + assert _resolve_C_gs(None, 2, 1.0, config) == 3.14 + + def test_cgs_rule_path(self): + """With no explicit value, the (A, b) rule yields cgs_rule(lam, A, b).""" + A, b, lam = 0.5, 2.0, 3.0 + config = SimpleNamespace( + prqpe_C_gs=None, + prqpe_cgs_rule=(A, b), + prqpe_cgs_max_qubits=14, + ) + result = _resolve_C_gs(None, 2, lam, config) + assert result == pytest.approx(cgs_rule(lam, A, b)) + assert result == pytest.approx(A * lam ** b) + + def test_small_system_auto_path(self, sample_ham): + """A small system with nothing set diagonalizes for a positive C_gs.""" + pauli_ham = build_representation(_qubit_operator_from_qhat(sample_ham)) + lam = float(sorted_weights(pauli_ham).sum()) + config = SimpleNamespace( + prqpe_C_gs=None, + prqpe_cgs_rule=None, + prqpe_cgs_max_qubits=14, + ) + result = _resolve_C_gs(pauli_ham, sample_ham.num_qubits(), lam, config) + assert isinstance(result, float) + assert result > 0 + + def test_too_large_without_config_raises(self, sample_ham): + """Beyond prqpe_cgs_max_qubits, with nothing else set, raises ValueError.""" + pauli_ham = build_representation(_qubit_operator_from_qhat(sample_ham)) + config = SimpleNamespace( + prqpe_C_gs=None, + prqpe_cgs_rule=None, + prqpe_cgs_max_qubits=1, # below the 2-qubit system size + ) + with pytest.raises(ValueError, match="C_gs required"): + _resolve_C_gs(pauli_ham, 2, 1.0, config) diff --git a/analysis/tests/test_prqpe_integration.py b/analysis/tests/test_prqpe_integration.py new file mode 100644 index 0000000..f0cb6bc --- /dev/null +++ b/analysis/tests/test_prqpe_integration.py @@ -0,0 +1,107 @@ +""" +End-to-end integration tests for the prqpe estimator through the real QHAT seams. + +Unlike test_prqpe_adapter.py, this module exercises the genuine pipeline objects +and therefore imports the full qhat stack (hamiltonian.py, analysis.py, +algorithm.py, unitary.py). It is intended to be run separately, with the heavy +dependencies (pyLIQTR, qualtran, pyscf, ...) available. + +A real Hamiltonian is built from data_pauli.json using the same loader pattern +as test_pauli_hamiltonian.py. The test then walks the three integration seams: + +- the unitary "none" no-op returns the Hamiltonian unchanged; +- the algorithm "QPE: partially randomized" method builds a PRQPEAlgorithm carrier; +- analyze_algorithm with resource_estimator="prqpe" produces a prqpe estimate. +""" + +import os +import pytest + +from qhat.analysis.config_types import ( + AlgorithmConfiguration, + AnalysisConfiguration, + GeneralConfiguration, + GeneralConfigurationUser, + HamiltonianConfiguration, + UnitaryConfiguration, +) +from qhat.analysis.hamiltonian import load_pauli +from qhat.analysis.unitary import encode_as_unitary +from qhat.analysis.algorithm import PRQPEAlgorithm, build_algorithm +from qhat.analysis.analysis import analyze_algorithm + + +# ================================================================================== +# Shared fixtures +# ================================================================================== + +@pytest.fixture(scope="session") +def general_config(): + """Create a single GeneralConfiguration for the entire test session.""" + user_config = GeneralConfigurationUser() + return GeneralConfiguration(user_config) + + +@pytest.fixture +def real_hamiltonian(general_config): + """Load the real 4-qubit Hamiltonian from data_pauli.json.""" + test_file = os.path.join( + os.path.dirname(os.path.abspath(__file__)), + "data_pauli.json", + ) + config_ham = HamiltonianConfiguration() + config_ham.load_pauli_strings(test_file) + return load_pauli(config_ham) + + +# ================================================================================== +# Test: full prqpe integration through the real seams +# ================================================================================== + +class TestPrqpeIntegration: + """Exercise the unitary, algorithm, and analysis seams end to end.""" + + def test_unitary_none_returns_hamiltonian_unchanged(self, real_hamiltonian): + """The 'none' unitary method is a no-op: it returns the Hamiltonian.""" + config_unitary = UnitaryConfiguration() + config_unitary.encode_none() + assert config_unitary.method == "none" + + result = encode_as_unitary(config_unitary, real_hamiltonian, None) + assert result is real_hamiltonian + + def test_build_algorithm_returns_prqpe_carrier(self, real_hamiltonian): + """'QPE: partially randomized' builds a PRQPEAlgorithm carrier.""" + config_algorithm = AlgorithmConfiguration() + config_algorithm.method = "QPE: partially randomized" + config_algorithm.overlap = 1.0 + + carrier = build_algorithm(config_algorithm, real_hamiltonian, None) + + assert isinstance(carrier, PRQPEAlgorithm) + assert carrier.method == "partially_randomized" + assert carrier.overlap == 1.0 + assert carrier.randomizer == "rte" + + def test_analyze_algorithm_produces_prqpe_estimate(self, real_hamiltonian): + """analyze_algorithm dispatches to prqpe and yields a positive estimate.""" + config_algorithm = AlgorithmConfiguration() + config_algorithm.method = "QPE: partially randomized" + config_algorithm.overlap = 1.0 + carrier = build_algorithm(config_algorithm, real_hamiltonian, None) + + config_analysis = AnalysisConfiguration() + config_analysis.resource_estimator = "prqpe" + config_analysis.prqpe_target_precision = 0.0016 + # Auto C_gs: the 4-qubit system is well below prqpe_cgs_max_qubits. + + results = analyze_algorithm( + config_analysis, carrier, hamiltonian=real_hamiltonian + ) + + assert "resource_estimates" in results + estimate = results["resource_estimates"] + assert estimate["estimator"] == "prqpe" + assert estimate["toffoli_count"] > 0 + assert estimate["logical_qubits"] > 0 + assert estimate["C_gs"] > 0 diff --git a/analysis/unitary.py b/analysis/unitary.py index a23204c..d63e671 100644 --- a/analysis/unitary.py +++ b/analysis/unitary.py @@ -198,5 +198,7 @@ def encode_as_unitary( raise NotImplementedError() elif config_unitary.method.lower() in ("ramped trotter"): return encode_ramped_trotter(config_unitary, hamiltonian, tevol_hbar) + elif config_unitary.method.lower() in ("none", "analytic"): + return hamiltonian else: raise ValueError(f"Invalid unitary encoding method \"{config_unitary.method}\".") diff --git a/prqpe/__init__.py b/prqpe/__init__.py new file mode 100644 index 0000000..6b6ff28 --- /dev/null +++ b/prqpe/__init__.py @@ -0,0 +1,17 @@ +"""Resource estimation for ground-state phase estimation with partially +randomized product formulas. The public entry point is :func:`estimate_resources`. +""" + +from __future__ import annotations + +from .api import ( + ResourceEstimate, + SplitEstimate, + estimate_resources, +) + +__all__ = [ + "estimate_resources", + "ResourceEstimate", + "SplitEstimate", +] diff --git a/prqpe/api.py b/prqpe/api.py new file mode 100644 index 0000000..a61cfa1 --- /dev/null +++ b/prqpe/api.py @@ -0,0 +1,244 @@ +"""Public entry point: estimate phase-estimation resources from term weights. + +``estimate_resources`` takes the cost-relevant *data* about a Hamiltonian -- its +sorted Pauli weights ``|h_l|``, the system qubit count, and the ground-state +Trotter constant ``C_gs`` -- and returns a :class:`ResourceEstimate`. Weight +reduction, the fermion-to-qubit mapping, and the ``C_gs`` measurement all happen +upstream. + +The estimator sweeps the deterministic/randomized split ``L_D`` (the number of +largest-weight terms compiled deterministically), jointly optimizes the Trotter +step and randomization scale at each split, and keeps the whole sweep as a +Pareto curve alongside the cheapest point. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Literal + +import numpy as np + +from . import costs, schedule + +Method = Literal["deterministic", "randomized", "partially_randomized"] +Randomizer = Literal["rte", "qdrift"] + +#: Default commuting-group length for Hamming-weight phasing. +DEFAULT_COMMUTING_GROUP_SIZE = 10 + +#: qDRIFT draws half as many rotations as RTE for the same accuracy. +_RANDOMIZER_FACTOR = {"rte": 1.0, "qdrift": 0.5} + + +# ---------------------------------------------------------------------------- +# Result containers. +# ---------------------------------------------------------------------------- + +@dataclass(frozen=True) +class SplitEstimate: + """Cost of one deterministic/randomized split -- one point on the sweep.""" + + L_D: int + lambda_R: float + delta: float + kappa: float + eps_qpe: float + toffoli_total: float + max_toffoli_per_circuit: float + num_circuits: int + logical_qubits: int + + +@dataclass(frozen=True) +class ResourceEstimate: + """The cheapest split, plus the full Pareto curve over ``L_D``.""" + + toffoli_count: float + logical_qubits: int + max_toffoli_per_circuit: float + num_circuits: int # one guiding-state preparation per circuit + method: str + metadata: dict = field(default_factory=dict) + pareto: list[SplitEstimate] = field(default_factory=list) + + +# ---------------------------------------------------------------------------- +# Split selection. +# ---------------------------------------------------------------------------- + +def _candidate_splits(L: int, n_points: int = 200) -> np.ndarray: + """Values of ``L_D`` to evaluate: log-spaced over ``[1, L]`` plus ``0`` and ``L``. + + The terms are sorted by descending weight, so ``L_D`` interpolates between + the fully randomized (``L_D = 0``) and fully deterministic (``L_D = L``) + endpoints. A geometric spacing resolves the interior minimum without + evaluating every integer split. + """ + if L <= 1: + return np.array([0, L][: L + 1], dtype=np.int64) + interior = np.round(np.geomspace(1, L, n_points)).astype(np.int64) + return np.unique(np.concatenate([[0], interior, [L]])) + + +# ---------------------------------------------------------------------------- +# Per-split cost. +# ---------------------------------------------------------------------------- + +def _estimate_for_split(L_D: int, lam_R: float, n_system_qubits: int, lam: float, + C_gs: float, *, eps: float, order: int, K: int, + overlap: float, xi: float, randomizer: Randomizer, + eps_synth: float, n_stage: int = 2) -> tuple[SplitEstimate, dict]: + """Optimize the cost at one split and return ``(SplitEstimate, metadata)``. + + The randomized part has weight ``lam_R`` and dyadic-angle exponent ``J``; the + deterministic part has ``L_D`` terms whose synthesis cost ``G_det`` depends on + the chosen Trotter step ``delta``, which in turn depends on ``G_det`` through + the joint optimum -- a small fixed point resolved in a couple of iterations. + """ + randomizer_factor = _RANDOMIZER_FACTOR[randomizer] + + # Randomized part: dyadic-angle exponent and per-rotation Toffoli cost. + J = costs.phase_state_precision(lam_R if lam_R > 0.0 else lam, eps) + G_rand = costs.rotation_toffolis(J, K) + + # Deterministic part: per-rotation synthesis cost. + if L_D > 0: + delta_seed = costs.trotter_step_size(eps, C_gs, order) + eps_prime = costs.synthesis_precision_per_rotation( + delta_seed, eps_synth, n_stage, L_D) + G_det = costs.synthesis_toffolis(eps_prime) + for _ in range(2): + opt = costs.minimize_delta_kappa( + L_D, lam_R, eps, C_gs, G_det, G_rand, order=order, + n_stage=n_stage, randomizer_factor=randomizer_factor) + if opt.delta > 0.0: + eps_prime = costs.synthesis_precision_per_rotation( + opt.delta, eps_synth, n_stage, L_D) + G_det = costs.synthesis_toffolis(eps_prime) + opt = costs.minimize_delta_kappa( + L_D, lam_R, eps, C_gs, G_det, G_rand, order=order, + n_stage=n_stage, randomizer_factor=randomizer_factor) + else: + G_det = 0.0 + opt = costs.minimize_delta_kappa( + 0, lam_R, eps, C_gs, G_det, G_rand, order=order, + n_stage=n_stage, randomizer_factor=randomizer_factor) + + # Round count and circuit schedule (at least one round for any input). + if L_D == 0: + M = schedule.num_rounds_randomized(lam_R, opt.eps_qpe, xi) + else: + M = schedule.num_rounds(opt.eps_qpe, opt.delta, xi) + M = max(M, 0) + sched = schedule.rpe_schedule(M, opt.kappa, overlap) + + # Deepest single circuit. With no Trotter step the randomized arm reads it + # off the closed-form total; otherwise it is the depth-M circuit explicitly. + if L_D == 0: + max_per_circuit = costs.randomized_deepest_circuit_toffolis( + opt.toffoli_total) + else: + max_per_circuit = schedule.max_toffolis_per_circuit( + sched, G_det=G_det, G_rand=G_rand * randomizer_factor, L_D=L_D, + lam_R=lam_R, delta=opt.delta, kappa=opt.kappa, n_stage=n_stage) + + # Qubits: a randomized part needs the phasing register; a purely + # deterministic arm only needs the system and the Hadamard-test ancilla. + if lam_R > 0.0: + qubits = costs.logical_qubits(n_system_qubits, J, K) + else: + qubits = n_system_qubits + 1 + + estimate = SplitEstimate( + L_D=int(L_D), lambda_R=float(lam_R), delta=opt.delta, kappa=opt.kappa, + eps_qpe=opt.eps_qpe, toffoli_total=opt.toffoli_total, + max_toffoli_per_circuit=float(max_per_circuit), + num_circuits=sched.num_circuits, logical_qubits=int(qubits)) + meta = dict(M=M, J=J, K=K, G_det=G_det, G_rand=G_rand) + return estimate, meta + + +# ---------------------------------------------------------------------------- +# Entry point. +# ---------------------------------------------------------------------------- + +def estimate_resources(weights, n_system_qubits: int, C_gs: float, *, + target_precision: float = 0.0016, + method: Method = "partially_randomized", + representation: str = "pauli", + order: int = 2, overlap: float = 1.0, + xi: float | None = None, + randomizer: Randomizer = "rte", + synthesis_precision: float = 1e-4, + commuting_group_size: int | None = None) -> ResourceEstimate: + """Estimate Toffoli and qubit costs for ground-state phase estimation. + + Parameters + ---------- + weights : array of Pauli coefficients ``|h_l|`` (sorted internally, descending). + n_system_qubits : qubits carrying the Hamiltonian. + C_gs : ground-state Trotter constant, ``|E0 - E0~| = C_gs * delta^p``. + target_precision : energy precision ``eps`` (Hartree). + method : ``deterministic`` (L_D = L), ``randomized`` (L_D = 0), or + ``partially_randomized`` (the full ``L_D`` sweep). + representation : only ``pauli`` is implemented. + overlap, xi : guiding-state overlap and depth-reduction factor. + randomizer : ``rte`` (unbiased) or ``qdrift`` (half as many rotations). + synthesis_precision : ground-state error budget for rotation synthesis. + commuting_group_size : ``K`` for Hamming-weight phasing (default 10). + + Returns + ------- + ResourceEstimate + The cheapest split, with the full Pareto curve over ``L_D`` in ``pareto``. + """ + if representation != "pauli": + raise ValueError( + f"unknown representation {representation!r}; only 'pauli' is implemented") + if randomizer not in _RANDOMIZER_FACTOR: + raise ValueError(f"unknown randomizer {randomizer!r}") + + w = np.sort(np.abs(np.asarray(weights, dtype=float)))[::-1] + L = int(w.size) + if L == 0: + raise ValueError("weights is empty") + lam = float(w.sum()) + K = commuting_group_size or DEFAULT_COMMUTING_GROUP_SIZE + xi_value = schedule.resolve_depth_tradeoff(overlap, xi) + + if method == "deterministic": + ld_values = np.array([L], dtype=np.int64) + elif method == "randomized": + ld_values = np.array([0], dtype=np.int64) + elif method == "partially_randomized": + ld_values = _candidate_splits(L) + else: + raise ValueError(f"unknown method {method!r}") + + # Tail weight: lambda_R(L_D) = sum of the weights not treated deterministically. + tail = np.concatenate([np.cumsum(w[::-1])[::-1], [0.0]]) + + results = [ + _estimate_for_split( + int(L_D), float(tail[L_D]), n_system_qubits, lam, C_gs, + eps=target_precision, order=order, K=K, overlap=overlap, + xi=xi_value, randomizer=randomizer, eps_synth=synthesis_precision) + for L_D in ld_values + ] + + pareto = [est for est, _ in results] + best_est, best_meta = min(results, key=lambda r: r[0].toffoli_total) + + metadata = dict( + L_D=best_est.L_D, lambda_R=best_est.lambda_R, delta=best_est.delta, + M=best_meta["M"], kappa=best_est.kappa, xi=xi_value, + eps_qpe=best_est.eps_qpe, J=best_meta["J"], K=best_meta["K"], + G_det=best_meta["G_det"], G_rand=best_meta["G_rand"]) + + return ResourceEstimate( + toffoli_count=best_est.toffoli_total, + logical_qubits=best_est.logical_qubits, + max_toffoli_per_circuit=best_est.max_toffoli_per_circuit, + num_circuits=best_est.num_circuits, + method=method, metadata=metadata, pareto=pareto) diff --git a/prqpe/bliss/__init__.py b/prqpe/bliss/__init__.py new file mode 100644 index 0000000..0bd8232 --- /dev/null +++ b/prqpe/bliss/__init__.py @@ -0,0 +1,443 @@ +"""Hamiltonian-weight reduction for ground-state phase estimation. + +The gate cost scales with the Hamiltonian 1-norm ``lambda = sum_l |h_l|``, which +for an electronic Hamiltonian in spatial-orbital chemists'-notation integrals +``(h1, g)`` has a closed form (see :func:`lambda_pauli`). Two freedoms shrink it +without changing the physics on the target ``n``-electron sector: + +* an **orbital rotation** ``U`` (a basis change, so it preserves the spectrum), and +* a **block-invariant shift** by the number-operator deficit ``(n_hat - n)``, + which is identically zero on the ``n``-electron sector but cancels 1-norm on + the other particle-number sectors. + +The entry point :func:`reduce_hamiltonian_weight` jointly minimizes the +closed-form 1-norm over both freedoms with an annealed smoothing of ``|.|`` and +L-BFGS-B. The shift follows Loaiza and Izmaylov, Quantum Sci. Technol. 8, +035019 (2023). Depends only on ``numpy`` and ``scipy``. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field + +import numpy as np +import scipy.linalg as sla +from scipy.optimize import minimize + +__all__ = ["ReducedIntegrals", "reduce_hamiltonian_weight"] + +#: Default smoothing-parameter schedule for the annealed objective. Each stage +#: replaces ``|x|`` by ``sqrt(x^2 + mu^2) - mu`` and warm-starts from the +#: previous stage; the final ``mu = 0`` stage is the exact (nonsmooth) 1-norm. +DEFAULT_MU_SCHEDULE = (1e-1, 1e-2, 1e-3, 1e-5, 0.0) + + +# ---------------------------------------------------------------------------- +# Closed-form Pauli 1-norm. +# ---------------------------------------------------------------------------- + +def effective_one_body(h1: np.ndarray, g: np.ndarray) -> np.ndarray: + """The effective one-body matrix ``T`` folding the two-body mean field. + + With spatial-orbital chemists'-notation integrals ``g = (pq|rs)`` the + one-body part of the normal-ordered Hamiltonian is + + T = h1 + sum_r g[:, :, r, r] - 0.5 * sum_r g[:, r, r, :]. + """ + return (h1 + + np.einsum("pqrr->pq", g) + - 0.5 * np.einsum("prrq->pq", g)) + + +def lambda_pauli(h1: np.ndarray, g: np.ndarray) -> dict: + """Closed-form Jordan-Wigner / Bravyi-Kitaev 1-norm of ``(h1, g)``. + + Inputs are spatial-orbital, real, chemists'-notation integrals: ``h1`` is the + symmetric bare one-body matrix ``(N, N)`` and ``g = (pq|rs)`` is the + 8-fold-symmetric two-electron tensor ``(N, N, N, N)``. The 1-norm splits as + + lambda = lambda_T + lambda_V2 + lambda_V3, + + a one-body term ``lambda_T = sum |T|`` with ``T`` from + :func:`effective_one_body`, an opposite-spin two-body term + ``lambda_V2 = sum |g| / 4``, and a same-spin two-body term + ``lambda_V3 = sum_{p!=r, q!=s} |D| / 8`` with the antisymmetrized + ``D[p, q, r, s] = g[p, q, r, s] - g[p, s, r, q]``. Returns a dict with the + three parts and their ``total``. + """ + N = h1.shape[0] + T = effective_one_body(h1, g) + D = g - g.transpose(0, 3, 2, 1) + off = ~np.eye(N, dtype=bool) + mask = off[:, None, :, None] & off[None, :, None, :] # p != r AND q != s + lam_T = float(np.abs(T).sum()) + lam_V2 = float(np.abs(g).sum() / 4.0) + lam_V3 = float(np.abs(D)[mask].sum() / 8.0) + return { + "one_body": lam_T, + "opposite_spin": lam_V2, + "same_spin": lam_V3, + "total": lam_T + lam_V2 + lam_V3, + } + + +# ---------------------------------------------------------------------------- +# Integral transforms: orbital rotation and block-invariant shift. +# ---------------------------------------------------------------------------- + +def rotate_integrals(h1: np.ndarray, g: np.ndarray, + U: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Rotate ``(h1, g)`` into the orbital basis defined by orthogonal ``U``. + + ``U`` mixes spatial orbitals; the change of basis is + + h1_rot = U.T @ h1 @ U, + g_rot[a, b, c, d] = sum g[p, q, r, s] U[p, a] U[q, b] U[r, c] U[s, d]. + + The four-index contraction is staged as two einsums for cost. An orbital + rotation is a unitary change of basis and preserves the full spectrum, so the + core energy is unchanged (handled by the caller). + """ + h1_rot = U.T @ h1 @ U + g_rot = np.einsum("pqrs,pa,qb->abrs", g, U, U, optimize=True) + g_rot = np.einsum("abrs,rc,sd->abcd", g_rot, U, U, optimize=True) + return h1_rot, g_rot + + +def shifted_integrals(h1: np.ndarray, g: np.ndarray, s: np.ndarray, + nu: float, c: float, + n_elec: float) -> tuple[np.ndarray, np.ndarray]: + """Block-invariant shift of ``(h1, g)`` by symmetry-operator parameters. + + The shift adds ``Delta H`` built from the number-operator deficit + ``(n_hat - n)`` with a symmetric one-body operator ``s`` and scalars + ``nu, c``. Expanded into integral coefficients (``I`` the ``N x N`` identity, + ``n = n_elec``): + + gp = g + einsum('pq,rs->pqrs', s, I) + einsum('pq,rs->pqrs', I, s) + + nu * einsum('pq,rs->pqrs', I, I), + hp = h1 + (1 - n) * s + (nu * (0.5 - n) + c) * I. + + Returns ``(hp, gp)`` only; the accompanying core-energy constant + ``nu * n^2 / 2 - c * n`` is added by the caller. + """ + N = h1.shape[0] + eye = np.eye(N) + gp = (g + + np.einsum("pq,rs->pqrs", s, eye) + + np.einsum("pq,rs->pqrs", eye, s) + + nu * np.einsum("pq,rs->pqrs", eye, eye)) + hp = h1 + (1.0 - n_elec) * s + (nu * (0.5 - n_elec) + c) * eye + return hp, gp + + +def shift_core_constant(nu: float, c: float, n_elec: float) -> float: + """The core-energy constant added by the shift: ``nu * n^2 / 2 - c * n``.""" + return nu * n_elec ** 2 / 2.0 - c * n_elec + + +# ---------------------------------------------------------------------------- +# Smoothed objective and analytic gradient. +# ---------------------------------------------------------------------------- + +def _sabs(x: np.ndarray, mu: float) -> np.ndarray: + """Smoothed absolute value ``sqrt(x^2 + mu^2) - mu`` (equals ``|x|`` at mu=0).""" + if mu == 0.0: + return np.abs(x) + return np.sqrt(x * x + mu * mu) - mu + + +def _sgrad(x: np.ndarray, mu: float) -> np.ndarray: + """Derivative of :func:`_sabs`: ``x / sqrt(x^2 + mu^2)`` (sign at mu=0).""" + if mu == 0.0: + return np.sign(x) + return x / np.sqrt(x * x + mu * mu) + + +def _lambda_smooth_value(hp: np.ndarray, gp: np.ndarray, mu: float, + mask: np.ndarray) -> float: + """Smoothed 1-norm on shifted integrals ``(hp, gp)`` at smoothing ``mu``.""" + T = effective_one_body(hp, gp) + D = gp - gp.transpose(0, 3, 2, 1) + lam_T = _sabs(T, mu).sum() + lam_V2 = 0.25 * _sabs(gp, mu).sum() + lam_V3 = 0.125 * (_sabs(D, mu) * mask).sum() + return float(lam_T + lam_V2 + lam_V3) + + +def _lambda_smooth_grad_integrals(hp: np.ndarray, gp: np.ndarray, mu: float, + mask: np.ndarray + ) -> tuple[np.ndarray, np.ndarray]: + """Cotangents ``(dL/dhp, dL/dgp)`` of the smoothed 1-norm on ``(hp, gp)``.""" + N = hp.shape[0] + eye = np.eye(N) + T = effective_one_body(hp, gp) + D = gp - gp.transpose(0, 3, 2, 1) + + dT = _sgrad(T, mu) # dL/dT + dgp = 0.25 * _sgrad(gp, mu) # opposite-spin term contribution + + # one-body term backprop: T = h1 + einsum('pqrr->pq') - 0.5 einsum('prrq->pq') + dhp = dT.copy() + dgp += np.einsum("pq,rs->pqrs", dT, eye) # d g[p, q, r, r] += dT[p, q] + dgp += -0.5 * np.einsum("ps,qr->pqrs", dT, eye) # d g[p, r, r, q] += -0.5 dT[p, q] + + # same-spin term backprop: D = gp - gp.transpose(0, 3, 2, 1) + dD = 0.125 * (_sgrad(D, mu) * mask) + dgp += dD - dD.transpose(0, 3, 2, 1) + return dhp, dgp + + +# ---------------------------------------------------------------------------- +# Parameter packing. +# ---------------------------------------------------------------------------- + +@dataclass(frozen=True) +class _Layout: + """Flat-vector layout for the optimization parameters ``theta``. + + ``theta = concat([X_lower (if orbital), s_upper, [nu], [c]])`` with + ``X_lower`` the strict lower triangle of the antisymmetric rotation generator + and ``s_upper`` the upper triangle (incl. diagonal) of the symmetric shift + matrix. + """ + + N: int + optimize_orbitals: bool + + @property + def itril(self): + return np.tril_indices(self.N, -1) + + @property + def iu(self): + return np.triu_indices(self.N) + + @property + def nX(self): + return self.N * (self.N - 1) // 2 if self.optimize_orbitals else 0 + + @property + def nS(self): + return self.N * (self.N + 1) // 2 + + @property + def n_par(self): + return self.nX + self.nS + 2 + + def unpack(self, theta): + """``theta -> (X, s, nu, c)`` with antisymmetric ``X`` and symmetric ``s``.""" + N = self.N + theta = np.asarray(theta) + X = np.zeros((N, N), dtype=theta.dtype) + off = 0 + if self.optimize_orbitals: + X[self.itril] = theta[:self.nX] + X = X - X.T + off = self.nX + s = np.zeros((N, N), dtype=theta.dtype) + s[self.iu] = theta[off:off + self.nS] + s = s + s.T - np.diag(np.diag(s)) + nu = theta[off + self.nS] + c = theta[off + self.nS + 1] + return X, s, nu, c + + +# ---------------------------------------------------------------------------- +# The full pipeline (rotate then shift) and its gradient. +# ---------------------------------------------------------------------------- + +def _pipeline(theta, layout, h1, g, n_elec): + """Apply rotation then shift; return ``(hp, gp, U, X, s, nu, c)``.""" + N = layout.N + X, s, nu, c = layout.unpack(theta) + if layout.optimize_orbitals: + U = sla.expm(X) # X is already the antisymmetric generator (X - X^T) + h_r, g_r = rotate_integrals(h1, g, U) + else: + U = np.eye(N) + h_r, g_r = h1, g + hp, gp = shifted_integrals(h_r, g_r, s, nu, c, n_elec) + return hp, gp, U, X, s, nu, c + + +def _shift_block_grad(dhp, dgp, layout, n_elec): + """Backprop integral cotangents through the shift to ``(d s_upper, d nu, d c)``.""" + N = layout.N + eye = np.eye(N) + # hp = h_r + (1 - n) s + (nu (0.5 - n) + c) I + # gp = g_r + s (x) I + I (x) s + nu I (x) I + ds = ((1.0 - n_elec) * dhp + + np.einsum("pqrs,rs->pq", dgp, eye) + + np.einsum("pqrs,pq->rs", dgp, eye)) + # fold the lower triangle onto the upper (diagonal once) to match the symmetric parametrization + ds_sym = ds + ds.T - np.diag(np.diag(ds)) + d_s_upper = ds_sym[layout.iu] + d_nu = float((0.5 - n_elec) * np.trace(dhp) + + np.einsum("pqrs,pq,rs->", dgp, eye, eye)) + d_c = float(np.trace(dhp)) + return d_s_upper, d_nu, d_c + + +def _make_funcs(layout, h1, g, n_elec, mask): + """Build ``(value, grad)`` closures for the smoothed objective at a given mu. + + The shift-block gradient is analytic. The orbital-block gradient (the + expm pushforward) is obtained by a complex-step derivative on the rotation + parameters only, exact to machine precision. + """ + N = layout.N + + def value(theta, mu): + hp, gp, *_ = _pipeline(theta, layout, h1, g, n_elec) + return _lambda_smooth_value(hp, gp, mu, mask) + + def grad(theta, mu): + theta = np.asarray(theta, dtype=float) + hp, gp, U, X, s, nu, c = _pipeline(theta, layout, h1, g, n_elec) + dhp, dgp = _lambda_smooth_grad_integrals(hp, gp, mu, mask) + d_s_upper, d_nu, d_c = _shift_block_grad(dhp, dgp, layout, n_elec) + + if layout.optimize_orbitals: + d_X = np.zeros(layout.nX) + step = 1e-30 + base = theta.copy() + for k in range(layout.nX): + pert = base.astype(complex) + pert[k] += 1j * step + val = _value_complex(pert, layout, h1, g, n_elec, mu, mask) + d_X[k] = val.imag / step + return np.concatenate([d_X, d_s_upper, [d_nu], [d_c]]) + return np.concatenate([d_s_upper, [d_nu], [d_c]]) + + return value, grad + + +def _value_complex(theta, layout, h1, g, n_elec, mu, mask): + """Complex-step-friendly smoothed value (holomorphic sqrt, not np.abs). + + Used only for the rotation-block complex-step gradient, so mu is the + positive smoothing of the current anneal stage and sqrt is holomorphic. + """ + hp, gp, *_ = _pipeline(theta, layout, h1, g, n_elec) + T = effective_one_body(hp, gp) + D = gp - gp.transpose(0, 3, 2, 1) + sab = lambda z: np.sqrt(z * z + mu * mu) - mu + return sab(T).sum() + 0.25 * sab(gp).sum() + 0.125 * (sab(D) * mask).sum() + + +# ---------------------------------------------------------------------------- +# Result container. +# ---------------------------------------------------------------------------- + +@dataclass(frozen=True) +class ReducedIntegrals: + """Reduced integrals and provenance from :func:`reduce_hamiltonian_weight`.""" + + h1: np.ndarray # (N, N) reduced bare one-body integrals (rotated + shifted) + g: np.ndarray # (N, N, N, N) reduced chemists' (pq|rs), 8-fold symmetric + e_core: float # core energy incl. shift constant + n_elec: int # target electron number (the preserved sector) + lambda_reduced: float # final closed-form Pauli 1-norm at the optimum + lambda_base: float # 1-norm of the input integrals (no rotation/shift) + U: np.ndarray = field(default=None) # orbital rotation actually applied + s: np.ndarray = field(default=None) # symmetric shift matrix at optimum + nu: float = 0.0 # quadratic number shift at optimum + c: float = 0.0 # linear number shift at optimum + + +# ---------------------------------------------------------------------------- +# Public entry point. +# ---------------------------------------------------------------------------- + +def reduce_hamiltonian_weight( + integrals, + *, + optimize_orbitals: bool = True, + nu_nonneg: bool = False, + mu_schedule=DEFAULT_MU_SCHEDULE, + maxiter: int = 3000, + x0: np.ndarray = None, + rng_seed: int = 0, +) -> ReducedIntegrals: + """Minimize the closed-form Pauli 1-norm over orbital and shift freedoms. + + ``integrals`` exposes ``.h1`` ``(N, N)``, ``.g`` ``(N, N, N, N)``, ``.e_core`` + and ``.n_elec``. The 1-norm of :func:`lambda_pauli` is minimized over an + orbital rotation ``U = expm(X - X^T)`` and a block-invariant shift + ``(s, nu, c)`` on the rotated-then-shifted integrals; the nonsmooth ``|.|`` is + annealed through ``mu_schedule`` with L-BFGS-B (warm-started between stages), + the final stage being the exact 1-norm. The returned :class:`ReducedIntegrals` + carries the reduced spatial-orbital integrals and the shifted core energy + ``e_core + nu * n^2 / 2 - c * n``; mapping them to qubits is the caller's job. + + Parameters + ---------- + integrals : object with ``.h1``, ``.g``, ``.e_core``, ``.n_elec``. + optimize_orbitals : include the orbital rotation (else fit the shift only). + nu_nonneg : constrain the quadratic number-shift scalar ``nu >= 0``. + mu_schedule : decreasing smoothing parameters; the final value should be ``0``. + maxiter : maximum L-BFGS-B iterations per annealing stage. + x0 : explicit flat start vector of length ``n_par`` (overrides ``rng_seed``). + rng_seed : seed for the small random start perturbation; distinct seeds reach + distinct local optima of the nonconvex objective (``0`` = zero start). + """ + h1 = np.asarray(integrals.h1, dtype=float) + g = np.asarray(integrals.g, dtype=float) + e_core = float(integrals.e_core) + n_elec = int(integrals.n_elec) + N = h1.shape[0] + + layout = _Layout(N=N, optimize_orbitals=optimize_orbitals) + off = ~np.eye(N, dtype=bool) + mask = off[:, None, :, None] & off[None, :, None, :] + + lambda_base = lambda_pauli(h1, g)["total"] + value, grad = _make_funcs(layout, h1, g, n_elec, mask) + + def fun(theta, mu): + return value(theta, mu), grad(theta, mu) + + if x0 is not None: + x = np.asarray(x0, dtype=float).copy() + if x.shape[0] != layout.n_par: + raise ValueError("x0 has the wrong length for the parameter layout") + elif rng_seed: + # A small random perturbation of the zero start; the scale is tiny so the + # first eval is still near the base 1-norm while distinct seeds reach + # distinct local optima of the nonconvex objective. + rng = np.random.default_rng(rng_seed) + x = 1e-2 * rng.standard_normal(layout.n_par) + else: + x = np.zeros(layout.n_par) + + bounds = None + if nu_nonneg: + bounds = [(None, None)] * layout.n_par + bounds[layout.nX + layout.nS] = (0.0, None) # nu slot + + res = None + for mu in mu_schedule: + res = minimize( + fun, x, args=(mu,), jac=True, method="L-BFGS-B", bounds=bounds, + options=dict(maxiter=maxiter, maxfun=10 * maxiter, + ftol=1e-15, gtol=1e-12, maxcor=30), + ) + x = res.x + + # The achieved 1-norm at the optimum. + hp, gp, U, X, s, nu, c = _pipeline(x, layout, h1, g, n_elec) + lambda_reduced = lambda_pauli(hp, gp)["total"] + + return ReducedIntegrals( + h1=hp, + g=gp, + e_core=e_core + shift_core_constant(nu, c, n_elec), + n_elec=n_elec, + lambda_reduced=float(lambda_reduced), + lambda_base=float(lambda_base), + U=U, + s=s, + nu=float(nu), + c=float(c), + ) diff --git a/prqpe/cgs/__init__.py b/prqpe/cgs/__init__.py new file mode 100644 index 0000000..e6f7945 --- /dev/null +++ b/prqpe/cgs/__init__.py @@ -0,0 +1,117 @@ +"""Ground-state Trotter constant ``C_gs``. + +The cost core needs a single number, the ground-state Trotter constant ``C_gs``, +that sets how fast the second-order product-formula bias shrinks with the step +size: ``|E0 - E0_tilde(delta)| = C_gs * delta^p`` with ``p`` the formula order +(``p = 2`` for the symmetric Trotter step). This module estimates that constant +two ways: + +* :func:`estimate_cgs` -- the direct fit on a small system. For a grid of step + sizes it builds the exact Trotter unitary, extracts the effective ground-state + energy, and fits the bias power law in log-log space. +* :func:`cgs_rule` -- the extrapolation ``C_gs = A * lambda^b`` for systems too + large to diagonalize, with ``(A, b)`` fit offline from :func:`estimate_cgs` + values across small systems. + +The fit helper :func:`fit_cgs` is exposed for callers that already hold the bias +errors. Everything here depends only on ``numpy``/``scipy`` and the +:class:`~prqpe.representation.pauli.PauliHamiltonian` container; the direct fit +is small-systems-only because it diagonalizes dense operators. +""" + +from __future__ import annotations + +import numpy as np + +from ..representation.pauli import PauliHamiltonian +from .trotter_matrix import ( + effective_ground_energy, + ground_state, + trot2_matrix, +) + +#: Default geometric grid of step sizes for :func:`estimate_cgs`. +DEFAULT_DELTAS = np.geomspace(0.2, 0.02, 7) + +__all__ = [ + "estimate_cgs", + "fit_cgs", + "cgs_rule", + "DEFAULT_DELTAS", +] + + +def fit_cgs(errors, deltas) -> tuple[float, float, np.ndarray]: + """Fit the bias power law ``error = C_gs * delta^p`` in log-log space. + + A straight-line least-squares fit of ``log(error)`` against ``log(delta)`` + gives the exponent ``p`` (the slope) and ``C_gs`` (the exponentiated + intercept). + + Parameters + ---------- + errors : array_like + The Trotter bias ``|E0 - E0_tilde(delta)|`` at each step size. All + entries must be strictly positive. + deltas : array_like + The corresponding step sizes. + + Returns + ------- + tuple + ``(C_gs, p, residuals)``: the fitted constant, the fitted exponent, and + the fit residuals ``log(error) - (p * log(delta) + log(C_gs))``. + """ + deltas = np.asarray(deltas, dtype=float) + errors = np.asarray(errors, dtype=float) + if np.any(errors <= 0.0): + raise ValueError( + "non-positive Trotter bias encountered; widen the step-size grid" + ) + log_d = np.log(deltas) + log_e = np.log(errors) + slope, intercept = np.polyfit(log_d, log_e, 1) + residuals = log_e - (slope * log_d + intercept) + return float(np.exp(intercept)), float(slope), residuals + + +def estimate_cgs(ham: PauliHamiltonian, deltas=None) -> float: + """Estimate the ground-state Trotter constant ``C_gs`` of a Hamiltonian. + + The exact ground state is found by diagonalization. For each step size the + symmetric Trotter unitary ``S_2(delta)`` is built, its effective + ground-state energy is extracted, and the bias ``|E0 - E0_tilde(delta)|`` is + fit to ``C_gs * delta^p`` in log-log space. The returned value is the + fitted ``C_gs``. + + Parameters + ---------- + ham : PauliHamiltonian + The Hamiltonian. + deltas : array_like, optional + The step-size grid. Defaults to :data:`DEFAULT_DELTAS`, a geometric + grid small enough to sit in the asymptotic regime. + + Returns + ------- + float + The fitted ground-state Trotter constant ``C_gs``. + """ + deltas = DEFAULT_DELTAS if deltas is None else np.asarray(deltas, dtype=float) + e0, psi0 = ground_state(ham) + errors = np.empty(len(deltas), dtype=float) + for i, delta in enumerate(deltas): + U = trot2_matrix(ham, delta) + e0_tilde, _ = effective_ground_energy(U, delta, psi0) + errors[i] = abs(e0 - e0_tilde) + c_gs, _p, _residuals = fit_cgs(errors, deltas) + return c_gs + + +def cgs_rule(lam: float, A: float, b: float) -> float: + """Extrapolate ``C_gs = A * lambda**b`` from a fitted scaling rule. + + ``(A, b)`` is fit offline from :func:`estimate_cgs` values on small systems; + ``lam`` is the Hamiltonian 1-norm ``lambda = sum_l |h_l|``. + """ + return float(A) * float(lam) ** float(b) diff --git a/prqpe/cgs/trotter_matrix.py b/prqpe/cgs/trotter_matrix.py new file mode 100644 index 0000000..c2a53de --- /dev/null +++ b/prqpe/cgs/trotter_matrix.py @@ -0,0 +1,107 @@ +"""Dense second-order Trotter unitaries and exact ground-state extraction. + +The ground-state Trotter constant is read off the bias of the symmetric +second-order product formula. This module assembles that formula as an explicit +matrix and pulls the effective ground-state energy out of it by diagonalization: + +* :func:`trot2_matrix` -- the dense symmetric Trotter unitary ``S_2(delta)``; +* :func:`effective_ground_energy` -- the energy and overlap of the Trotter + eigenstate closest to the exact ground state; +* :func:`ground_state` -- the exact ground state ``(E0, |psi0>)`` of the + Hamiltonian. + +Every routine here diagonalizes a dense ``dim x dim`` operator, so it is meant +for small systems (``dim`` up to a few thousand). It depends only on ``numpy``, +``scipy``, and the :class:`~prqpe.representation.pauli.PauliHamiltonian` +container. +""" + +from __future__ import annotations + +import numpy as np + +from ..representation.pauli import PauliHamiltonian + + +def trot2_matrix(ham: PauliHamiltonian, delta: float) -> np.ndarray: + """The symmetric second-order Trotter unitary ``S_2(delta)`` as a dense array. + + For ``H = sum_l h_l P_l`` the symmetric (Strang) splitting is + + S_2(delta) = (prod_l exp(-i delta H_l / 2)) + (prod_l reversed exp(-i delta H_l / 2)), + + a half step through the terms in stored order followed by a half step back. + The product is built by left-multiplying the identity with each rotation via + :meth:`PauliHamiltonian.rotation_apply_matrix`. + + Parameters + ---------- + ham : PauliHamiltonian + The Hamiltonian whose Trotter step is assembled. + delta : float + The Trotter step size. + + Returns + ------- + numpy.ndarray + The dense ``(dim, dim)`` unitary. Building it costs ``O(L * dim^2)`` + work. + """ + M = np.eye(ham.dim, dtype=complex) + half = 0.5 * float(delta) + # Rightmost factor first: ascending half-sweep, then descending half-sweep. + for l in range(ham.L): + M = ham.rotation_apply_matrix(l, half * ham.coeffs[l], M) + for l in range(ham.L - 1, -1, -1): + M = ham.rotation_apply_matrix(l, half * ham.coeffs[l], M) + return M + + +def effective_ground_energy(U: np.ndarray, delta: float, gs) -> tuple[float, float]: + """Effective ground-state energy of a Trotter unitary ``U = S_2(delta)``. + + Writing ``U = exp(-i delta H_eff)``, the eigenphases of ``U`` are + ``delta`` times the effective energies. The eigenvector with the largest + overlap with the exact ground state ``gs`` is taken to be the Trotter image + of the true ground state, and its energy is read from the eigenphase as + ``E0_tilde = -arg(eigenvalue) / delta``. + + Parameters + ---------- + U : numpy.ndarray + A dense ``(dim, dim)`` Trotter unitary. + delta : float + The Trotter step that produced ``U``. + gs : array_like + The exact ground-state vector, length ``dim``. + + Returns + ------- + tuple of float + ``(E0_tilde, overlap_sq)``: the effective ground-state energy and the + squared overlap of the selected eigenvector with ``gs``. + """ + evals, evecs = np.linalg.eig(U) + overlaps = np.abs(np.asarray(gs, dtype=complex).conj() @ evecs) ** 2 + j = int(np.argmax(overlaps)) + e0_tilde = -np.angle(evals[j]) / float(delta) + return float(e0_tilde), float(overlaps[j]) + + +def ground_state(ham: PauliHamiltonian) -> tuple[float, np.ndarray]: + """The exact ground state ``(E0, |psi0>)`` of a Hamiltonian. + + Parameters + ---------- + ham : PauliHamiltonian + The Hamiltonian to diagonalize. + + Returns + ------- + tuple + ``(E0, psi0)``: the ground-state energy and the corresponding + complex eigenvector. + """ + evals, evecs = np.linalg.eigh(ham.to_dense()) + return float(evals[0]), np.ascontiguousarray(evecs[:, 0]).astype(complex) diff --git a/prqpe/costs.py b/prqpe/costs.py new file mode 100644 index 0000000..3d29312 --- /dev/null +++ b/prqpe/costs.py @@ -0,0 +1,298 @@ +"""Closed-form cost primitives for (partially) randomized phase estimation. + +Every quantity in this module is an analytic resource formula for ground-state +energy estimation with second-order product formulas that are *deterministic*, +*randomized*, or a *partial* mix of the two. + +The cost model has three layers, each a small group of functions below: + +1. Product-formula error control -- splitting the precision budget between + phase estimation and Trotter error, and the resulting Trotter step size. +2. Per-operation gate costs -- Toffolis per Pauli rotation (Hamming-weight + phasing with a shared catalyst phase state) and per arbitrary-angle + rotation synthesis. +3. The partially randomized total -- a closed-form gate bound at a fixed + deterministic/randomized split, jointly optimized over the Trotter step + ``delta`` and the randomization scale ``kappa``. + +All formulas take the *unnormalized* Hamiltonian 1-norm ``lam`` and the target +precision ``eps`` in the same (Hartree) units. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass + +import numpy as np + +# ---------------------------------------------------------------------------- +# Empirical robust-phase-estimation sample-schedule constants. +# +# Robust phase estimation runs rounds m = 0, ..., M and takes +# N_m = N_M + D*(M - m) samples per quadrature in round m. The numerically +# calibrated values for an exact eigenstate are N_M = 11 in the deepest round +# and slope D = 4. These two constants set the closed-form prefactors of the +# total-gate bound below. +# ---------------------------------------------------------------------------- + +RPE_LAST_ROUND_SAMPLES = 11.0 # N_M: samples per quadrature in the deepest round +RPE_SLOPE = 4.0 # D: extra samples per quadrature per round below the deepest + +#: Empirical depth constant: the deepest round reaches t_max with +#: eps_qpe ~= C_MAX_TIME * pi / t_max, so M = ceil(log2(C_MAX_TIME*pi/(eps*delta))). +C_MAX_TIME = 0.08 + +#: Depth constant in the closed-form gate bound: 2^M = C_MAX_GATE*pi/(delta*eps_qpe). +C_MAX_GATE = 0.1 + +#: Two T gates are counted as one Toffoli gate. +T_GATES_PER_TOFFOLI = 2.0 + + +def _gate_bound_prefactors(N_M: float = RPE_LAST_ROUND_SAMPLES, + D: float = RPE_SLOPE) -> tuple[float, float]: + """Closed-form prefactors (deterministic, randomized) of the total bound. + + Summing the per-round sample schedule against the per-round gate counts + gives ``2*(N_M + D)`` for the deterministic term and ``8*(N_M/3 + D/9)`` + for the randomized term. For (N_M, D) = (11, 4) these are 30 and 296/9. + """ + return 2.0 * (N_M + D), 8.0 * (N_M / 3.0 + D / 9.0) + + +# ---------------------------------------------------------------------------- +# 1. Product-formula error control. +# ---------------------------------------------------------------------------- + +def qpe_trotter_split(eps: float, order: int = 2) -> tuple[float, float]: + """Split the precision budget as eps^2 = eps_qpe^2 + eps_trotter^2. + + The phase-estimation error is unbiased and independent of the Trotter + error, so the two add in quadrature. Minimizing the total gate count + under this constraint puts ``eps_qpe = eps*sqrt(p/(p+1))`` and hence + ``eps_trotter = eps/sqrt(p+1)``. + """ + eps_qpe = eps * math.sqrt(order / (order + 1.0)) + eps_trotter = eps / math.sqrt(order + 1.0) + return eps_qpe, eps_trotter + + +def trotter_step_size(eps: float, C_gs: float, order: int = 2) -> float: + """Cost-optimal Trotter step ``delta = (eps / (C_gs*sqrt(p+1)))^(1/p)``. + + The ground-state Trotter bias obeys ``|E0 - E0~| = C_gs * delta^p``, so the + Trotter share ``eps_trotter = eps/sqrt(p+1)`` of the budget fixes + ``delta^p = eps / (C_gs*sqrt(p+1))``. Equivalently + ``delta = (eps/C_gs)^(1/p) * (1/(p+1))^(1/(2p))``; this is the step that is + self-consistent with ``eps^2 = eps_qpe^2 + eps_trotter^2`` and with the + per-rotation synthesis allocation in :func:`synthesis_precision_per_rotation`. + """ + return (eps / (C_gs * math.sqrt(order + 1.0))) ** (1.0 / order) + + +# ---------------------------------------------------------------------------- +# 2. Per-operation gate costs. +# ---------------------------------------------------------------------------- + +def phase_state_precision(lam: float, eps: float) -> int: + """Dyadic-angle exponent ``J = ceil(log2(lam / (pi*eps)))`` (at least one bit). + + The randomized product formula is free to choose its step so that every + Pauli rotation is over the same angle ``pi*2^-J``, implemented from a shared + ``J``-qubit catalyst phase state. ``J`` is fixed by requiring the step to be + ``O(eps/lam)``. When the weight is below the precision scale (``lam < pi*eps``) + a single bit suffices and the randomized contribution is negligible. + """ + return max(1, int(math.ceil(math.log2(lam / (math.pi * eps))))) + + +def rotation_toffolis(J: int, K: int) -> float: + """Toffolis per Pauli rotation with Hamming-weight phasing: ``1 + (J-2)/K``. + + A run of ``K`` mutually commuting Paulis is mapped to single-qubit Z + rotations and phased through the catalyst state: ``K-1`` Toffolis compute and + uncompute the Hamming weight and ``J-1`` Toffolis apply the phase, i.e. + ``((K-1)+(J-1))/K = 1 + (J-2)/K`` Toffolis per rotation. + """ + return 1.0 + (J - 2.0) / K + + +def phasing_ancillas(J: int, K: int) -> int: + """Ancilla qubits for Hamming-weight phasing: ``K + 2J - 2``. + + A ``J``-qubit catalyst phase state, ``K-1`` Hamming-weight register qubits, + and ``J-1`` adder ancillas: ``J + (K-1) + (J-1) = K + 2J - 2``. + """ + return K + 2 * J - 2 + + +def logical_qubits(n_system: int, J: int, K: int) -> int: + """Total logical qubits: system + one phase-estimation ancilla + phasing. + + ``n_system + 1 + (K + 2J - 2)``. The single extra ancilla carries the + Hadamard test; the remaining ``K + 2J - 2`` are the Hamming-weight phasing + register, adder ancillas, and catalyst phase state of :func:`phasing_ancillas`. + """ + return n_system + 1 + phasing_ancillas(J, K) + + +def synthesis_toffolis(eps_prime: float) -> float: + """Toffolis to synthesize one arbitrary-angle rotation to error ``eps'``. + + A rotation costs ``1.14*log2(1/eps') + 9.2`` T gates, i.e. half that many + Toffolis at two T per Toffoli. + """ + t_gates = 1.14 * math.log2(1.0 / eps_prime) + 9.2 + return t_gates / T_GATES_PER_TOFFOLI + + +def synthesis_precision_per_rotation(delta: float, eps_synth: float, + n_stage: int, n_rot: int) -> float: + """Per-rotation synthesis error ``eps' = delta*eps_synth/(n_stage*n_rot)``. + + Imperfect synthesis of every rotation perturbs the effective Hamiltonian by + at most ``n_stage*n_rot*eps'/delta``; keeping the induced ground-state error + below ``eps_synth`` fixes ``eps'`` as above. ``n_rot`` is the number of + rotations per Trotter stage (the deterministic term count for a Pauli + Hamiltonian). + """ + return delta * eps_synth / (n_stage * n_rot) + + +# ---------------------------------------------------------------------------- +# 3. The partially randomized total. +# ---------------------------------------------------------------------------- + +@dataclass(frozen=True) +class SplitOptimum: + """Cost-optimal operating point at a fixed deterministic/randomized split.""" + + toffoli_total: float + delta: float + kappa: float + eps_qpe: float + + +def optimal_kappa(A: float, B: float) -> float: + """Minimizer of ``G(kappa) = (A + B*kappa) * e^(2/kappa)`` over ``kappa > 0``. + + ``kappa`` scales the randomized rotation budget ``r = kappa*lam_R^2*delta^2*4^m`` + and sets the damping ``B_damp <= e^(1/kappa)`` that inflates every sample + count by ``e^(2/kappa)``. Setting ``dG/dkappa = 0`` gives the quadratic + ``B*kappa^2 - 2*B*kappa - 2*A = 0`` with positive root + ``kappa* = 1 + sqrt(1 + 2A/B)``. ``A`` is the deterministic prefactor + (no randomized part -> ``kappa* -> inf``); ``A = 0`` recovers the fully + randomized optimum ``kappa* = 2``. + """ + if B == 0.0: + return math.inf + return 1.0 + math.sqrt(1.0 + 2.0 * A / B) + + +def partial_toffoli_bound(L_D: int, lam_R: float, delta: float, eps_qpe: float, + kappa: float, G_det: float, G_rand: float, + n_stage: int = 2, randomizer_factor: float = 1.0, + N_M: float = RPE_LAST_ROUND_SAMPLES, + D: float = RPE_SLOPE) -> float: + """Closed-form total Toffoli bound at a fixed split and operating point. + + Summing ``sum_m 2*N_m*(G_det*n_stage*L_D*2^(m-1) + G_rand*kappa*lam_R^2*delta^2*4^m)`` + over the robust-phase-estimation rounds, with ``N_m = e^(2/kappa)*(N_M + D*(M-m))`` + and ``2^M = C_MAX_GATE*pi/(delta*eps_qpe)``, gives + + G = 2*(N_M+D) * G_det*n_stage*L_D * e^(2/kappa) * C_MAX_GATE*pi/(delta*eps_qpe) + + 8*(N_M/3+D/9) * G_rand*kappa*e^(2/kappa) * (C_MAX_GATE*pi*lam_R)^2/eps_qpe^2. + + The deterministic term already includes the factor-2 reduction from halving + the deterministic evolution time (``2^(m-1)`` rather than ``2^m``). ``L_D = 0`` + drops the deterministic term; ``lam_R = 0`` (``kappa = inf``) drops the + randomized term. ``randomizer_factor`` scales the randomized rotation count + (``1`` for RTE, ``1/2`` for qDRIFT). + """ + det_pre, rand_pre = _gate_bound_prefactors(N_M, D) + inflation = 1.0 if math.isinf(kappa) else math.exp(2.0 / kappa) + + det = 0.0 + if L_D > 0 and delta > 0.0: + det = (det_pre * G_det * n_stage * L_D * inflation + * C_MAX_GATE * math.pi / (delta * eps_qpe)) + rand = 0.0 + if lam_R > 0.0: + rand = (rand_pre * randomizer_factor * G_rand * kappa * inflation + * (C_MAX_GATE * math.pi * lam_R) ** 2 / eps_qpe ** 2) + return det + rand + + +def minimize_delta_kappa(L_D: int, lam_R: float, eps: float, C_gs: float, + G_det: float, G_rand: float, *, order: int = 2, + n_stage: int = 2, randomizer_factor: float = 1.0, + n_grid: int = 400, + N_M: float = RPE_LAST_ROUND_SAMPLES, + D: float = RPE_SLOPE) -> SplitOptimum: + """Jointly optimize the total bound over ``(delta, kappa)`` at fixed ``L_D``. + + Substituting ``eps_qpe = sqrt(eps^2 - C_gs^2*delta^4)`` couples the Trotter + step to the phase-estimation budget. Parameterize ``x = C_gs*delta^2/eps`` + in ``(0, 1)``, so ``delta = sqrt(x*eps/C_gs)`` and ``eps_qpe = eps*sqrt(1-x^2)``; + at each grid point ``kappa`` is the closed-form root of :func:`optimal_kappa`, + and the global minimum over the grid is returned. + + The two degenerate splits are exact closed forms (no grid): + + * ``L_D = 0`` (fully randomized): ``kappa = 2``, ``eps_qpe = eps``, + ``delta`` immaterial. + * ``lam_R = 0`` (fully deterministic): ``kappa = inf`` and the optimum sits + at ``x = 1/sqrt(3)``, i.e. ``eps_qpe = eps*sqrt(2/3)``. + """ + if order != 2: + raise NotImplementedError( + "the partially randomized cost model is implemented for order 2") + + if L_D == 0: + total = partial_toffoli_bound(0, lam_R, 0.0, eps, 2.0, G_det, G_rand, + n_stage=n_stage, + randomizer_factor=randomizer_factor, + N_M=N_M, D=D) + return SplitOptimum(total, delta=0.0, kappa=2.0, eps_qpe=eps) + + if lam_R == 0.0: + x = 1.0 / math.sqrt(3.0) + delta = math.sqrt(x * eps / C_gs) + eps_qpe = eps * math.sqrt(1.0 - x ** 2) + total = partial_toffoli_bound(L_D, 0.0, delta, eps_qpe, math.inf, + G_det, G_rand, n_stage=n_stage, + randomizer_factor=randomizer_factor, + N_M=N_M, D=D) + return SplitOptimum(total, delta=delta, kappa=math.inf, eps_qpe=eps_qpe) + + # Grid in x, denser near 0, with the analytic fully-deterministic point + # x = 1/sqrt(3) included so the interior optimum is never missed. + x = np.linspace(1e-4, 1.0 - 1e-6, n_grid) ** 2 + x = np.sort(np.append(x, 1.0 / math.sqrt(3.0))) + delta = np.sqrt(x * eps / C_gs) + eps_qpe = eps * np.sqrt(1.0 - x ** 2) + + det_pre, rand_pre = _gate_bound_prefactors(N_M, D) + A = det_pre * G_det * n_stage * L_D * C_MAX_GATE * math.pi / (delta * eps_qpe) + B = (rand_pre * randomizer_factor * G_rand + * (C_MAX_GATE * math.pi * lam_R) ** 2 / eps_qpe ** 2) + kappa = 1.0 + np.sqrt(1.0 + 2.0 * A / B) + total = (A + B * kappa) * np.exp(2.0 / kappa) + + total = np.where(np.isfinite(total), total, np.inf) + i = int(np.argmin(total)) + return SplitOptimum(float(total[i]), float(delta[i]), float(kappa[i]), + float(eps_qpe[i])) + + +def randomized_deepest_circuit_toffolis(total_randomized_toffolis: float, + N_M: float = RPE_LAST_ROUND_SAMPLES, + D: float = RPE_SLOPE) -> float: + """Toffolis in the deepest randomized circuit, from the randomized total. + + The schedule sums to ``r_tot = 8*(N_M/3 + D/9) * r_max``, so the deepest + single circuit carries ``total / (8*(N_M/3 + D/9))`` of the Toffolis. + """ + _, rand_pre = _gate_bound_prefactors(N_M, D) + return total_randomized_toffolis / rand_pre diff --git a/prqpe/representation/__init__.py b/prqpe/representation/__init__.py new file mode 100644 index 0000000..b64c8c2 --- /dev/null +++ b/prqpe/representation/__init__.py @@ -0,0 +1,132 @@ +"""Hamiltonian-to-cost-input adapter. + +The cost core consumes only the descending-sorted weights ``|h_l|``, the system +qubit count, and the ground-state Trotter constant. This module owns the step that +turns a qubit Hamiltonian into the operator container that produces those weights: + +* :func:`build_representation` -- convert a qubit Hamiltonian into a + :class:`~prqpe.representation.pauli.PauliHamiltonian`. +* :func:`truncate_small_terms` -- drop the smallest-weight terms within a fixed + weight budget. +* :func:`sorted_weights` -- the descending ``|h_l|`` array the cost core reads. + +Only :func:`build_representation` touches an external qubit-Hamiltonian package, +and it imports it lazily. +""" + +from __future__ import annotations + +import numpy as np + +from .pauli import PauliHamiltonian + +#: Default cumulative weight budget for :func:`truncate_small_terms`. +DEFAULT_WEIGHT_THRESHOLD = 1e-3 + +__all__ = [ + "PauliHamiltonian", + "build_representation", + "truncate_small_terms", + "sorted_weights", + "DEFAULT_WEIGHT_THRESHOLD", +] + + +def build_representation(hamiltonian, representation: str = "pauli" + ) -> PauliHamiltonian: + """Build a :class:`PauliHamiltonian` from a qubit Hamiltonian. + + Parameters + ---------- + hamiltonian : PauliHamiltonian or qubit operator + Either an already-built :class:`PauliHamiltonian` (returned unchanged) or + an openfermion ``QubitOperator``. + representation : str + Only ``"pauli"`` is implemented. + + Returns + ------- + PauliHamiltonian + The Pauli-string Hamiltonian, with any identity term dropped (a constant + energy shift carries no resource cost and is reinstated classically). + """ + if representation != "pauli": + raise ValueError( + f"unknown representation {representation!r}; only 'pauli' is implemented") + + if isinstance(hamiltonian, PauliHamiltonian): + return hamiltonian + return _pauli_from_qubit_operator(hamiltonian) + + +def _pauli_from_qubit_operator(qop, tol: float = 1e-12) -> PauliHamiltonian: + """Convert an openfermion ``QubitOperator`` to a :class:`PauliHamiltonian`. + + The letter-to-bitmask map follows :mod:`prqpe.representation.pauli`. The + identity term is dropped and terms with negligible real part are skipped. + """ + from openfermion import count_qubits + + n_qubits = count_qubits(qop) + xs, zs, cs = [], [], [] + for term, coeff in qop.terms.items(): + coeff = complex(coeff) + if abs(coeff.imag) > 1e-9: + raise ValueError(f"non-real Pauli coefficient {coeff} for {term}") + if not term: # identity term -> classical constant, no cost + continue + if abs(coeff.real) < tol: + continue + x = z = 0 + for q, letter in term: + bit = 1 << (n_qubits - 1 - q) + if letter in ("X", "Y"): + x |= bit + if letter in ("Z", "Y"): + z |= bit + xs.append(x) + zs.append(z) + cs.append(coeff.real) + return PauliHamiltonian(n_qubits, xs, zs, cs) + + +def truncate_small_terms(ham: PauliHamiltonian, + weight_threshold: float = DEFAULT_WEIGHT_THRESHOLD + ) -> PauliHamiltonian: + """Drop the smallest-weight terms within a cumulative weight budget. + + Working upward from the smallest ``|h_l|``, terms are dropped while the + *cumulative* dropped weight stays within ``weight_threshold``. The first + term that would push the cumulative dropped weight past the budget is kept, + along with everything larger; the kept Hamiltonian therefore differs from the + input by at most ``weight_threshold`` in 1-norm. + + Parameters + ---------- + ham : PauliHamiltonian + The Hamiltonian to truncate. + weight_threshold : float + Maximum total dropped weight (same units as the coefficients). + + Returns + ------- + PauliHamiltonian + A new Hamiltonian with the kept terms, ordered by decreasing ``|h_l|``. + """ + ordered = ham.sorted_by_weight() + weights = np.abs(ordered.coeffs) + # Cumulative weight accumulated from the smallest term upward. + cum_from_small = np.cumsum(weights[::-1]) + # Number of trailing (smallest) terms whose cumulative weight stays within + # the budget: searchsorted from the right finds the largest such count. + n_drop = int(np.searchsorted(cum_from_small, weight_threshold, side="right")) + n_keep = ordered.L - n_drop + return ordered.subset(np.arange(n_keep)) + + +def sorted_weights(ham: PauliHamiltonian) -> np.ndarray: + """The descending ``|h_l|`` array the cost core consumes. + + Equivalent to ``np.abs(ham.sorted_by_weight().coeffs)``. + """ + return np.abs(ham.sorted_by_weight().coeffs) diff --git a/prqpe/representation/pauli.py b/prqpe/representation/pauli.py new file mode 100644 index 0000000..093cf1e --- /dev/null +++ b/prqpe/representation/pauli.py @@ -0,0 +1,154 @@ +"""Real-coefficient Pauli-string Hamiltonian container. + +A Hamiltonian is stored as a weighted sum of Pauli strings, + + H = sum_l h_l P_l, h_l real, P_l a tensor product of {I, X, Y, Z}, + +with 1-norm ``lambda = sum_l |h_l|``. + +Each Pauli string is held as a pair of integer bitmasks ``(x, z)`` with the phase +convention + + P(x, z) = i^popcount(x AND z) * X^x Z^z + +(the X factor applied after the Z factor). This is Hermitian and self-inverse, and +maps the letter ``Y`` on a qubit to ``x_q = z_q = 1`` since ``Y = i X Z``. Bit +``n_qubits - 1 - q`` of a computational-basis index belongs to qubit ``q`` -- qubit 0 +is the most significant bit -- so the dense and sparse matrices built here agree +entrywise with the standard sparse-operator ordering used by qubit-Hamiltonian +toolkits. + +Action on a basis state ``|b>``:: + + P(x, z) |b> = i^popcount(x & z) * (-1)^popcount(b & z) |b XOR x>, + +hence ``(P v)[c] = i^popcount(x & z) * (-1)^popcount((c ^ x) & z) * v[c ^ x]``: a bit +permutation followed by a diagonal phase, vectorized over terms below. + +Pauli rotations use ``P^2 = I``:: + + exp(-i phi P) |v> = cos(phi) |v> - i sin(phi) P |v>. + +This module depends only on ``numpy`` and ``scipy``. +""" + +from __future__ import annotations + +import numpy as np +import scipy.sparse as sp + + +def popcount(a): + """Number of set bits, elementwise, for a nonnegative integer array.""" + a = np.asarray(a, dtype=np.uint64) + if hasattr(np, "bitwise_count"): + return np.bitwise_count(a).astype(np.int64) + out = np.zeros(a.shape, dtype=np.int64) + x = a.copy() + while np.any(x): + out += (x & np.uint64(1)).astype(np.int64) + x >>= np.uint64(1) + return out + + +class PauliHamiltonian: + """A real-coefficient Pauli-string Hamiltonian ``H = sum_l h_l P_l``. + + Parameters + ---------- + n_qubits : int + Number of qubits the Hamiltonian acts on. + x_masks, z_masks : integer arrays of shape ``(L,)`` + The X and Z bitmasks of each Pauli string, in the convention of the + module docstring. + coeffs : float array of shape ``(L,)`` + The real coefficients ``h_l``. + """ + + def __init__(self, n_qubits, x_masks, z_masks, coeffs): + self.n_qubits = int(n_qubits) + self.dim = 1 << self.n_qubits + self.x_masks = np.asarray(x_masks, dtype=np.int64).ravel() + self.z_masks = np.asarray(z_masks, dtype=np.int64).ravel() + coeffs = np.asarray(coeffs) + if np.iscomplexobj(coeffs): + if np.max(np.abs(coeffs.imag), initial=0.0) > 1e-10: + raise ValueError("Pauli coefficients must be real") + coeffs = coeffs.real + self.coeffs = np.asarray(coeffs, dtype=np.float64).ravel() + if not (len(self.x_masks) == len(self.z_masks) == len(self.coeffs)): + raise ValueError("mask/coefficient length mismatch") + self._perms = None + self._phases = None + + # ------------------------------------------------------------------ basic + @property + def L(self): + """Number of Pauli terms.""" + return len(self.coeffs) + + @property + def lambda_weight(self): + """The 1-norm ``lambda = sum_l |h_l|``.""" + return float(np.abs(self.coeffs).sum()) + + def subset(self, idx): + """A new Hamiltonian holding the terms ``idx`` (order preserved).""" + idx = np.asarray(idx) + return PauliHamiltonian( + self.n_qubits, self.x_masks[idx], self.z_masks[idx], self.coeffs[idx] + ) + + def sorted_by_weight(self): + """A copy with terms ordered by decreasing ``|h_l|`` (stable sort).""" + order = np.argsort(-np.abs(self.coeffs), kind="stable") + return self.subset(order) + + # ------------------------------------------------------------------ tables + def _tables(self): + """Per-term permutation and diagonal phase. + + Returns ``(perms, phases)`` of shape ``(L, dim)`` such that + ``(P_l v)[c] = phases[l, c] * v[perms[l, c]]``: ``perms`` is the + column permutation ``c XOR x_l`` and ``phases`` carries the head phase + ``i^popcount(x_l & z_l)`` times the basis sign ``(-1)^popcount(c & z_l)``. + Cached after the first call. + """ + if self._perms is None: + c = np.arange(self.dim, dtype=np.int64) + self._perms = c[None, :] ^ self.x_masks[:, None] + head = 1j ** (popcount(self.x_masks & self.z_masks) % 4) + signs = 1.0 - 2.0 * (popcount(self._perms & self.z_masks[:, None]) % 2) + self._phases = head[:, None] * signs + return self._perms, self._phases + + # -------------------------------------------------------- matrix operators + def pauli_apply_matrix(self, l, M): + """Left-multiply a dense ``(dim, dim)`` matrix: ``P_l @ M``.""" + perms, phases = self._tables() + return phases[l][:, None] * M[perms[l], :] + + def rotation_apply_matrix(self, l, phi, M): + """``exp(-i phi P_l) @ M`` for dense ``M`` (assembles product-formula unitaries).""" + return np.cos(phi) * M - 1j * np.sin(phi) * self.pauli_apply_matrix(l, M) + + def term_matrix(self, l): + """The sparse matrix of ``P_l``.""" + perms, phases = self._tables() + data = phases[l].copy() + rows = np.arange(self.dim) + return sp.csr_matrix((data, (rows, perms[l])), shape=(self.dim, self.dim)) + + def to_sparse(self): + """``H = sum_l h_l P_l`` as a scipy CSR matrix (small ``n_qubits`` only).""" + perms, phases = self._tables() + data = (self.coeffs[:, None] * phases).ravel() + rows = np.tile(np.arange(self.dim), self.L) + cols = perms.ravel() + return sp.coo_matrix( + (data, (rows, cols)), shape=(self.dim, self.dim) + ).tocsr() + + def to_dense(self): + """``H`` as a dense complex array (small ``n_qubits`` only).""" + return self.to_sparse().toarray() diff --git a/prqpe/schedule.py b/prqpe/schedule.py new file mode 100644 index 0000000..d062101 --- /dev/null +++ b/prqpe/schedule.py @@ -0,0 +1,129 @@ +"""Robust-phase-estimation round and sample schedule. + +Phase estimation runs rounds ``m = 0, ..., M`` at evolution time ``t = 2^m``, +taking ``N_m`` samples per quadrature in round ``m``. This module owns the +round count ``M``, the sample schedule ``N_m = N_M + D*(M - m)`` (inflated by the +randomized-formula damping), the resulting circuit count, and the per-circuit +Toffoli cost from which the deepest circuit is read off. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass + +from .costs import C_MAX_TIME, RPE_LAST_ROUND_SAMPLES, RPE_SLOPE + +#: Lowest ground-state overlap for which robust phase estimation is guaranteed. +MIN_OVERLAP = 4.0 - 2.0 * math.sqrt(3.0) # ~= 0.536 + + +def resolve_depth_tradeoff(overlap: float = 1.0, xi: float | None = None) -> float: + """Return a valid depth-reduction factor ``xi`` in ``(0, 1]``. + + A smaller ``xi`` shortens the deepest circuit (``t_max ~ xi/eps``) at the + cost of more circuits. Robust phase estimation admits depth reduction only + when the guiding-state overlap is high: ``xi`` must exceed + ``arcsin((1-overlap)/overlap)``, which is ``0`` for an exact eigenstate + (``overlap = 1``). ``xi = None`` means no reduction and returns ``1.0``. + """ + if xi is None: + return 1.0 + if not 0.0 < xi <= 1.0: + raise ValueError(f"xi must lie in (0, 1], got {xi}") + if overlap <= MIN_OVERLAP: + raise ValueError( + f"overlap must exceed {MIN_OVERLAP:.3f} for robust phase estimation") + if overlap < 1.0: + xi_min = math.asin((1.0 - overlap) / overlap) + if xi <= xi_min: + raise ValueError( + f"xi must exceed arcsin((1-overlap)/overlap) = {xi_min:.4f}") + return xi + + +def _overlap_inflation(overlap: float) -> float: + """Sample-count inflation from imperfect overlap, normalized to ``1`` at ``eta=1``. + + The per-round sample count scales as ``beta^-2`` with the signal headroom + ``beta = overlap*(1 + sin(pi/3)) - 1``. The empirical schedule constants are + calibrated at ``overlap = 1`` (``beta = sin(pi/3)``), so the relative + inflation for general overlap is ``(beta(1)/beta(overlap))^2``. + """ + if overlap >= 1.0: + return 1.0 + if overlap <= MIN_OVERLAP: + raise ValueError( + f"overlap must exceed {MIN_OVERLAP:.3f} for robust phase estimation") + headroom = overlap * (1.0 + math.sin(math.pi / 3.0)) - 1.0 + return (math.sin(math.pi / 3.0) / headroom) ** 2 + + +def num_rounds(eps_qpe: float, delta: float, xi: float = 1.0) -> int: + """Round count for a Trotter-stepped arm: ``M = ceil(log2(C*pi*xi/(eps_qpe*delta)))``. + + The deepest round reaches ``t_max = 2^M`` with ``eps_qpe ~= C_MAX_TIME*pi/t_max`` + in units of the Trotter step ``delta``; the depth-reduction factor ``xi`` + shrinks ``M`` when the guiding-state overlap is high. + """ + return int(math.ceil(math.log2(C_MAX_TIME * math.pi * xi / (eps_qpe * delta)))) + + +def num_rounds_randomized(lam_R: float, eps_qpe: float, xi: float = 1.0) -> int: + """Round count for the fully randomized arm: ``M = ceil(log2(xi*lam_R/eps_qpe))``. + + With no Trotter step the deepest round estimates the normalized Hamiltonian + ``H/lam_R`` to precision ``eps_qpe/lam_R`` at time ``t_max = 2^M``. + """ + return int(math.ceil(math.log2(xi * lam_R / eps_qpe))) + + +@dataclass(frozen=True) +class RPESchedule: + """Sample schedule ``N_m = ceil(inflation*(N_M + D*(M - m)))`` over rounds ``0..M``.""" + + M: int + inflation: float + N_M: float = RPE_LAST_ROUND_SAMPLES + D: float = RPE_SLOPE + + def samples(self, m: int) -> int: + """Samples per quadrature in round ``m`` (decreasing toward the deepest).""" + return int(math.ceil(self.inflation * (self.N_M + self.D * (self.M - m)))) + + @property + def num_circuits(self) -> int: + """Total Hadamard-test circuits ``sum_m 2*N_m`` (two quadratures per round).""" + return int(2 * sum(self.samples(m) for m in range(self.M + 1))) + + +def rpe_schedule(M: int, kappa: float, overlap: float = 1.0) -> RPESchedule: + """Build the sample schedule for ``M`` rounds at randomization scale ``kappa``. + + The randomized product formula damps the signal by ``B <= e^(1/kappa)``, so + every sample count is inflated by ``B^2 = e^(2/kappa)``; a fully + deterministic arm (``kappa = inf``) has no inflation. Imperfect overlap + adds a further ``(beta(1)/beta(overlap))^2`` factor. + """ + damping = 1.0 if math.isinf(kappa) else math.exp(2.0 / kappa) + return RPESchedule(M=M, inflation=damping * _overlap_inflation(overlap)) + + +def round_circuit_toffolis(m: int, *, G_det: float, G_rand: float, L_D: int, + lam_R: float, delta: float, kappa: float, + n_stage: int = 2) -> float: + """Toffolis in the single depth-``m`` circuit (no sample factor). + + The deterministic part applies ``n_stage*L_D`` rotations over ``2^(m-1)`` + halved Trotter steps; the randomized part draws ``kappa*lam_R^2*delta^2*4^m`` + Pauli rotations. Each kind costs ``G_det`` and ``G_rand`` Toffolis. + """ + det = G_det * n_stage * L_D * 2.0 ** (m - 1) if L_D > 0 else 0.0 + rand = (G_rand * kappa * lam_R ** 2 * delta ** 2 * 4.0 ** m + if lam_R > 0.0 and delta > 0.0 else 0.0) + return det + rand + + +def max_toffolis_per_circuit(schedule: RPESchedule, **circuit) -> float: + """Toffolis in the deepest circuit, round ``m = M``.""" + return round_circuit_toffolis(schedule.M, **circuit) diff --git a/prqpe/tests/__init__.py b/prqpe/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/prqpe/tests/test_bliss.py b/prqpe/tests/test_bliss.py new file mode 100644 index 0000000..74f0434 --- /dev/null +++ b/prqpe/tests/test_bliss.py @@ -0,0 +1,278 @@ +"""First-principles tests for the Hamiltonian-weight reduction. + +The two load-bearing properties of the reduction are checked from scratch on +small synthetic integrals built with the full 8-fold permutation symmetry of +real chemists'-notation two-electron integrals: + +* the closed-form Pauli 1-norm strictly decreases (a reducible case), and +* the energy on the target electron-number sector is preserved to machine + precision (the block-invariant shift is zero on that sector). + +Energy preservation is verified by assembling the second-quantized Hamiltonian +as a dense Jordan-Wigner matrix in pure ``numpy`` and comparing the +``n``-electron block of the original and reduced operators. Everything here +depends only on ``numpy`` and ``scipy``. +""" + +from collections import namedtuple + +import numpy as np +import scipy.linalg as sla + +from qhat.prqpe.bliss import ( + ReducedIntegrals, + lambda_pauli, + reduce_hamiltonian_weight, + rotate_integrals, + shifted_integrals, + shift_core_constant, +) + +Integrals = namedtuple("Integrals", "h1 g e_core n_elec") + + +# ---------------------------------------------------------------------------- +# Synthetic integral construction. +# ---------------------------------------------------------------------------- + +def synthetic_integrals(N=4, seed=0, e_core=0.7, n_elec=4): + """Random symmetric ``h1`` and 8-fold-symmetric chemists' ``g``.""" + rng = np.random.default_rng(seed) + h1 = rng.standard_normal((N, N)) + h1 = (h1 + h1.T) / 2.0 + g = rng.standard_normal((N, N, N, N)) + g = sum( + g.transpose(t) for t in [ + (0, 1, 2, 3), (1, 0, 2, 3), (0, 1, 3, 2), (1, 0, 3, 2), + (2, 3, 0, 1), (3, 2, 0, 1), (2, 3, 1, 0), (3, 2, 1, 0), + ] + ) / 8.0 + return Integrals(h1, g, e_core, n_elec) + + +def test_g_has_eightfold_symmetry(): + """The synthetic ``g`` carries the permutation symmetries it should.""" + _, g, _, _ = synthetic_integrals() + assert np.allclose(g, g.transpose(1, 0, 2, 3)) # (pq|rs) = (qp|rs) + assert np.allclose(g, g.transpose(0, 1, 3, 2)) # (pq|rs) = (pq|sr) + assert np.allclose(g, g.transpose(2, 3, 0, 1)) # (pq|rs) = (rs|pq) + # the antisymmetrization index pattern is NOT a symmetry of g: + assert not np.allclose(g, g.transpose(0, 3, 2, 1)) + + +# ---------------------------------------------------------------------------- +# Pure-numpy Jordan-Wigner machinery for the energy-preservation test. +# ---------------------------------------------------------------------------- + +_I2 = np.eye(2) +_Z = np.array([[1.0, 0.0], [0.0, -1.0]]) +# Annihilation operator on one qubit (JW lowering): |1><0|. +_A = np.array([[0.0, 1.0], [0.0, 0.0]]) + + +def _kron_list(mats): + out = np.array([[1.0]]) + for m in mats: + out = np.kron(out, m) + return out + + +def jw_annihilation(p, n_modes): + """Jordan-Wigner annihilation operator on mode ``p`` of ``n_modes``.""" + mats = [_Z] * p + [_A] + [_I2] * (n_modes - p - 1) + return _kron_list(mats) + + +def to_dense(hh, gg, cc, n_spatial, n_alpha_beta=2): + """Dense Jordan-Wigner matrix of ``cc*I + one-body + two-body``. + + Spin orbitals are ordered ``(p, sigma)`` with ``sigma`` the fast index, i.e. + mode ``2*p + sigma``. The two-body term uses the chemists'-notation + contraction with the standard ``0.5`` prefactor and operator ordering + ``a^dag_{p s1} a^dag_{r s2} a_{s s2} a_{q s1}``. + """ + n_modes = n_spatial * n_alpha_beta + dim = 1 << n_modes + a = [jw_annihilation(m, n_modes) for m in range(n_modes)] + ad = [op.conj().T for op in a] + + H = cc * np.eye(dim) + for p in range(n_spatial): + for q in range(n_spatial): + if hh[p, q] == 0.0: + continue + for sig in range(n_alpha_beta): + H = H + hh[p, q] * (ad[2 * p + sig] @ a[2 * q + sig]) + + for p in range(n_spatial): + for q in range(n_spatial): + for r in range(n_spatial): + for s in range(n_spatial): + val = gg[p, q, r, s] + if val == 0.0: + continue + for s1 in range(n_alpha_beta): + for s2 in range(n_alpha_beta): + op = (ad[2 * p + s1] @ ad[2 * r + s2] + @ a[2 * s + s2] @ a[2 * q + s1]) + H = H + 0.5 * val * op + return H + + +def number_diag(n_spatial, n_alpha_beta=2): + """Diagonal of the total number operator over the JW computational basis.""" + n_modes = n_spatial * n_alpha_beta + dim = 1 << n_modes + idx = np.arange(dim) + counts = np.zeros(dim, dtype=int) + for m in range(n_modes): + bit = (idx >> (n_modes - 1 - m)) & 1 + counts += bit + return counts + + +# ---------------------------------------------------------------------------- +# (a) The 1-norm strictly decreases. +# ---------------------------------------------------------------------------- + +def test_lambda_decreases_shift_only(): + """Block-invariant shift alone strictly shrinks the 1-norm.""" + ints = synthetic_integrals(seed=0) + lam0 = lambda_pauli(ints.h1, ints.g)["total"] + res = reduce_hamiltonian_weight(ints, optimize_orbitals=False) + assert res.lambda_reduced < lam0 - 1e-3 + assert res.lambda_base == lam0 + + +def test_lambda_decreases_joint(): + """Joint orbital rotation + shift shrinks the 1-norm at least as much.""" + ints = synthetic_integrals(seed=0) + lam0 = lambda_pauli(ints.h1, ints.g)["total"] + shift_only = reduce_hamiltonian_weight(ints, optimize_orbitals=False) + joint = reduce_hamiltonian_weight(ints, optimize_orbitals=True) + assert joint.lambda_reduced < lam0 - 1e-3 + assert joint.lambda_reduced <= shift_only.lambda_reduced + 1e-6 + + +def test_zero_theta_matches_base_exactly(): + """At ``theta = 0`` (no rotation, ``s = nu = c = 0``) lambda equals the base.""" + ints = synthetic_integrals(seed=2) + hp, gp = shifted_integrals(ints.h1, ints.g, np.zeros((4, 4)), 0.0, 0.0, + ints.n_elec) + assert np.allclose(hp, ints.h1) + assert np.allclose(gp, ints.g) + assert lambda_pauli(hp, gp)["total"] == lambda_pauli(ints.h1, ints.g)["total"] + + +# ---------------------------------------------------------------------------- +# (b) Energy on the n-electron sector is preserved. +# ---------------------------------------------------------------------------- + +def _sector_block(H, n_diag, n_elec): + sector = np.isclose(n_diag, n_elec) + return H[np.ix_(sector, sector)] + + +def test_shift_preserves_sector_with_random_params(): + """A random shift leaves the whole ``n``-electron block unchanged.""" + N = 3 + ints = synthetic_integrals(N=N, seed=4, e_core=0.3, n_elec=N) + rng = np.random.default_rng(7) + s = rng.standard_normal((N, N)) + s = (s + s.T) / 2.0 + nu = float(rng.standard_normal()) + c = float(rng.standard_normal()) + + hp, gp = shifted_integrals(ints.h1, ints.g, s, nu, c, ints.n_elec) + e_core_p = ints.e_core + shift_core_constant(nu, c, ints.n_elec) + + H0 = to_dense(ints.h1, ints.g, ints.e_core, N) + H1 = to_dense(hp, gp, e_core_p, N) + n_diag = number_diag(N) + + blk0 = _sector_block(H0, n_diag, ints.n_elec) + blk1 = _sector_block(H1, n_diag, ints.n_elec) + assert np.abs(blk0 - blk1).max() < 1e-9 + e0 = np.linalg.eigvalsh(blk0)[0] + e1 = np.linalg.eigvalsh(blk1)[0] + assert abs(e0 - e1) < 1e-9 + + +def test_rotation_preserves_full_spectrum(): + """An orbital rotation is a basis change and preserves the whole spectrum.""" + N = 3 + ints = synthetic_integrals(N=N, seed=5, e_core=0.0, n_elec=N) + rng = np.random.default_rng(9) + X = rng.standard_normal((N, N)) + U = sla.expm(X - X.T) + h_r, g_r = rotate_integrals(ints.h1, ints.g, U) + + H0 = to_dense(ints.h1, ints.g, ints.e_core, N) + Hr = to_dense(h_r, g_r, ints.e_core, N) + assert np.allclose(np.sort(np.linalg.eigvalsh(H0)), + np.sort(np.linalg.eigvalsh(Hr)), atol=1e-9) + + +def test_optimizer_output_preserves_sector_energy(): + """The actual optimizer output preserves the ``n``-electron spectrum. + + The joint optimizer also rotates the orbital basis, so the sector block is + written in a different basis and is not entrywise identical; the physical + invariant is the sector spectrum (and in particular the ground energy). + """ + N = 3 + ints = synthetic_integrals(N=N, seed=6, e_core=1.1, n_elec=N) + res = reduce_hamiltonian_weight(ints, optimize_orbitals=True) + + H0 = to_dense(ints.h1, ints.g, ints.e_core, N) + H1 = to_dense(res.h1, res.g, res.e_core, N) + n_diag = number_diag(N) + + blk0 = _sector_block(H0, n_diag, ints.n_elec) + blk1 = _sector_block(H1, n_diag, ints.n_elec) + spec0 = np.sort(np.linalg.eigvalsh(blk0)) + spec1 = np.sort(np.linalg.eigvalsh(blk1)) + assert abs(spec0[0] - spec1[0]) < 1e-9 # ground energy + assert np.allclose(spec0, spec1, atol=1e-9) # whole sector spectrum + + +# ---------------------------------------------------------------------------- +# Consistency and the BLISS-only path. +# ---------------------------------------------------------------------------- + +def test_fast_matches_dense_on_random_transforms(): + """The internal closed form matches :func:`lambda_pauli` on rotated+shifted.""" + N = 4 + ints = synthetic_integrals(N=N, seed=8, n_elec=N) + rng = np.random.default_rng(11) + for _ in range(5): + X = rng.standard_normal((N, N)) * 0.3 + U = sla.expm(X - X.T) + h_r, g_r = rotate_integrals(ints.h1, ints.g, U) + s = rng.standard_normal((N, N)) + s = (s + s.T) / 2.0 + nu = float(rng.standard_normal()) + c = float(rng.standard_normal()) + hp, gp = shifted_integrals(h_r, g_r, s, nu, c, ints.n_elec) + lam = lambda_pauli(hp, gp)["total"] + # re-evaluate the three parts independently as a cross-check + parts = lambda_pauli(hp, gp) + assert abs(parts["total"] + - (parts["one_body"] + parts["opposite_spin"] + + parts["same_spin"])) < 1e-12 * (abs(lam) + 1.0) + + +def test_orbital_opt_false_uses_identity_rotation(): + """With ``optimize_orbitals=False`` the applied rotation is the identity.""" + ints = synthetic_integrals(seed=0) + res = reduce_hamiltonian_weight(ints, optimize_orbitals=False) + assert np.allclose(res.U, np.eye(ints.h1.shape[0])) + assert isinstance(res, ReducedIntegrals) + + +def test_nu_nonneg_bound_respected(): + """The ``nu >= 0`` restriction is honored at the optimum.""" + ints = synthetic_integrals(seed=0) + res = reduce_hamiltonian_weight(ints, optimize_orbitals=False, + nu_nonneg=True) + assert res.nu >= -1e-9 diff --git a/prqpe/tests/test_cgs.py b/prqpe/tests/test_cgs.py new file mode 100644 index 0000000..aeda48b --- /dev/null +++ b/prqpe/tests/test_cgs.py @@ -0,0 +1,150 @@ +"""First-principles tests for the ground-state Trotter constant estimator. + +The benchmark is a small transverse-field Ising chain, whose coupling and field +terms do not commute, so the symmetric second-order Trotter step carries a +nonzero ground-state bias. The expected behavior is derived from the formula +order, not from any tabulated value: the bias scales as ``delta^2``, so the +fitted exponent must be close to two and the fitted constant positive and +finite. Everything here is ``numpy``-only. +""" + +import numpy as np +import pytest +import scipy.linalg as sla + +from qhat.prqpe.cgs import cgs_rule, estimate_cgs, fit_cgs +from qhat.prqpe.cgs.trotter_matrix import ( + effective_ground_energy, + ground_state, + trot2_matrix, +) +from qhat.prqpe.representation.pauli import PauliHamiltonian + + +def transverse_field_ising(n_qubits, J=1.0, g=0.6): + """A transverse-field Ising chain ``H = -J sum ZZ - g sum X`` (open chain). + + The ``ZZ`` couplings and the ``X`` fields do not commute, so the symmetric + Trotter step has a nonzero ground-state bias. + """ + xs, zs, cs = [], [], [] + for q in range(n_qubits - 1): + zz = (1 << (n_qubits - 1 - q)) | (1 << (n_qubits - 1 - (q + 1))) + xs.append(0) + zs.append(zz) + cs.append(-J) + for q in range(n_qubits): + xs.append(1 << (n_qubits - 1 - q)) + zs.append(0) + cs.append(-g) + return PauliHamiltonian(n_qubits, xs, zs, cs) + + +# ---------------------------------------------------------------------------- +# The bias power law: second-order Trotter gives p ~= 2. +# ---------------------------------------------------------------------------- + +def test_estimated_exponent_is_about_two(): + # On a clean step-size grid the fitted exponent of |E0 - E0~| = C_gs*delta^p + # must be close to the second-order value p = 2. + ham = transverse_field_ising(4) + e0, psi0 = ground_state(ham) + deltas = np.geomspace(0.2, 0.02, 7) + errors = np.array( + [abs(e0 - effective_ground_energy(trot2_matrix(ham, d), d, psi0)[0]) + for d in deltas] + ) + _c_gs, p, _res = fit_cgs(errors, deltas) + assert abs(p - 2.0) < 0.15 + + +def test_estimate_cgs_is_positive_and_finite(): + c_gs = estimate_cgs(transverse_field_ising(4)) + assert np.isfinite(c_gs) + assert c_gs > 0.0 + + +def test_estimate_cgs_matches_manual_fit(): + # The convenience entry point reproduces the explicit fit on its own grid. + ham = transverse_field_ising(3) + e0, psi0 = ground_state(ham) + from qhat.prqpe.cgs import DEFAULT_DELTAS + + errors = np.array( + [abs(e0 - effective_ground_energy(trot2_matrix(ham, d), d, psi0)[0]) + for d in DEFAULT_DELTAS] + ) + c_gs, _p, _res = fit_cgs(errors, DEFAULT_DELTAS) + assert estimate_cgs(ham) == pytest.approx(c_gs) + + +# ---------------------------------------------------------------------------- +# Trotter-unitary and effective-energy sanity. +# ---------------------------------------------------------------------------- + +def test_trot2_matrix_is_unitary(): + ham = transverse_field_ising(3) + U = trot2_matrix(ham, 0.1) + np.testing.assert_allclose(U.conj().T @ U, np.eye(ham.dim), atol=1e-12) + + +def test_trot2_matrix_is_time_reversal_symmetric(): + # The symmetric Strang step obeys S_2(-delta) = S_2(delta)^H. + ham = transverse_field_ising(3) + for delta in (0.2, 0.1, 0.05): + forward = trot2_matrix(ham, delta) + backward = trot2_matrix(ham, -delta) + np.testing.assert_allclose(forward, backward.conj().T, atol=1e-12) + + +def test_trot2_matrix_local_error_is_third_order(): + # ||S_2(delta) - exp(-i delta H)||_2 scales as delta^3 for the symmetric step. + ham = transverse_field_ising(3) + H = ham.to_dense() + deltas = np.geomspace(0.2, 0.025, 5) + errors = np.array( + [np.linalg.norm(trot2_matrix(ham, d) - sla.expm(-1j * d * H), 2) + for d in deltas] + ) + slope = np.polyfit(np.log(deltas), np.log(errors), 1)[0] + assert abs(slope - 3.0) < 0.2 + + +def test_effective_energy_converges_to_exact_as_step_shrinks(): + # The Trotter bias vanishes with the step, so E0~ -> E0 as delta -> 0 and + # the selected eigenstate's overlap with the exact ground state -> 1. + ham = transverse_field_ising(3) + e0, psi0 = ground_state(ham) + e_coarse, ov_coarse = effective_ground_energy(trot2_matrix(ham, 0.2), 0.2, psi0) + e_fine, ov_fine = effective_ground_energy(trot2_matrix(ham, 0.01), 0.01, psi0) + assert abs(e_fine - e0) < abs(e_coarse - e0) + assert abs(e_fine - e0) < 1e-4 + assert ov_fine > ov_coarse + assert ov_fine == pytest.approx(1.0, abs=1e-6) + + +def test_ground_state_matches_dense_diagonalization(): + ham = transverse_field_ising(4) + e0, psi0 = ground_state(ham) + evals = np.linalg.eigvalsh(ham.to_dense()) + assert e0 == pytest.approx(evals[0]) + # psi0 is a unit eigenvector at energy E0. + Hpsi = ham.to_dense() @ psi0 + assert np.vdot(psi0, psi0) == pytest.approx(1.0) + np.testing.assert_allclose(Hpsi, e0 * psi0, atol=1e-8) + + +# ---------------------------------------------------------------------------- +# Extrapolation rule. +# ---------------------------------------------------------------------------- + +def test_cgs_rule_is_the_power_law(): + assert cgs_rule(10.0, 0.5, 1.3) == pytest.approx(0.5 * 10.0 ** 1.3) + assert cgs_rule(64.0, 0.02, 1.0) == pytest.approx(0.02 * 64.0) + + +def test_cgs_rule_scales_with_lambda(): + # Doubling lambda multiplies C_gs by 2**b. + base = cgs_rule(50.0, 0.3, 1.2) + doubled = cgs_rule(100.0, 0.3, 1.2) + assert doubled / base == pytest.approx(2.0 ** 1.2) diff --git a/prqpe/tests/test_costs.py b/prqpe/tests/test_costs.py new file mode 100644 index 0000000..339f830 --- /dev/null +++ b/prqpe/tests/test_costs.py @@ -0,0 +1,175 @@ +"""Paper-reproduction anchors and closed-form property tests for the cost core. + +The reference points are the worked FeMoco numbers and the equations they come +from. The FeMoco system has weight ``lambda = 405`` and is estimated to +precision ``eps = 0.0016`` Hartree with commuting-group length ``K = 10``. +""" + +import math + +import numpy as np +import pytest + +from qhat.prqpe import costs +from qhat.prqpe.costs import ( + minimize_delta_kappa, + optimal_kappa, + partial_toffoli_bound, + phase_state_precision, + qpe_trotter_split, + rotation_toffolis, + trotter_step_size, +) + +EPS = 0.0016 +LAM_FEMOCO = 405.0 +K = 10 + + +# Closed-form totals for the two degenerate splits, used here as cross-check +# oracles. They are exactly the L_D=0 / lam_R=0 cases of minimize_delta_kappa. +def _randomized_total(lam, eps, G_rand): + return minimize_delta_kappa(0, lam, eps, 0.1, G_det=0.0, G_rand=G_rand).toffoli_total + + +def _deterministic_total(L, eps, C_gs, G_det): + return minimize_delta_kappa(L, 0.0, eps, C_gs, G_det=G_det, G_rand=2.5).toffoli_total + + +# ---------------------------------------------------------------------------- +# FeMoco paper anchors (weight-independent: derivable from lambda and eps). +# ---------------------------------------------------------------------------- + +def test_femoco_dyadic_exponent(): + # J = ceil(log2(lambda / (pi*eps))). + assert phase_state_precision(LAM_FEMOCO, EPS) == 17 + + +def test_femoco_toffolis_per_rotation(): + # 1 + (J-2)/K = 1 + 15/10 = 2.5. + assert rotation_toffolis(17, K) == pytest.approx(2.5) + + +def test_dyadic_exponent_floored_at_one_bit(): + # A randomized weight below the precision scale still needs one bit, keeping + # the per-rotation cost positive. + assert phase_state_precision(1e-9, EPS) == 1 + assert rotation_toffolis(phase_state_precision(1e-9, EPS), K) > 0.0 + + +def test_femoco_phasing_ancillas(): + # K + 2J - 2 = 10 + 34 - 2 = 42. + assert costs.phasing_ancillas(17, K) == 42 + + +def test_femoco_logical_qubits(): + # 152 system qubits + 1 Hadamard-test ancilla + 42 phasing ancillas. + assert costs.logical_qubits(152, 17, K) == 195 + + +def test_femoco_fully_randomized_total(): + # The fully randomized arm depends only on lambda; its value sets the scale + # the partially randomized estimate improves on (paper: ~2.8e12 Toffoli, + # the upper method in the FeMoco comparison figure). + G_rand = rotation_toffolis(17, K) + total = _randomized_total(LAM_FEMOCO, EPS, G_rand) + assert total == pytest.approx(2.83e12, rel=0.05) + + +def test_fully_randomized_matches_rte_rotation_bookkeeping(): + # Independent cross-check: the App. E.2 randomized bookkeeping counts + # r = 16.3 lambda^2/eps^2 Pauli rotations, so G_rand*r Toffolis. This must + # agree with the schedule-summed fully randomized total to ~10%. + G_rand = rotation_toffolis(17, K) + rotation_count = 16.3 * LAM_FEMOCO ** 2 / EPS ** 2 + bookkeeping = G_rand * rotation_count + closed_form = _randomized_total(LAM_FEMOCO, EPS, G_rand) + assert closed_form == pytest.approx(bookkeeping, rel=0.1) + + +# ---------------------------------------------------------------------------- +# Error-budget and Trotter-step consistency. +# ---------------------------------------------------------------------------- + +def test_qpe_trotter_split_adds_in_quadrature(): + eps_qpe, eps_trot = qpe_trotter_split(EPS, order=2) + assert eps_qpe ** 2 + eps_trot ** 2 == pytest.approx(EPS ** 2) + assert eps_qpe == pytest.approx(EPS * math.sqrt(2.0 / 3.0)) + + +@pytest.mark.parametrize("C_gs", [1e-3, 0.1, 1.4]) +def test_trotter_step_is_self_consistent(C_gs): + # The Trotter step must satisfy C_gs*delta^p = eps_trotter exactly, so that + # eps^2 = eps_qpe^2 + eps_trotter^2 holds at the optimum. + delta = trotter_step_size(EPS, C_gs, order=2) + _, eps_trot = qpe_trotter_split(EPS, order=2) + assert C_gs * delta ** 2 == pytest.approx(eps_trot) + + +# ---------------------------------------------------------------------------- +# kappa optimization. +# ---------------------------------------------------------------------------- + +@pytest.mark.parametrize("A,B", [(1.0, 1.0), (10.0, 1.0), (0.5, 3.0)]) +def test_optimal_kappa_solves_quadratic(A, B): + kappa = optimal_kappa(A, B) + # dG/dkappa = 0 <=> B*kappa^2 - 2*B*kappa - 2*A = 0. + assert B * kappa ** 2 - 2.0 * B * kappa - 2.0 * A == pytest.approx(0.0, abs=1e-9) + + +def test_optimal_kappa_no_randomized_part_is_infinite(): + assert math.isinf(optimal_kappa(1.0, 0.0)) + + +def test_optimal_kappa_no_deterministic_part_is_two(): + # A = 0 recovers the RTE-alone optimum r = 2 t^2. + assert optimal_kappa(0.0, 5.0) == pytest.approx(2.0) + + +# ---------------------------------------------------------------------------- +# Split endpoints: the optimizer reduces to the closed-form arms. +# ---------------------------------------------------------------------------- + +def test_fully_randomized_endpoint(): + G_rand = rotation_toffolis(17, K) + opt = minimize_delta_kappa(0, LAM_FEMOCO, EPS, 0.1, G_det=0.0, G_rand=G_rand) + assert opt.kappa == pytest.approx(2.0) + assert opt.eps_qpe == pytest.approx(EPS) + assert opt.toffoli_total == pytest.approx( + _randomized_total(LAM_FEMOCO, EPS, G_rand)) + + +def test_fully_deterministic_endpoint(): + L, G_det = 5000, 23.0 + opt = minimize_delta_kappa(L, 0.0, EPS, 0.1, G_det=G_det, G_rand=2.5) + assert math.isinf(opt.kappa) + assert opt.eps_qpe == pytest.approx(EPS * math.sqrt(2.0 / 3.0)) + assert opt.toffoli_total == pytest.approx( + _deterministic_total(L, EPS, 0.1, G_det)) + + +def test_bound_matches_optimizer_at_returned_point(): + # Re-evaluating the closed-form bound at the optimizer's (delta, kappa) + # reproduces its reported total. + opt = minimize_delta_kappa(2000, 40.0, EPS, 0.1, G_det=23.0, G_rand=2.5) + recomputed = partial_toffoli_bound( + 2000, 40.0, opt.delta, opt.eps_qpe, opt.kappa, 23.0, 2.5) + assert recomputed == pytest.approx(opt.toffoli_total) + + +# ---------------------------------------------------------------------------- +# Asymptotic scaling with target precision. +# ---------------------------------------------------------------------------- + +def test_fully_randomized_scales_as_inverse_eps_squared(): + G_rand = 2.5 + coarse = _randomized_total(LAM_FEMOCO, EPS, G_rand) + fine = _randomized_total(LAM_FEMOCO, EPS / 2.0, G_rand) + assert fine / coarse == pytest.approx(4.0) + + +def test_fully_deterministic_scales_as_eps_to_the_minus_three_halves(): + L, G_det = 5000, 23.0 + coarse = _deterministic_total(L, EPS, 0.1, G_det) + fine = _deterministic_total(L, EPS / 2.0, 0.1, G_det) + assert fine / coarse == pytest.approx(2.0 ** 1.5) diff --git a/prqpe/tests/test_estimate.py b/prqpe/tests/test_estimate.py new file mode 100644 index 0000000..8578e8e --- /dev/null +++ b/prqpe/tests/test_estimate.py @@ -0,0 +1,151 @@ +"""Integration and property tests for ``estimate_resources`` through the public API.""" + +import math + +import numpy as np +import pytest + +from qhat.prqpe import estimate_resources +from qhat.prqpe.costs import ( + minimize_delta_kappa, + phase_state_precision, + rotation_toffolis, +) + +EPS = 0.0016 +K = 10 + + +# Closed-form totals for the degenerate splits (the L_D=0 / lam_R=0 cases of +# minimize_delta_kappa), used here as cross-check oracles for the endpoints. +def _randomized_total(lam, eps, G_rand): + return minimize_delta_kappa(0, lam, eps, 0.1, G_det=0.0, G_rand=G_rand).toffoli_total + + +def _deterministic_total(L, eps, C_gs, G_det): + return minimize_delta_kappa(L, 0.0, eps, C_gs, G_det=G_det, G_rand=2.5).toffoli_total + + +def representative_weights(lam: float, n_terms: int, seed: int = 0) -> np.ndarray: + """A heavy-tailed positive weight array summing to ``lam``. + + Stands in for a molecular Pauli spectrum: a few large coefficients and a + long tail of small ones. Not any specific molecule. + """ + rng = np.random.default_rng(seed) + w = rng.pareto(1.5, size=n_terms) + 1e-3 + return np.sort(w)[::-1] * (lam / w.sum()) + + +# ---------------------------------------------------------------------------- +# Endpoint methods reduce to the closed-form arms. +# ---------------------------------------------------------------------------- + +def test_randomized_method_matches_closed_form(): + w = representative_weights(405.0, 3000) + est = estimate_resources(w, 152, C_gs=0.1, method="randomized") + G_rand = rotation_toffolis(phase_state_precision(405.0, EPS), K) + assert est.toffoli_count == pytest.approx( + _randomized_total(405.0, EPS, G_rand)) + assert est.metadata["J"] == 17 + assert est.metadata["kappa"] == pytest.approx(2.0) + assert est.logical_qubits == 195 # 152 + 1 + 42 + assert len(est.pareto) == 1 + + +def test_deterministic_method_has_no_randomized_part(): + w = representative_weights(405.0, 3000) + est = estimate_resources(w, 152, C_gs=0.1, method="deterministic") + assert est.metadata["lambda_R"] == 0.0 + assert math.isinf(est.metadata["kappa"]) + assert est.logical_qubits == 153 # system + Hadamard-test ancilla only + assert est.toffoli_count == pytest.approx( + _deterministic_total(len(w), EPS, 0.1, est.metadata["G_det"])) + + +# ---------------------------------------------------------------------------- +# The partially randomized sweep. +# ---------------------------------------------------------------------------- + +def test_partial_sweep_returns_full_pareto_curve(): + w = representative_weights(405.0, 4000) + est = estimate_resources(w, 152, C_gs=0.1, method="partially_randomized") + assert len(est.pareto) > 10 + L_D = [p.L_D for p in est.pareto] + assert L_D == sorted(L_D) + assert L_D[0] == 0 and L_D[-1] == len(w) + + +def test_partial_optimum_beats_both_endpoints(): + w = representative_weights(405.0, 4000) + est = estimate_resources(w, 152, C_gs=0.1, method="partially_randomized") + randomized = est.pareto[0].toffoli_total # L_D = 0 + deterministic = est.pareto[-1].toffoli_total # L_D = L + assert est.toffoli_count <= randomized + assert est.toffoli_count <= deterministic + + +def test_partial_endpoints_match_closed_forms(): + w = representative_weights(405.0, 4000) + est = estimate_resources(w, 152, C_gs=0.1, method="partially_randomized") + G_rand = rotation_toffolis(phase_state_precision(405.0, EPS), K) + assert est.pareto[0].toffoli_total == pytest.approx( + _randomized_total(405.0, EPS, G_rand)) + assert est.pareto[-1].lambda_R == 0.0 + + +def test_sweep_is_unimodal(): + # Cost as a function of L_D falls to a single interior minimum and rises. + w = representative_weights(405.0, 4000) + est = estimate_resources(w, 152, C_gs=0.1, method="partially_randomized") + totals = np.array([p.toffoli_total for p in est.pareto]) + signs = np.sign(np.diff(totals)) + signs = signs[signs != 0] + sign_changes = int(np.sum(np.diff(signs) != 0)) + assert sign_changes <= 1 # down-then-up: one change at the minimum + + +def test_partial_improves_on_fully_randomized_at_scale(): + # A FeMoco-scale problem: the optimum strictly improves on full randomization. + w = representative_weights(405.0, 4000) + est = estimate_resources(w, 152, C_gs=0.1, method="partially_randomized") + G_rand = rotation_toffolis(17, K) + assert est.toffoli_count < _randomized_total(405.0, EPS, G_rand) + assert est.metadata["L_D"] not in (0, len(w)) # genuinely partial + + +# ---------------------------------------------------------------------------- +# Reporting and validation. +# ---------------------------------------------------------------------------- + +def test_metadata_is_complete(): + w = representative_weights(405.0, 2000) + est = estimate_resources(w, 152, C_gs=0.1) + for key in ("L_D", "lambda_R", "delta", "M", "kappa", "xi", + "eps_qpe", "J", "K", "G_det", "G_rand"): + assert key in est.metadata + + +def test_qdrift_is_cheaper_than_rte(): + w = representative_weights(405.0, 3000) + rte = estimate_resources(w, 152, C_gs=0.1, method="randomized", randomizer="rte") + qdrift = estimate_resources(w, 152, C_gs=0.1, method="randomized", + randomizer="qdrift") + assert qdrift.toffoli_count == pytest.approx(0.5 * rte.toffoli_count) + + +def test_unknown_representation_rejected(): + w = representative_weights(405.0, 100) + with pytest.raises(ValueError): + estimate_resources(w, 152, C_gs=0.1, representation="double_factorized") + + +def test_unknown_method_rejected(): + w = representative_weights(405.0, 100) + with pytest.raises(ValueError): + estimate_resources(w, 152, C_gs=0.1, method="nonsense") + + +def test_empty_weights_rejected(): + with pytest.raises(ValueError): + estimate_resources([], 152, C_gs=0.1) diff --git a/prqpe/tests/test_representation.py b/prqpe/tests/test_representation.py new file mode 100644 index 0000000..f9e0dca --- /dev/null +++ b/prqpe/tests/test_representation.py @@ -0,0 +1,265 @@ +"""Correctness tests for the Pauli representation layer. + +These verify the operator container against hand-computed dense matrices in the +stated ``(x, z)``-bitmask convention, the weight bookkeeping (``lambda``, sorting), +and the cumulative-budget truncation -- including its off-by-one +boundary. The openfermion conversion path is exercised only when openfermion is +installed; everything else is ``numpy``-only. +""" + +import numpy as np +import pytest + +from qhat.prqpe.representation import ( + build_representation, + sorted_weights, + truncate_small_terms, +) +from qhat.prqpe.representation.pauli import PauliHamiltonian, popcount + +# Single-qubit Pauli matrices, qubit 0 = the only (most significant) bit. +I2 = np.eye(2, dtype=complex) +X = np.array([[0, 1], [1, 0]], dtype=complex) +Y = np.array([[0, -1j], [1j, 0]], dtype=complex) +Z = np.array([[1, 0], [0, -1]], dtype=complex) + + +# ---------------------------------------------------------------------------- +# Single-qubit dense matrices against the (x, z) convention. +# ---------------------------------------------------------------------------- + +def test_single_qubit_x(): + # x = 1, z = 0 -> P = i^0 X^1 Z^0 = X. + ham = PauliHamiltonian(1, [1], [0], [1.0]) + np.testing.assert_allclose(ham.to_dense(), X) + + +def test_single_qubit_z(): + # x = 0, z = 1 -> P = i^0 X^0 Z^1 = Z. + ham = PauliHamiltonian(1, [0], [1], [1.0]) + np.testing.assert_allclose(ham.to_dense(), Z) + + +def test_single_qubit_y(): + # x = z = 1 -> P = i^popcount(1) X Z = i X Z = Y. + ham = PauliHamiltonian(1, [1], [1], [1.0]) + np.testing.assert_allclose(ham.to_dense(), Y) + + +def test_single_qubit_identity(): + # x = z = 0 -> P = I. + ham = PauliHamiltonian(1, [0], [0], [1.0]) + np.testing.assert_allclose(ham.to_dense(), I2) + + +def test_coefficient_scales_matrix(): + ham = PauliHamiltonian(1, [0], [1], [-2.5]) + np.testing.assert_allclose(ham.to_dense(), -2.5 * Z) + + +# ---------------------------------------------------------------------------- +# Two-qubit term, qubit 0 = MSB (left tensor factor). +# ---------------------------------------------------------------------------- + +def test_two_qubit_xz(): + # X on qubit 0 (bit 1) tensor Z on qubit 1 (bit 0): x = 0b10 = 2, z = 0b01 = 1. + # popcount(2 & 1) = 0, so the head phase is i^0 = 1 and P = X (x) Z (x) Z. + ham = PauliHamiltonian(2, [2], [1], [1.0]) + np.testing.assert_allclose(ham.to_dense(), np.kron(X, Z)) + + +def test_two_qubit_yy(): + # Y on both qubits: x = 0b11 = 3, z = 0b11 = 3, head phase i^popcount(3) = i^2 = -1. + ham = PauliHamiltonian(2, [3], [3], [1.0]) + np.testing.assert_allclose(ham.to_dense(), np.kron(Y, Y)) + + +def test_two_qubit_sum_is_hermitian(): + ham = PauliHamiltonian(2, [2, 3], [1, 3], [0.5, -0.75]) + dense = ham.to_dense() + np.testing.assert_allclose(dense, 0.5 * np.kron(X, Z) - 0.75 * np.kron(Y, Y)) + np.testing.assert_allclose(dense, dense.conj().T) + + +def test_to_sparse_matches_to_dense(): + ham = PauliHamiltonian(2, [2, 3, 0], [1, 3, 2], [0.5, -0.75, 1.25]) + np.testing.assert_allclose(ham.to_sparse().toarray(), ham.to_dense()) + + +# ---------------------------------------------------------------------------- +# popcount and the phase head. +# ---------------------------------------------------------------------------- + +def test_popcount_matches_python(): + vals = np.array([0, 1, 2, 3, 7, 255, 1023], dtype=np.int64) + expected = np.array([bin(int(v)).count("1") for v in vals]) + np.testing.assert_array_equal(popcount(vals), expected) + + +# ---------------------------------------------------------------------------- +# pauli_apply_matrix and rotation_apply_matrix against explicit matrices. +# ---------------------------------------------------------------------------- + +def test_pauli_apply_matrix_left_multiplies(): + rng = np.random.default_rng(1) + ham = PauliHamiltonian(2, [2, 3], [1, 3], [1.0, 1.0]) + M = rng.standard_normal((4, 4)) + 1j * rng.standard_normal((4, 4)) + for l in range(ham.L): + P = ham.term_matrix(l).toarray() + np.testing.assert_allclose(ham.pauli_apply_matrix(l, M), P @ M) + + +def test_rotation_apply_matrix_matches_expm(): + from scipy.linalg import expm + + rng = np.random.default_rng(2) + ham = PauliHamiltonian(2, [2, 3], [1, 3], [1.0, -1.0]) + M = rng.standard_normal((4, 4)) + 1j * rng.standard_normal((4, 4)) + phi = 0.37 + for l in range(ham.L): + P = ham.term_matrix(l).toarray() + expected = expm(-1j * phi * P) @ M + np.testing.assert_allclose(ham.rotation_apply_matrix(l, phi, M), expected, + atol=1e-12) + + +# ---------------------------------------------------------------------------- +# Weight bookkeeping: lambda, sorting. +# ---------------------------------------------------------------------------- + +def test_lambda_is_sum_abs_coeffs(): + coeffs = [0.5, -2.0, 1.5, -0.25] + ham = PauliHamiltonian(2, [0, 1, 2, 3], [0, 1, 2, 3], coeffs) + assert ham.lambda_weight == pytest.approx(sum(abs(c) for c in coeffs)) + + +def test_sorted_by_weight_is_descending(): + coeffs = np.array([0.5, -2.0, 1.5, -0.25, 3.0]) + ham = PauliHamiltonian(3, [0, 1, 2, 3, 4], [0, 1, 2, 3, 4], coeffs) + ordered = ham.sorted_by_weight() + w = np.abs(ordered.coeffs) + assert np.all(np.diff(w) <= 0) + np.testing.assert_array_equal(w, np.sort(np.abs(coeffs))[::-1]) + + +def test_sorted_weights_helper_contract(): + coeffs = np.array([0.5, -2.0, 1.5, -0.25, 3.0]) + ham = PauliHamiltonian(3, [0, 1, 2, 3, 4], [0, 1, 2, 3, 4], coeffs) + np.testing.assert_array_equal( + sorted_weights(ham), np.abs(ham.sorted_by_weight().coeffs)) + assert np.all(np.diff(sorted_weights(ham)) <= 0) + + +def test_subset_preserves_order(): + ham = PauliHamiltonian(2, [0, 1, 2, 3], [3, 2, 1, 0], [1.0, 2.0, 3.0, 4.0]) + sub = ham.subset([3, 0]) + np.testing.assert_array_equal(sub.coeffs, [4.0, 1.0]) + np.testing.assert_array_equal(sub.x_masks, [3, 0]) + np.testing.assert_array_equal(sub.z_masks, [0, 3]) + + +# ---------------------------------------------------------------------------- +# truncate_small_terms: cumulative-budget rule and its boundary. +# ---------------------------------------------------------------------------- + +def test_truncate_drops_smallest_within_budget(): + # Tail weights 1e-4 + 2e-4 + 3e-4 = 6e-4 <= 1e-3; adding the next (5e-4) + # would exceed the budget, so exactly three terms are dropped. + coeffs = np.array([1.0, 0.5, 5e-4, 3e-4, 2e-4, 1e-4]) + ham = PauliHamiltonian(3, np.arange(6), np.arange(6), coeffs) + out = truncate_small_terms(ham, weight_threshold=1e-3) + np.testing.assert_array_equal(np.abs(out.coeffs), [1.0, 0.5, 5e-4]) + assert out.L == 3 + + +def test_truncate_boundary_exactly_at_budget_drops_term(): + # Cumulative dropped weight hits the budget exactly: side="right" keeps it + # within budget, so the term at the boundary is dropped. + coeffs = np.array([1.0, 6e-4, 4e-4]) # tail 4e-4, then 4e-4 + 6e-4 = 1e-3 + ham = PauliHamiltonian(2, [0, 1, 2], [0, 1, 2], coeffs) + out = truncate_small_terms(ham, weight_threshold=1e-3) + # 4e-4 dropped (cum 4e-4 <= 1e-3); 6e-4 dropped (cum 1e-3 <= 1e-3); 1.0 kept. + np.testing.assert_array_equal(np.abs(out.coeffs), [1.0]) + assert out.L == 1 + + +def test_truncate_off_by_one_just_over_budget_keeps_term(): + # The same shape, but the second-smallest is nudged so the cumulative drop + # would be 1.0001e-3 > 1e-3: that term must be kept (off-by-one boundary). + coeffs = np.array([1.0, 6.01e-4, 4e-4]) # 4e-4 + 6.01e-4 = 1.001e-3 > 1e-3 + ham = PauliHamiltonian(2, [0, 1, 2], [0, 1, 2], coeffs) + out = truncate_small_terms(ham, weight_threshold=1e-3) + np.testing.assert_array_equal(np.abs(out.coeffs), [1.0, 6.01e-4]) + assert out.L == 2 + + +def test_truncate_drops_nothing_when_budget_below_smallest(): + coeffs = np.array([1.0, 0.5, 0.01]) + ham = PauliHamiltonian(2, [0, 1, 2], [0, 1, 2], coeffs) + out = truncate_small_terms(ham, weight_threshold=1e-3) + assert out.L == ham.L + np.testing.assert_array_equal(np.abs(out.coeffs), [1.0, 0.5, 0.01]) + + +def test_truncate_total_dropped_weight_within_budget(): + rng = np.random.default_rng(7) + coeffs = rng.standard_normal(50) * 1e-3 + ham = PauliHamiltonian(6, np.arange(50), np.arange(50), coeffs) + budget = 5e-3 + out = truncate_small_terms(ham, weight_threshold=budget) + dropped = ham.lambda_weight - out.lambda_weight + assert dropped <= budget + 1e-15 + assert np.all(np.diff(np.abs(out.coeffs)) <= 0) + + +def test_truncate_output_is_sorted_descending(): + coeffs = np.array([0.5, 1.0, 2e-4, 0.3]) + ham = PauliHamiltonian(2, [0, 1, 2, 3], [0, 1, 2, 3], coeffs) + out = truncate_small_terms(ham, weight_threshold=1e-3) + assert np.all(np.diff(np.abs(out.coeffs)) <= 0) + + +# ---------------------------------------------------------------------------- +# build_representation: pass-through and the unimplemented arm. +# ---------------------------------------------------------------------------- + +def test_build_representation_passthrough(): + ham = PauliHamiltonian(2, [2, 3], [1, 3], [0.5, -0.75]) + assert build_representation(ham) is ham + assert build_representation(ham, representation="pauli") is ham + + +def test_build_representation_unknown_raises(): + ham = PauliHamiltonian(1, [1], [0], [1.0]) + with pytest.raises(ValueError): + build_representation(ham, representation="qubitization") + + +# ---------------------------------------------------------------------------- +# openfermion conversion path (skipped if openfermion is absent). +# ---------------------------------------------------------------------------- + +def test_build_representation_from_qubit_operator(): + openfermion = pytest.importorskip("openfermion") + from openfermion import QubitOperator, get_sparse_operator + + qop = (0.5 * QubitOperator("X0 Z1") + - 0.75 * QubitOperator("Y0 Y1") + + 1.3 * QubitOperator("")) # identity term -> dropped + ham = build_representation(qop) + assert ham.n_qubits == 2 + assert ham.L == 2 # identity dropped + + # Matches openfermion's own sparse operator on the non-identity part. + nonconst = 0.5 * QubitOperator("X0 Z1") - 0.75 * QubitOperator("Y0 Y1") + ref = get_sparse_operator(nonconst, n_qubits=2).toarray() + np.testing.assert_allclose(ham.to_dense(), ref, atol=1e-12) + + +def test_qubit_operator_lambda_matches_coeffs(): + pytest.importorskip("openfermion") + from openfermion import QubitOperator + + qop = 0.3 * QubitOperator("X0") - 1.2 * QubitOperator("Z0 Z1") + ham = build_representation(qop) + assert ham.lambda_weight == pytest.approx(0.3 + 1.2) diff --git a/prqpe/tests/test_schedule.py b/prqpe/tests/test_schedule.py new file mode 100644 index 0000000..18e8d3b --- /dev/null +++ b/prqpe/tests/test_schedule.py @@ -0,0 +1,117 @@ +"""Property tests for the robust-phase-estimation schedule.""" + +import math + +import pytest + +from qhat.prqpe import schedule +from qhat.prqpe.schedule import ( + RPESchedule, + num_rounds, + num_rounds_randomized, + resolve_depth_tradeoff, + round_circuit_toffolis, + rpe_schedule, +) + + +# ---------------------------------------------------------------------------- +# Depth-reduction factor xi. +# ---------------------------------------------------------------------------- + +def test_xi_none_means_no_reduction(): + assert resolve_depth_tradeoff(overlap=1.0, xi=None) == 1.0 + + +def test_xi_passes_through_for_exact_eigenstate(): + assert resolve_depth_tradeoff(overlap=1.0, xi=0.1) == 0.1 + + +def test_xi_out_of_range_rejected(): + with pytest.raises(ValueError): + resolve_depth_tradeoff(overlap=1.0, xi=1.5) + + +def test_xi_below_overlap_floor_rejected(): + # For overlap < 1, xi must exceed arcsin((1-overlap)/overlap). + with pytest.raises(ValueError): + resolve_depth_tradeoff(overlap=0.9, xi=1e-4) + + +def test_low_overlap_rejected(): + with pytest.raises(ValueError): + resolve_depth_tradeoff(overlap=0.4, xi=0.5) + + +# ---------------------------------------------------------------------------- +# Round count. +# ---------------------------------------------------------------------------- + +def test_num_rounds_is_ceil_log2(): + eps_qpe, delta = 1.3e-3, 0.02 + expected = math.ceil(math.log2(0.08 * math.pi / (eps_qpe * delta))) + assert num_rounds(eps_qpe, delta) == expected + + +def test_xi_reduces_round_count(): + eps_qpe, delta = 1.3e-3, 0.02 + assert num_rounds(eps_qpe, delta, xi=0.1) <= num_rounds(eps_qpe, delta, xi=1.0) + + +def test_randomized_round_count_matches_paper_femoco(): + # FeMoco fully randomized at xi = 0.1: M = ceil(log2(xi*lambda/eps)) = 15. + assert num_rounds_randomized(405.0, 0.0016, xi=0.1) == 15 + + +# ---------------------------------------------------------------------------- +# Sample schedule. +# ---------------------------------------------------------------------------- + +def test_samples_decrease_toward_deepest_round(): + sched = RPESchedule(M=10, inflation=1.0) + counts = [sched.samples(m) for m in range(sched.M + 1)] + assert all(a >= b for a, b in zip(counts, counts[1:])) + assert sched.samples(sched.M) == 11 # N_M at inflation 1 + + +def test_num_circuits_is_twice_the_sample_sum(): + sched = RPESchedule(M=8, inflation=math.e) + assert sched.num_circuits == 2 * sum(sched.samples(m) for m in range(9)) + + +def test_rte_inflation_is_e(): + # RTE optimal kappa = 2 inflates samples by e^(2/2) = e. + sched = rpe_schedule(M=5, kappa=2.0, overlap=1.0) + assert sched.inflation == pytest.approx(math.e) + + +def test_deterministic_has_no_inflation(): + sched = rpe_schedule(M=5, kappa=math.inf, overlap=1.0) + assert sched.inflation == pytest.approx(1.0) + + +def test_lower_overlap_inflates_samples(): + high = rpe_schedule(M=5, kappa=2.0, overlap=1.0) + low = rpe_schedule(M=5, kappa=2.0, overlap=0.8) + assert low.inflation > high.inflation + + +# ---------------------------------------------------------------------------- +# Per-circuit cost. +# ---------------------------------------------------------------------------- + +def test_deterministic_depth_uses_halving(): + # The deterministic part runs 2^(m-1) halved Trotter steps. + circuit = dict(G_det=23.0, G_rand=2.5, L_D=100, lam_R=0.0, + delta=0.02, kappa=math.inf) + assert round_circuit_toffolis(3, **circuit) == pytest.approx( + 23.0 * 2 * 100 * 2.0 ** 2) + + +def test_deepest_circuit_is_round_M(): + sched = RPESchedule(M=6, inflation=math.e) + circuit = dict(G_det=23.0, G_rand=2.5, L_D=500, lam_R=30.0, + delta=0.02, kappa=4.0) + deepest = schedule.max_toffolis_per_circuit(sched, **circuit) + assert deepest == round_circuit_toffolis(6, **circuit) + assert deepest >= round_circuit_toffolis(5, **circuit) diff --git a/pyproject.toml b/pyproject.toml index 99771a1..d6cf4c3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,5 +28,5 @@ dependencies = [ ] [tool.setuptools] -packages = ["qhat", "qhat.analysis", "qhat.common", "qhat.hamiltonian_generator"] -package-dir = {"qhat" = ".", "qhat.analysis" = "analysis", "qhat.common" = "common", "qhat.hamiltonian_generator" = "hamiltonian_generator"} +packages = ["qhat", "qhat.analysis", "qhat.common", "qhat.hamiltonian_generator", "qhat.prqpe", "qhat.prqpe.representation", "qhat.prqpe.cgs", "qhat.prqpe.bliss"] +package-dir = {"qhat" = ".", "qhat.analysis" = "analysis", "qhat.common" = "common", "qhat.hamiltonian_generator" = "hamiltonian_generator", "qhat.prqpe" = "prqpe", "qhat.prqpe.representation" = "prqpe/representation", "qhat.prqpe.cgs" = "prqpe/cgs", "qhat.prqpe.bliss" = "prqpe/bliss"}