diff --git a/analysis/README.md b/analysis/README.md index 0b1f6aaf2..0fa02cb85 100644 --- a/analysis/README.md +++ b/analysis/README.md @@ -320,10 +320,49 @@ qubits alone does not determine the window, so set either `algorithm.energy_erro e.g. `num_precision_qubits = 3` with `num_confidence_qubits = 3`, or `num_phase_qubits = 6` with `num_confidence_qubits = 3`. Setting all three is allowed if they are consistent; any other combination is an error. "rectangular" accepts `algorithm.num_phase_qubits` alone (and uses -only the total if the split is also given). When sizing from the error, the number of -confidence qubits is chosen from the paper's Theorem 1 (a bound on the failure probability -averaged over phases); this rule is expected to change. `num_precision_qubits` and -`num_confidence_qubits` apply only to "QPE: QHAT flexible". +only the total if the split is also given). `num_precision_qubits` and `num_confidence_qubits` +apply only to "QPE: QHAT flexible". + +#### Sizing tapered QPE for several runs + +The initial state usually overlaps the target eigenstate (e.g. the ground state) only partly, so a +single run of QPE often returns some other eigenvalue. For the tapered `algorithm.ancilla_prep` +choices, QHAT therefore plans N runs of the circuit, of which the K-th smallest energy estimate is +reported. This uses the multi-run analysis by Yigit Subasi and Alan Bu (unpublished), +implemented in `common/tapered_qpe_sizing.py`. Two options control it: + +- `algorithm.overlap`: the assumed overlap (probability) of the initial state with the target + eigenstate, with 0 < overlap <= 1; default 0.5. The overlap is with an eigenstate of the + implemented unitary (e.g. the Trotterized one); its error from the exact time evolution belongs + in the separate unitary error budget. +- `algorithm.force_lowest_estimate`: if true, report the smallest estimate (K = 1); default false, + which lets the sizing choose K. + +When sizing from `algorithm.energy_error` and `algorithm.probability_of_failure`, the precision +qubits are the fewest that meet the energy error, and the confidence qubits, N and K are those of +lowest cost N x 2^(number of phase qubits) for which the reported estimate fails with probability +at most `algorithm.probability_of_failure`. The log reports the probability of failure for one +run and for the N runs, and a table of the N-run failure if the true overlap were 0.01, 0.1, 0.2, +..., 1.0. With an overlap of 1, a single run suffices and the sizing reduces to one run's failure +bound. For example, energy error 0.2 (a phase error of 0.118 turns) and probability of failure +1e-3 with the default overlap 0.5 give 3 precision and 4 confidence qubits, run 14 times +("dpss"), reporting the 2nd-smallest estimate; a single run fails with probability 0.5. + +The circuit and all analyses (resource estimates, simulations) are for a **single** run; nothing +is multiplied by N. If no plan costs at most 1e18, QHAT logs an error and falls back to the +single-run sizing (as if the overlap were 1). If qubit counts are given, the circuit uses them, +and if `algorithm.probability_of_failure` is set the log reports the fewest runs that meet it; +the energy error bound then comes from the number of precision qubits. + +The probabilities of failure use the paper's Theorem 1, a bound on one run's failure averaged over +where the true phase lies relative to the QPE grid (times 4 for "truncated dpss" and +"bandwidth-optimized"). The paper's random phase shift (Sec. VI) would make this a worst-case +bound, but QHAT does not apply it. + +The analysis also assumes the energy window leaves at least the energy error of margin at both +ends: above the highest eigenvalue, so that a correct estimate of it cannot wrap around to a low +energy, and below the target, so that its own estimates cannot wrap around to a high one. The +driver does not yet guarantee this margin. ### Analyzing an Algorithm diff --git a/analysis/algorithm.py b/analysis/algorithm.py index 3d4df645f..61bc824b0 100644 --- a/analysis/algorithm.py +++ b/analysis/algorithm.py @@ -40,9 +40,16 @@ from qhat.common.flexible_qpe import FlexibleQPE from qhat.common.hermiticity import require_real_scalar from qhat.common.qpe_window_state import precision_bits, RectangularWindowState +from qhat.common.tapered_qpe_sizing import ( + COST_CEILING, + fewest_runs, + multi_run_failure, + size_tapered_qpe, +) from qhat.common.tapered_window_state import ( BandwidthOptimizedWindowState, DPSSWindowState, + TaperedWindowStateBase, TruncatedDPSSWindowState, ) @@ -183,6 +190,156 @@ def build_qpe_qualtran_textbook( } INVERSE_QFTS = {"textbook": None} +# Multi-run sizing of tapered windows (see qhat.common.tapered_qpe_sizing). +DEFAULT_OVERLAP = 0.5 +OVERLAP_TABLE = (0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0) +MULTI_RUN_OPTIONS = ("overlap", "force_lowest_estimate") +SINGLE_RUN_NOTE = "-- NOTE: the analyses below report on a single run of the algorithm." + + +def warn_unused_multi_run_options(config_algorithm: AlgorithmConfiguration, where: str): + for name in MULTI_RUN_OPTIONS: + if getattr(config_algorithm, name) is not None: + logger.warning( + f"algorithm.{name} is ignored with {where}: it applies only to \"QPE: QHAT " + "flexible\" with a tapered ancilla_prep (\"dpss\", \"truncated dpss\" or " + "\"bandwidth-optimized\") and algorithm.probability_of_failure.") + + +def multi_run_options(config_algorithm: AlgorithmConfiguration): + """Check algorithm.overlap and algorithm.force_lowest_estimate and fill in their defaults. + + Returns `(overlap, overlap_is_default, force_lowest_estimate, force_is_default)`, e.g. + `(0.5, True, False, True)` when neither is set. + """ + overlap = config_algorithm.overlap + overlap_is_default = overlap is None + if overlap_is_default: + overlap = DEFAULT_OVERLAP + else: + overlap = require_real_scalar("algorithm.overlap", overlap) + if not 0 < overlap <= 1: + raise ValueError(f"algorithm.overlap must be in (0, 1], got {overlap}.") + force = config_algorithm.force_lowest_estimate + force_is_default = force is None + if force_is_default: + force = False + elif not isinstance(force, bool): + raise TypeError( + f"algorithm.force_lowest_estimate must be True or False, got {force!r}.") + return overlap, overlap_is_default, force, force_is_default + + +def log_multi_run_options(overlap, overlap_is_default, force, force_is_default): + def source(is_default): + return "default" if is_default else "user-specified" + logger.info(f"-- overlap with target eigenstate = {overlap:g} ({source(overlap_is_default)})") + if force: + logger.info("-- report the smallest of N runs (K = 1): forced by " + f"algorithm.force_lowest_estimate ({source(force_is_default)})") + else: + logger.info("-- report the K-th smallest of N runs: K chosen by sizing " + f"({source(force_is_default)})") + + +def smallest(k: int) -> str: + """'smallest', '2nd-smallest', '3rd-smallest', '11th-smallest', ...""" + if k == 1: + return "smallest" + suffix = "th" if k % 100 in (11, 12, 13) else {1: "st", 2: "nd", 3: "rd"}.get(k % 10, "th") + return f"{k}{suffix}-smallest" + + +def log_multi_run_plan(window_cls, plan, error_bound: str): + """Log the failure probability for 1 run and for the plan's N runs, then for N runs at the + overlaps in OVERLAP_TABLE. `error_bound` is e.g. "energy error <= 0.2".""" + runs, which = plan.num_runs, smallest(plan.quantile) + times = "1 time" if runs == 1 else f"{runs} times" + logger.info("-- if you run the circuit 1 time:") + logger.info(f"-- probability of failure <= {plan.single_run_failure:.2e}, {error_bound}") + logger.info(f"-- if you run the circuit {times} and take the {which} estimate:") + logger.info(f"-- probability of failure <= {plan.multi_run_failure:.2e}, {error_bound}") + logger.info(f"-- if the true overlap differs from {plan.overlap:g} (same {runs} runs, {which} " + f"estimate; {error_bound} in every row):") + logger.info("-- overlap probability of failure") + for overlap in OVERLAP_TABLE: + failure = multi_run_failure(window_cls, plan.num_precision_qubits, + plan.num_confidence_qubits, runs, plan.quantile, overlap) + logger.info(f"-- {overlap:<9.2f} <= {failure:.2e}") + + +def counts_text(num_precision, num_confidence): + return (f"{num_precision} precision + {num_confidence} confidence qubits = " + f"{num_precision + num_confidence} phase qubits") + + +def size_tapered_window(config_algorithm, window_cls, phase_error, probability_of_failure): + """`(num_precision_qubits, num_confidence_qubits)` from the multi-run sizing, logging the + plan. If no plan exists, logs an error and falls back on `window_cls.from_requirements` (a + single run, as if the overlap were 1).""" + overlap, overlap_is_default, force, force_is_default = multi_run_options(config_algorithm) + log_multi_run_options(overlap, overlap_is_default, force, force_is_default) + plan = size_tapered_qpe(window_cls, phase_error, probability_of_failure, overlap, force) + if plan is None: + fallback = window_cls.from_requirements(phase_error, probability_of_failure) + logger.error( + f"Tapered QPE sizing found no plan meeting probability of failure " + f"{probability_of_failure} at overlap {overlap:g} with cost N x 2^(phase qubits) <= " + f"{COST_CEILING:g}. FALLBACK: sizing for a single run with overlap 1; the number of " + "runs and the overlap table are skipped.") + logger.info("-- FALLBACK tapered QPE sizing: " + counts_text( + fallback.num_precision_qubits, fallback.num_confidence_qubits)) + logger.info(SINGLE_RUN_NOTE) + return fallback.num_precision_qubits, fallback.num_confidence_qubits + logger.info("-- tapered QPE sizing: " + + counts_text(plan.num_precision_qubits, plan.num_confidence_qubits)) + log_multi_run_plan(window_cls, plan, f"energy error <= {config_algorithm.energy_error:.3g}") + logger.info(SINGLE_RUN_NOTE) + return plan.num_precision_qubits, plan.num_confidence_qubits + + +def report_given_tapered_counts(config_algorithm, unitary, window_cls, window): + """Log the multi-run plan for user-given tapered qubit counts: the fewest runs meeting + algorithm.probability_of_failure at the overlap. The error bound is the window's + 2^-(num_precision_qubits + 1) turns, as an energy if the unitary can convert it.""" + num_precision, num_confidence = window.num_precision_qubits, window.num_confidence_qubits + overlap, overlap_is_default, force, force_is_default = multi_run_options(config_algorithm) + phase_bound = 2.0**-(num_precision + 1) + error_bound = f"phase error <= {phase_bound:.3g} turns" + if hasattr(unitary, "phase_error_from_energy_error"): + energy_bound = phase_bound / unitary.phase_error_from_energy_error(1.0) + error_bound = f"energy error <= {energy_bound:.3g}" + energy_error = config_algorithm.energy_error + if energy_error is not None and energy_bound > energy_error * (1 + 1e-12): + logger.warning( + f"num_precision_qubits = {num_precision} bounds the energy error by " + f"{energy_bound:.3g}, more than algorithm.energy_error = {energy_error}.") + + logger.info("-- tapered QPE: user-specified " + counts_text(num_precision, num_confidence)) + probability_of_failure = config_algorithm.probability_of_failure + if probability_of_failure is None: + warn_unused_multi_run_options(config_algorithm, "no algorithm.probability_of_failure") + logger.info("-- algorithm.probability_of_failure is not set: skipping the number of runs " + "and the overlap table.") + logger.info(SINGLE_RUN_NOTE) + return + if not 0 < probability_of_failure < 1: + raise ValueError( + f"algorithm.probability_of_failure must be in (0, 1), got {probability_of_failure}.") + log_multi_run_options(overlap, overlap_is_default, force, force_is_default) + plan = fewest_runs(window_cls, num_precision, num_confidence, probability_of_failure, overlap, + force) + if plan is None: + logger.error( + f"No number of runs of the tapered QPE with {num_precision} precision and " + f"{num_confidence} confidence qubits meets probability of failure " + f"{probability_of_failure} at overlap {overlap:g} with cost N x 2^(phase qubits) <= " + f"{COST_CEILING:g}; the number of runs and the overlap table are skipped.") + else: + log_multi_run_plan(window_cls, plan, error_bound) + logger.info(SINGLE_RUN_NOTE) + + def build_qpe_qhat_flexible( config_algorithm: AlgorithmConfiguration, unitary): @@ -196,8 +353,12 @@ def build_qpe_qhat_flexible( if qft_inv_name not in INVERSE_QFTS: raise ValueError(f"Invalid QPE qft_inv \"{config_algorithm.qft_inv}\".") - components = dict(ancilla_prep=ANCILLA_PREPS[ancilla_prep_name], - qft_inv=INVERSE_QFTS[qft_inv_name]) + window_cls = ANCILLA_PREPS[ancilla_prep_name] + tapered = issubclass(window_cls, TaperedWindowStateBase) + if not tapered: + warn_unused_multi_run_options( + config_algorithm, f"algorithm.ancilla_prep = \"{ancilla_prep_name}\"") + components = dict(ancilla_prep=window_cls, qft_inv=INVERSE_QFTS[qft_inv_name]) counts = dict(num_phase_qubits=config_algorithm.num_phase_qubits, num_precision_qubits=config_algorithm.num_precision_qubits, num_confidence_qubits=config_algorithm.num_confidence_qubits) @@ -212,6 +373,15 @@ def build_qpe_qhat_flexible( f"algorithm.ancilla_prep = \"{ancilla_prep_name}\" with {given}: {error} Or set " "algorithm.energy_error and algorithm.probability_of_failure, and no qubit " "counts.") from error + if tapered: + report_given_tapered_counts(config_algorithm, unitary, window_cls, qpe.ancilla_prep) + elif tapered: + phase_error, probability_of_failure = qpe_requirements(config_algorithm, unitary) + num_precision, num_confidence = size_tapered_window( + config_algorithm, window_cls, phase_error, probability_of_failure) + qpe = FlexibleQPE.from_qubit_counts( + unitary, num_precision_qubits=num_precision, num_confidence_qubits=num_confidence, + **components) else: phase_error, probability_of_failure = qpe_requirements(config_algorithm, unitary) qpe = FlexibleQPE.from_requirements( @@ -285,6 +455,9 @@ def build_algorithm( logger.info("Beginning to construct quantum algorithm.") + if config_algorithm.method.lower() != "qpe: qhat flexible": + warn_unused_multi_run_options( + config_algorithm, f"algorithm.method = \"{config_algorithm.method}\"") if config_algorithm.method.lower() in ("qpe: qualtran textbook",): return build_qpe_qualtran_textbook(config_algorithm, unitary) elif config_algorithm.method.lower() in ("qpe: qhat flexible",): diff --git a/analysis/config_types.py b/analysis/config_types.py index 96fde62f8..5715ff425 100644 --- a/analysis/config_types.py +++ b/analysis/config_types.py @@ -200,6 +200,11 @@ def __init__(self): # Components for "QPE: QHAT flexible"; `None` yields default options self.ancilla_prep = None self.qft_inv = None + # Multi-run sizing for tapered ancilla_prep choices ("QPE: QHAT flexible" only): assumed + # overlap of the initial state with the target eigenstate (default 0.5), and whether to + # report the smallest of the runs' estimates (default False: the sizing picks which) + self.overlap = None + self.force_lowest_estimate = None def _generate_TOML_table(self): table = tomlkit.table() table["method"] = self.method @@ -210,6 +215,8 @@ def _generate_TOML_table(self): self.save_if_present(table, "energy_error") self.save_if_present(table, "ancilla_prep") self.save_if_present(table, "qft_inv") + self.save_if_present(table, "overlap") + self.save_if_present(table, "force_lowest_estimate") return table # ------------------------------------------------------------------------------------------------- diff --git a/analysis/examples/all_analyses.config b/analysis/examples/all_analyses.config index a8565caa7..d29c4e134 100644 --- a/analysis/examples/all_analyses.config +++ b/analysis/examples/all_analyses.config @@ -136,6 +136,12 @@ algorithm.method = "time evolution" # algorithm.ancilla_prep = "truncated dpss" # algorithm.num_precision_qubits = 3 # success = estimate within 2^-(3+1) turns # algorithm.num_confidence_qubits = 3 # sets the probability of failure (~1e-8 here) +# +# Tapered windows are sized for N runs of the circuit, reporting the K-th smallest estimate, so +# that the reported estimate meets probability_of_failure. The log gives N, K and the failure for +# other overlaps; all analyses are for a single run. +# algorithm.overlap = 0.5 # assumed overlap with the target eigenstate +# algorithm.force_lowest_estimate = False # True: report the smallest estimate (K = 1) # ================================================================================================= # ANALYSIS CONFIGURATION diff --git a/analysis/tests/test_flexible_qpe_algorithm.py b/analysis/tests/test_flexible_qpe_algorithm.py index 1431fda4b..431271770 100644 --- a/analysis/tests/test_flexible_qpe_algorithm.py +++ b/analysis/tests/test_flexible_qpe_algorithm.py @@ -5,10 +5,12 @@ import functools import itertools +import logging import math import numpy as np import pytest +import tomlkit from qualtran.bloqs.phase_estimation import TextbookQPE from qhat.analysis.algorithm import build_algorithm, qpe_energy_window @@ -23,6 +25,7 @@ from qhat.common.flexible_qpe import FlexibleQPE from qhat.common.pauli_string_evolution import PauliStringEvolution from qhat.common.qpe_window_state import precision_bits +from qhat.common.tapered_qpe_sizing import fewest_runs, per_run_qpe_failure, size_tapered_qpe from qhat.common.tapered_window_state import ( BandwidthOptimizedWindowState, DPSSWindowState, TruncatedDPSSWindowState) from qhat.common.trotter_flattened import Trotterization @@ -155,15 +158,22 @@ def test_rectangular_accepts_a_split(self, unitary): assert build_algorithm(config, unitary) == FlexibleQPE(unitary, 5) @pytest.mark.parametrize("name, cls", TAPERED.items()) - def test_from_energy_error_and_probability_of_failure(self, trotter, name, cls): + @pytest.mark.parametrize("overlap", [None, 0.3, 1.0]) + def test_from_energy_error_and_probability_of_failure(self, trotter, name, cls, overlap): dE, p_fail = 0.05, 1e-6 - config = make_config(ancilla_prep=name, energy_error=dE, probability_of_failure=p_fail) + config = make_config(ancilla_prep=name, energy_error=dE, probability_of_failure=p_fail, + overlap=overlap) phase_error = dE * trotter.time / (2 * math.pi) qpe = build_algorithm(config, trotter) - assert qpe.ancilla_prep == cls.from_requirements(phase_error, p_fail) - # One fewer precision bit than the rectangular window, and far fewer confidence bits. + plan = size_tapered_qpe(cls, phase_error, p_fail, 0.5 if overlap is None else overlap) + assert qpe == FlexibleQPE( + trotter, plan.num_phase_qubits, + cls(plan.num_precision_qubits, plan.num_confidence_qubits)) + # One fewer precision bit than the rectangular window. assert qpe.ancilla_prep.num_precision_qubits == precision_bits(phase_error) - 1 - assert qpe.num_ancilla_qubits < textbook_qubits(phase_error, p_fail) + if overlap == 1.0: # a single run: the single-run rule, with far fewer confidence bits + assert qpe.ancilla_prep == cls.from_requirements(phase_error, p_fail) + assert qpe.num_ancilla_qubits < textbook_qubits(phase_error, p_fail) @pytest.mark.parametrize("method, qpe_class", [ @@ -361,7 +371,8 @@ def test_simulation_reads_exact_phases_from_pauli_file( @pytest.mark.parametrize("ancilla_prep", ["dpss", "truncated dpss", "bandwidth-optimized"]) def test_tapered_qpe_meets_requirements_through_pauli_file(tmp_path, ancilla_prep): - """End to end with a tapered window sized from (energy_error, probability_of_failure). + """End to end with a tapered window sized from (energy_error, probability_of_failure) for an + initial state that is the target eigenstate (overlap 1), so a single run. U has off-grid eigenphases on the computational basis states (see `test_simulation_reads_exact_phases_from_pauli_file`). Here energy_error 0.2 is a phase error @@ -376,7 +387,7 @@ def test_tapered_qpe_meets_requirements_through_pauli_file(tmp_path, ancilla_pre settings = dict(time=2 * math.pi / energy_scale, num_trotter_steps=3) qpe = qpe_from_pauli_file(pauli_file, ancilla_prep=ancilla_prep, energy_error=dE, - probability_of_failure=p_fail, **settings) + probability_of_failure=p_fail, overlap=1.0, **settings) num_phase = qpe.num_ancilla_qubits assert qpe.ancilla_prep.num_precision_qubits == 3 and num_phase in (7, 8) # Guard: a rectangular window of the same size misses the target, so the taper matters. @@ -399,8 +410,253 @@ def test_tapered_qpe_meets_requirements_through_pauli_file(tmp_path, ancilla_pre assert (max(failures) <= p_fail) == meets_target +# The diagonal test problem: H = -1.7 diag(PHASES), t = 2 pi / 1.7, so energies are -1.7 x PHASES +# (in [-1.7, 0)) and the ground state is basis state 3, at -1.48. energy_error 0.2 is a phase error +# of 0.118 turns. +PHASES = [0.1234, 0.377, 0.6021, 0.871] +ENERGY_SCALE = 1.7 +SETTINGS = dict(time=2 * math.pi / ENERGY_SCALE, num_trotter_steps=3) + + +@pytest.fixture +def pauli_file(tmp_path): + path = tmp_path / "hamiltonian.txt" + write_diagonal_pauli_file(path, PHASES, ENERGY_SCALE) + return path + + +def direct_multi_run_failure(n, k, p_under, p_over): + """P(the k-th smallest of n estimates is an under- or overestimate), as binomial sums: at + least k underestimates, or at least n - k + 1 overestimates.""" + def at_least(j_min, q): + return sum(math.comb(n, j) * q**j * (1 - q)**(n - j) for j in range(j_min, n + 1)) + return at_least(k, p_under) + at_least(n - k + 1, p_over) + + +@pytest.mark.parametrize("ancilla_prep, cls", TAPERED.items()) +def test_tapered_qpe_multi_run_plan_through_pauli_file(tmp_path, pauli_file, ancilla_prep, cls): + """End to end at the default overlap 0.5: the circuit is sized for N runs, reporting the K-th + smallest energy estimate, with energy_error 0.2 and p_fail 1e-3 (3 precision and 4 confidence + qubits; 14 runs for the DPSS, 15 otherwise; K = 2). + + The initial state has weight 0.5 on the ground state and 0.2, 0.3 on two excited states. From + the simulated distribution of one run's energy estimate, the N-run failure is computed + directly; it must be within the sizing's bound and p_fail. + """ + dE, p_fail = 0.2, 1e-3 + qpe = qpe_from_pauli_file(pauli_file, ancilla_prep=ancilla_prep, energy_error=dE, + probability_of_failure=p_fail, **SETTINGS) + plan = size_tapered_qpe(cls, dE / ENERGY_SCALE, p_fail, 0.5) + assert qpe.ancilla_prep == cls(3, 4) + assert (plan.num_precision_qubits, plan.num_confidence_qubits) == (3, 4) + assert plan.num_runs == (14 if cls is DPSSWindowState else 15) and plan.quantile == 2 + + m = qpe.num_ancilla_qubits + input_file = str(tmp_path / "input.npy") + save_state(input_file, np.kron(np.eye(2**m)[0], np.sqrt([0.2, 0.3, 0.0, 0.5]))) + result = numerical_simulation(input_file, qpe, compute_unitary_matrix(qpe)) + out = load_state(result["simulations"][0]["output_file"]) + probs = np.sum(np.abs(out.reshape(2**m, -1))**2, axis=1) + + # U = exp(-i H t): register value k decodes to the energy in [-1.7, 0) congruent to + # -1.7 k / 2^m. + estimates = -ENERGY_SCALE + np.mod(-ENERGY_SCALE * np.arange(2**m) / 2**m, ENERGY_SCALE) + ground = -ENERGY_SCALE * PHASES[3] + p_under = probs[estimates < ground - dE].sum() + p_over = probs[estimates > ground + dE].sum() + u = per_run_qpe_failure(cls, 3, 4) + assert p_under <= u and 1 - p_under - p_over >= 0.5 * (1 - u) + failure = direct_multi_run_failure(plan.num_runs, plan.quantile, p_under, p_over) + assert failure <= plan.multi_run_failure <= p_fail + # Guard: a single run fails about half the time. + assert p_over > 0.45 + + +@pytest.fixture +def qhat_log(caplog): + """caplog for qhat's info-level records (the "qhat" logger does not propagate to root).""" + qhat_logger = logging.getLogger("qhat") + level = qhat_logger.level + qhat_logger.addHandler(caplog.handler) + qhat_logger.setLevel(logging.INFO) + caplog.handler.setLevel(logging.INFO) + yield caplog + qhat_logger.removeHandler(caplog.handler) + qhat_logger.setLevel(level) + + +def algorithm_messages(caplog, level=logging.INFO): + return [r.getMessage() for r in caplog.records + if r.name == "qhat.analysis.algorithm" and r.levelno == level] + + +TABLE_HEADER = "-- overlap probability of failure" +NOTE = "-- NOTE: the analyses below report on a single run of the algorithm." + + +class TestMultiRunSizing: + """The tapered path of "QPE: QHAT flexible" with the multi-run sizing, through a Pauli file.""" + + def test_default_log_block(self, pauli_file, qhat_log): + qpe = qpe_from_pauli_file(pauli_file, ancilla_prep="dpss", energy_error=0.2, + probability_of_failure=1e-3, **SETTINGS) + assert qpe.ancilla_prep == DPSSWindowState(3, 4) + info = algorithm_messages(qhat_log) + start = info.index("-- overlap with target eigenstate = 0.5 (default)") + table = [f"-- {p:<9.2f} <= {f}" for p, f in zip( + [0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0], + ["9.92e-01", "5.85e-01", "1.98e-01", "4.75e-02", "8.11e-03", "9.19e-04", "6.56e-05", + "8.40e-06", "6.80e-06", "6.79e-06", "6.79e-06"])] + assert info[start:] == [ + "-- overlap with target eigenstate = 0.5 (default)", + "-- report the K-th smallest of N runs: K chosen by sizing (default)", + "-- tapered QPE sizing: 3 precision + 4 confidence qubits = 7 phase qubits", + "-- if you run the circuit 1 time:", + "-- probability of failure <= 5.00e-01, energy error <= 0.2", + "-- if you run the circuit 14 times and take the 2nd-smallest estimate:", + "-- probability of failure <= 9.19e-04, energy error <= 0.2", + "-- if the true overlap differs from 0.5 (same 14 runs, 2nd-smallest estimate; " + "energy error <= 0.2 in every row):", + TABLE_HEADER, + *table, + NOTE, + ] + assert not algorithm_messages(qhat_log, logging.WARNING) + + def test_user_specified_options(self, pauli_file, qhat_log): + qpe = qpe_from_pauli_file(pauli_file, ancilla_prep="truncated dpss", energy_error=0.2, + probability_of_failure=1e-3, overlap=0.9, + force_lowest_estimate=True, **SETTINGS) + plan = size_tapered_qpe(TruncatedDPSSWindowState, 0.2 / ENERGY_SCALE, 1e-3, 0.9, True) + assert plan.quantile == 1 + assert qpe.ancilla_prep == TruncatedDPSSWindowState(plan.num_precision_qubits, + plan.num_confidence_qubits) + info = algorithm_messages(qhat_log) + assert "-- overlap with target eigenstate = 0.9 (user-specified)" in info + assert ("-- report the smallest of N runs (K = 1): forced by " + "algorithm.force_lowest_estimate (user-specified)") in info + assert (f"-- if you run the circuit {plan.num_runs} times and take the smallest " + "estimate:") in info + + def test_force_false_is_user_specified(self, pauli_file, qhat_log): + qpe_from_pauli_file(pauli_file, ancilla_prep="dpss", energy_error=0.2, + probability_of_failure=1e-3, force_lowest_estimate=False, **SETTINGS) + assert ("-- report the K-th smallest of N runs: K chosen by sizing (user-specified)" + in algorithm_messages(qhat_log)) + + @pytest.mark.parametrize("overlap, error", [ + (0.0, ValueError), (-0.2, ValueError), (1.5, ValueError), (math.nan, ValueError), + (0.5 + 0.1j, ValueError), ("0.5", TypeError), (True, TypeError), + ]) + def test_invalid_overlap_raises(self, pauli_file, overlap, error): + with pytest.raises(error, match="overlap"): + qpe_from_pauli_file(pauli_file, ancilla_prep="dpss", energy_error=0.2, + probability_of_failure=1e-3, overlap=overlap, **SETTINGS) + + @pytest.mark.parametrize("force", ["yes", 1, None]) + def test_invalid_force_lowest_estimate_raises(self, pauli_file, force): + config = dict(ancilla_prep="dpss", energy_error=0.2, probability_of_failure=1e-3, + force_lowest_estimate=force) + if force is None: # unset is fine + qpe_from_pauli_file(pauli_file, **config, **SETTINGS) + return + with pytest.raises(TypeError, match="force_lowest_estimate"): + qpe_from_pauli_file(pauli_file, **config, **SETTINGS) + + def test_overlap_one_is_the_single_run_rule(self, pauli_file, qhat_log): + qpe = qpe_from_pauli_file(pauli_file, ancilla_prep="bandwidth-optimized", + energy_error=0.2, probability_of_failure=1e-3, overlap=1.0, + **SETTINGS) + assert qpe.ancilla_prep == BandwidthOptimizedWindowState.from_requirements( + 0.2 / ENERGY_SCALE, 1e-3) + assert "-- if you run the circuit 1 time and take the smallest estimate:" in \ + algorithm_messages(qhat_log) + + def test_fallback_when_no_plan(self, pauli_file, qhat_log): + qpe = qpe_from_pauli_file(pauli_file, ancilla_prep="dpss", energy_error=0.2, + probability_of_failure=1e-3, overlap=1e-17, **SETTINGS) + assert qpe.ancilla_prep == DPSSWindowState.from_requirements(0.2 / ENERGY_SCALE, 1e-3) + errors = algorithm_messages(qhat_log, logging.ERROR) + assert len(errors) == 1 and "FALLBACK" in errors[0] and "1e-17" in errors[0] + info = algorithm_messages(qhat_log) + assert ("-- FALLBACK tapered QPE sizing: 3 precision + 4 confidence qubits = 7 phase " + "qubits") in info + assert TABLE_HEADER not in info and info[-1] == NOTE + + @pytest.mark.parametrize("probability_of_failure", [1e-3, None]) + def test_given_counts(self, pauli_file, qhat_log, probability_of_failure): + qpe = qpe_from_pauli_file(pauli_file, ancilla_prep="dpss", num_precision_qubits=3, + num_confidence_qubits=5, overlap=0.4, + probability_of_failure=probability_of_failure, **SETTINGS) + assert qpe == qpe_from_pauli_file(pauli_file, ancilla_prep="dpss", + num_phase_qubits=8, num_confidence_qubits=5, **SETTINGS) + info = algorithm_messages(qhat_log) + assert ("-- tapered QPE: user-specified 3 precision + 5 confidence qubits = 8 phase " + "qubits") in info + assert info[-1] == NOTE + warnings = algorithm_messages(qhat_log, logging.WARNING) + if probability_of_failure is None: + assert TABLE_HEADER not in info + assert any("skipping the number of runs" in message for message in info) + assert len(warnings) == 1 and "algorithm.overlap is ignored" in warnings[0] + return + plan = fewest_runs(DPSSWindowState, 3, 5, probability_of_failure, 0.4) + # The bound comes from the precision qubits: 2^-4 turns is 1.7 / 16 = 0.106 in energy. + assert (f"-- probability of failure <= {plan.multi_run_failure:.2e}, energy error " + "<= 0.106") in info + assert (f"-- if you run the circuit {plan.num_runs} times and take the " + f"{'smallest' if plan.quantile == 1 else '2nd-smallest'} estimate:") in info + assert plan.quantile in (1, 2) + assert TABLE_HEADER in info and not warnings + + @pytest.mark.parametrize("energy_error, warns", [(0.2, False), (0.1, True)]) + def test_given_counts_warn_if_energy_error_unmet(self, pauli_file, qhat_log, energy_error, + warns): + # 3 precision qubits bound the energy error by 0.106. + qpe_from_pauli_file(pauli_file, ancilla_prep="dpss", num_precision_qubits=3, + num_confidence_qubits=4, energy_error=energy_error, + probability_of_failure=1e-3, **SETTINGS) + warnings = algorithm_messages(qhat_log, logging.WARNING) + assert len(warnings) == warns + if warns: + assert "0.106" in warnings[0] and "energy_error = 0.1" in warnings[0] + + def test_given_counts_with_no_plan(self, pauli_file, qhat_log): + # One confidence qubit: no bound on a run's failure, so no number of runs works. + qpe = qpe_from_pauli_file(pauli_file, ancilla_prep="dpss", num_precision_qubits=3, + num_confidence_qubits=1, probability_of_failure=1e-3, + **SETTINGS) + assert qpe.ancilla_prep == DPSSWindowState(3, 1) + errors = algorithm_messages(qhat_log, logging.ERROR) + assert len(errors) == 1 and "No number of runs" in errors[0] + info = algorithm_messages(qhat_log) + assert TABLE_HEADER not in info and info[-1] == NOTE + + @pytest.mark.parametrize("config", [ + dict(ancilla_prep="rectangular", energy_error=0.2, probability_of_failure=1e-3), + dict(num_phase_qubits=4), + dict(method="QPE: qualtran textbook", num_phase_qubits=4), + ], ids=["rectangular from requirements", "rectangular from counts", "textbook"]) + def test_unused_options_warn(self, pauli_file, qhat_log, config): + qpe_from_pauli_file(pauli_file, overlap=0.3, force_lowest_estimate=True, **config, + **SETTINGS) + warnings = algorithm_messages(qhat_log, logging.WARNING) + assert len(warnings) == 2 + assert "algorithm.overlap is ignored" in warnings[0] + assert "algorithm.force_lowest_estimate is ignored" in warnings[1] + assert NOTE not in algorithm_messages(qhat_log) + + class TestConfigurationAndResources: + def test_toml_saves_multi_run_options(self): + config = make_config(ancilla_prep="dpss", overlap=0.25, force_lowest_estimate=True) + table = tomlkit.parse(tomlkit.dumps({"algorithm": config._generate_TOML_table()})) + assert table["algorithm"]["overlap"] == 0.25 + assert table["algorithm"]["force_lowest_estimate"] is True + table = make_config(ancilla_prep="dpss")._generate_TOML_table() + assert "overlap" not in table and "force_lowest_estimate" not in table + def test_toml_saves_flexible_components(self): table = make_config(ancilla_prep="rectangular", qft_inv="textbook")._generate_TOML_table() assert table["ancilla_prep"] == "rectangular" diff --git a/common/tapered_qpe_sizing.py b/common/tapered_qpe_sizing.py new file mode 100644 index 000000000..30444926b --- /dev/null +++ b/common/tapered_qpe_sizing.py @@ -0,0 +1,274 @@ +""" +Sizing tapered QPE (see `qhat.common.tapered_window_state`) when the initial state only partly +overlaps the target eigenstate: how many precision and confidence qubits to use, how many times N +to run the circuit, and which estimate to report (the K-th smallest of the N). This implements the +multi-run analysis by Yigit Subasi and Alan Bu (unpublished). + +One run returns one phase estimate. With per-run QPE failure u and overlap p between the initial +state and the target eigenstate, the estimate is + +- accepted (within 2^-(num_precision_qubits + 1) turns of the target phase) with probability at + least p (1 - u); +- an underestimate with probability at most u. Nothing lies below the target, so only a failed + QPE gives an underestimate; +- otherwise an overestimate. + +So one run fails with probability at most 1 - p (1 - u). Reporting the K-th smallest of N runs +fails if that estimate is an under- or overestimate; `multi_run_failure` bounds this over every +split of the failure probability allowed above. + +Example: phase error 0.2/1.7 turns, probability of failure 1e-3, overlap 0.5, DPSS window: 3 +precision and 4 confidence qubits; one run fails with probability up to 0.50, and the 2nd-smallest +of 14 runs fails with probability up to 9.2e-4. + +The per-run failure u is `failure_bound_factor * dpss_average_failure_bound`, an average over the +offset of the true phase from the QPE grid. The paper's random phase shift (arXiv:2403.18927, +Sec. VI) would make it a worst-case bound; qhat does not apply that shift. +""" + +import math +from typing import Callable, Optional, Tuple, Type + +import attrs +from scipy.special import betainc, betaincc + +from qhat.common.qpe_window_state import precision_bits +from qhat.common.tapered_window_state import TaperedWindowStateBase, dpss_average_failure_bound + +# Confidence qubits tried by `size_tapered_qpe`. +MIN_CONFIDENCE_QUBITS = 2 +MAX_CONFIDENCE_QUBITS = 12 +# Largest cost N * 2^(number of phase qubits) considered. +COST_CEILING = 1e18 +# With no answer yet, stop trying larger K after this many in a row fail to get closer. +_QUANTILE_PATIENCE = 24 +# Larger K tried after the first one that works, in case one of them needs fewer runs. +_QUANTILE_GUARD = 2 + + +@attrs.frozen +class TaperedQPESizing: + """Qubit counts and run plan. E.g. 3 precision and 4 confidence qubits, run 14 times and report + the 2nd-smallest estimate (`quantile` 2): `single_run_failure` 0.50, `multi_run_failure` + 9.2e-4 at `overlap` 0.5. + """ + + num_precision_qubits: int + num_confidence_qubits: int + num_runs: int + quantile: int + single_run_failure: float + multi_run_failure: float + overlap: float + force_lowest_estimate: bool + + @property + def num_phase_qubits(self) -> int: + return self.num_precision_qubits + self.num_confidence_qubits + + +def per_run_qpe_failure(window_cls: Type[TaperedWindowStateBase], num_precision_qubits: int, + num_confidence_qubits: int) -> float: + """u = `failure_bound_factor * dpss_average_failure_bound(...)`, at most 1. E.g. the DPSS with + 3 precision and 4 confidence qubits: 2.7e-4.""" + return min(1.0, window_cls.failure_bound_factor + * dpss_average_failure_bound(num_precision_qubits, num_confidence_qubits)) + + +def multi_run_failure(window_cls: Type[TaperedWindowStateBase], num_precision_qubits: int, + num_confidence_qubits: int, num_runs: int, quantile: int, + overlap: float) -> float: + """Bound on the probability that the `quantile`-th smallest of `num_runs` estimates misses the + target phase by more than 2^-(num_precision_qubits + 1) turns. E.g. the DPSS with 3 + 4 + qubits, 14 runs, quantile 2: 9.2e-4 at overlap 0.5, 6.8e-6 at overlap 1. + """ + u = per_run_qpe_failure(window_cls, num_precision_qubits, num_confidence_qubits) + return _worst_failure(num_runs, quantile, u, overlap * (1 - u)) + + +def fewest_runs(window_cls: Type[TaperedWindowStateBase], num_precision_qubits: int, + num_confidence_qubits: int, probability_of_failure: float, overlap: float, + force_lowest_estimate: bool = False) -> Optional[TaperedQPESizing]: + """The fewest runs N (and the best K for that N; K = 1 if `force_lowest_estimate`) that meet + `probability_of_failure` with these qubit counts, or None if none does within + `COST_CEILING`. + """ + _check_targets(probability_of_failure, overlap) + return _plan(window_cls, num_precision_qubits, num_confidence_qubits, probability_of_failure, + overlap, force_lowest_estimate, COST_CEILING) + + +def size_tapered_qpe(window_cls: Type[TaperedWindowStateBase], phase_error: float, + probability_of_failure: float, overlap: float, + force_lowest_estimate: bool = False) -> Optional[TaperedQPESizing]: + """The cheapest plan meeting `phase_error` (turns) and `probability_of_failure` (for the + reported estimate of N runs), or None if none costs at most `COST_CEILING`. + + Precision qubits: the fewest with 2^-(num_precision_qubits + 1) <= phase_error. Confidence + qubits (MIN_CONFIDENCE_QUBITS to MAX_CONFIDENCE_QUBITS), N and K minimize the cost + N * 2^(number of phase qubits); ties go to fewer qubits. K = 1 if `force_lowest_estimate`. + """ + _check_targets(probability_of_failure, overlap) + num_precision = precision_bits(phase_error) - 1 + best, limit = None, COST_CEILING + for num_confidence in range(MIN_CONFIDENCE_QUBITS, MAX_CONFIDENCE_QUBITS + 1): + if 2.0**(num_precision + num_confidence) > limit: + break + plan = _plan(window_cls, num_precision, num_confidence, probability_of_failure, overlap, + force_lowest_estimate, limit) + if plan is None: + continue + cost = plan.num_runs * 2.0**plan.num_phase_qubits + if best is None or cost < limit: + best, limit = plan, cost + return best + + +def _check_targets(probability_of_failure: float, overlap: float): + if not 0 < probability_of_failure < 1: + raise ValueError( + f"probability of failure must be in (0, 1), got {probability_of_failure}.") + if not 0 < overlap <= 1: + raise ValueError(f"overlap must be in (0, 1], got {overlap}.") + + +def _plan(window_cls, num_precision, num_confidence, target, overlap, force_lowest_estimate, + cost_limit) -> Optional[TaperedQPESizing]: + u = per_run_qpe_failure(window_cls, num_precision, num_confidence) + hit = overlap * (1 - u) # smallest probability that a run is accepted + if hit <= 0: + return None + cap = int(cost_limit // 2.0**(num_precision + num_confidence)) + # Even K = 1 fails if every run overestimates, which has probability up to (1 - hit)^N. + min_runs = math.log(target) / math.log1p(-hit) + if min_runs > cap: + return None + found = _fewest_runs(u, hit, target, force_lowest_estimate, + max(1, math.ceil(min_runs - 1e-9)), cap) + if found is None: + return None + runs, quantile = found + return TaperedQPESizing(num_precision, num_confidence, runs, quantile, 1 - hit, + _worst_failure(runs, quantile, u, hit), overlap, + force_lowest_estimate) + + +def _failure(n: int, k: int, under: float, not_over: float) -> float: + """P(the k-th smallest of n estimates is an under- or overestimate) when each run independently + underestimates with probability `under` and overestimates with probability 1 - `not_over`. + + Underestimate: at least k runs underestimate. Overestimate: at most k - 1 runs do not + overestimate. Both terms are computed directly, so a small failure keeps its precision. + """ + low = betainc(k, n - k + 1, under) + high = betaincc(k, n - k + 1, not_over) + return min(1.0, float(low + high)) + + +def _worst_failure(n: int, k: int, u: float, hit: float) -> float: + """`_failure` at the worse of the two extreme splits of the single-run failure 1 - hit: all + overestimates, or underestimates u and the rest overestimates. The worst split is one of the + two.""" + return max(_failure(n, k, 0.0, hit), _failure(n, k, u, min(1.0, u + hit))) + + +def _first(pred: Callable[[int], bool], lo: int, hi: int) -> int: + """Smallest x in [lo, hi] with pred(x), for pred false then true, and pred(hi) true.""" + while lo < hi: + mid = (lo + hi) // 2 + if pred(mid): + hi = mid + else: + lo = mid + 1 + return lo + + +def _last(pred: Callable[[int], bool], lo: int, hi: int) -> int: + """Largest x in [lo, hi] with pred(x), for pred true then false, and pred(lo) true.""" + while lo < hi: + mid = (lo + hi + 1) // 2 + if pred(mid): + lo = mid + else: + hi = mid - 1 + return lo + + +def _argmin(f: Callable[[int], float], lo: int, hi: int) -> int: + """x in [lo, hi] minimizing f, for f falling then rising.""" + while hi - lo > 2: + m1, m2 = lo + (hi - lo) // 3, hi - (hi - lo) // 3 + if f(m1) <= f(m2): + hi = m2 + else: + lo = m1 + return min(range(lo, hi + 1), key=f) + + +def _runs_for_quantile(k: int, u: float, hit: float, target: float, n_min: int, + cap: int) -> Tuple[Optional[int], float]: + """Fewest runs n in [max(k, n_min), cap] with `_worst_failure(n, k, ...) <= target`, or None; + also the smallest worst failure seen, to tell when larger k stop helping. + + The all-overestimate split only improves with more runs, which allows n >= n_over. In the other + split, more runs mean more chances to underestimate: its underestimate term alone allows + n <= n_under, and between n_over and n_under its failure falls, then rises. + """ + lo = max(k, n_min) + if lo > cap: + return None, 1.0 + if _failure(cap, k, 0.0, hit) > target: + return None, _failure(cap, k, 0.0, hit) + n_over = _first(lambda n: _failure(n, k, 0.0, hit) <= target, lo, cap) + split = lambda n: _failure(n, k, u, min(1.0, u + hit)) + if _failure(n_over, k, u, 1.0) > target: + return None, _worst_failure(n_over, k, u, hit) + n_under = _last(lambda n: _failure(n, k, u, 1.0) <= target, n_over, cap) + best = _argmin(split, n_over, n_under) + if split(best) > target: + return None, split(best) + return _first(lambda n: split(n) <= target, n_over, best), split(best) + + +def _fewest_runs(u: float, hit: float, target: float, force_lowest_estimate: bool, n_min: int, + cap: int) -> Optional[Tuple[int, int]]: + """(N, K) with the fewest runs N in [n_min, cap] meeting `target`, and the K with the smallest + failure at that N (K = 1 if `force_lowest_estimate`); None if no N works.""" + if force_lowest_estimate: + runs, _ = _runs_for_quantile(1, u, hit, target, n_min, cap) + return None if runs is None else (runs, 1) + + best = None + nearest, stale, extra = math.inf, 0, 0 + k = 1 + # Runs for quantile k number at least k, so once k reaches the best N nothing can beat it. + while k <= cap and (best is None or k < best[0]): + # The all-overestimate failure rises with k, so if it misses at cap, so do larger k. + if _failure(cap, k, 0.0, hit) > target: + break + runs, seen = _runs_for_quantile(k, u, hit, target, n_min, cap) + if runs is not None and (best is None or runs < best[0]): + best = (runs, k) + if best is None: + nearest, stale = (seen, 0) if seen < nearest else (nearest, stale + 1) + if stale >= _QUANTILE_PATIENCE: + break + else: + extra += 1 + if extra > _QUANTILE_GUARD: + break + k += 1 + if best is None: + return None + runs, k = best + return runs, _best_quantile(runs, k, u, hit, target) + + +def _best_quantile(n: int, k_found: int, u: float, hit: float, target: float) -> int: + """K with the smallest worst failure at n runs. Searched where each split's dominant term + meets the target (k_found does), so the failure falls, then rises.""" + k_lo = _first(lambda k: _failure(n, k, u, 1.0) <= target, 1, k_found) + k_hi = _last(lambda k: _failure(n, k, 0.0, hit) <= target, k_found, n) + worst = lambda k: _worst_failure(n, k, u, hit) + k_best = _argmin(worst, k_lo, k_hi) + return k_best if worst(k_best) < worst(k_found) else k_found diff --git a/common/tests/test_tapered_qpe_sizing.py b/common/tests/test_tapered_qpe_sizing.py new file mode 100644 index 000000000..771e89069 --- /dev/null +++ b/common/tests/test_tapered_qpe_sizing.py @@ -0,0 +1,309 @@ +""" +Tests for the multi-run sizing of tapered QPE (`qhat.common.tapered_qpe_sizing`). + +Reference values for the default problem come from the reference implementation of the +multi-run analysis by Yigit Subasi and Alan Bu. At very large N that implementation loses digits; +there the values are checked against 60-digit mpmath instead. +""" + +import itertools +import math + +import mpmath +import numpy as np +import pytest +from scipy.stats import binom + +from qhat.common.tapered_qpe_sizing import ( + COST_CEILING, + fewest_runs, + MAX_CONFIDENCE_QUBITS, + MIN_CONFIDENCE_QUBITS, + multi_run_failure, + per_run_qpe_failure, + size_tapered_qpe, +) +from qhat.common.tapered_window_state import ( + BandwidthOptimizedWindowState, + DPSSWindowState, + TruncatedDPSSWindowState, +) + +TAPERED = [DPSSWindowState, TruncatedDPSSWindowState, BandwidthOptimizedWindowState] +TAPERED_IDS = ["dpss", "truncated", "bandwidth-optimized"] + +# Energy error 0.2 on a spectrum scaled by 1.7: 3 precision qubits. +PHASE_ERROR = 0.2 / 1.7 +# Phase error for 10 precision qubits, as qhat's driver sets it (phase scale factor 1.01). +FINE_PHASE_ERROR = 1 / (1.01 * 2**10) + + +def direct_failure(n, k, under, over): + """P(the k-th smallest of n runs is an under- or overestimate), summed term by term.""" + def at_least(j_min, q): + return sum(math.comb(n, j) * q**j * (1 - q)**(n - j) for j in range(j_min, n + 1)) + # Underestimate: at least k runs underestimate. Overestimate: at least n - k + 1 overestimate. + return at_least(k, under) + at_least(n - k + 1, over) + + +def mpmath_failure(window_cls, num_precision, num_confidence, n, k, overlap): + """`multi_run_failure` in 60-digit arithmetic.""" + with mpmath.workdps(60): + u = mpmath.mpf(per_run_qpe_failure(window_cls, num_precision, num_confidence)) + hit = mpmath.mpf(overlap) * (1 - u) + + def failure(under, not_over): + return (mpmath.betainc(k, n - k + 1, 0, under, regularized=True) + + mpmath.betainc(k, n - k + 1, not_over, 1, regularized=True)) + return max(failure(0, hit), failure(u, u + hit)) + + +class TestMultiRunFailure: + @pytest.mark.parametrize("n, k", [(1, 1), (5, 2), (14, 2), (14, 9), (30, 30)]) + @pytest.mark.parametrize("overlap", [0.05, 0.5, 1.0]) + def test_matches_direct_sum(self, n, k, overlap): + u = per_run_qpe_failure(DPSSWindowState, 1, 3) + hit = overlap * (1 - u) + expected = max(direct_failure(n, k, 0.0, 1 - hit), + direct_failure(n, k, u, 1 - u - hit)) + assert multi_run_failure(DPSSWindowState, 1, 3, n, k, overlap) == pytest.approx( + expected, rel=1e-10, abs=1e-15) + + @pytest.mark.parametrize("n, k", [(7, 1), (14, 2), (25, 6)]) + def test_worst_split_is_an_endpoint(self, n, k): + # The adversary may split the single-run failure 1 - hit into underestimates (at most u) + # and overestimates however it likes; no split does worse than the two extremes. + u = per_run_qpe_failure(DPSSWindowState, 0, 3) + overlap = 0.4 + hit = overlap * (1 - u) + bound = multi_run_failure(DPSSWindowState, 0, 3, n, k, overlap) + for under in np.linspace(0, u, 21): + assert direct_failure(n, k, under, 1 - hit - under) <= bound * (1 + 1e-12) + + def test_single_run_at_overlap_one_is_the_qpe_failure(self): + u = per_run_qpe_failure(DPSSWindowState, 3, 4) + assert multi_run_failure(DPSSWindowState, 3, 4, 1, 1, 1.0) == pytest.approx(u, rel=1e-12) + + def test_overlap_table(self): + # The plan sized for overlap 0.5 (3 + 4 qubits, 2nd-smallest of 14 runs), if the true + # overlap is different. + overlaps = [0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] + expected = [9.92e-1, 5.85e-1, 1.98e-1, 4.75e-2, 8.11e-3, 9.19e-4, 6.56e-5, 8.40e-6, + 6.80e-6, 6.79e-6, 6.79e-6] + got = [multi_run_failure(DPSSWindowState, 3, 4, 14, 2, p) for p in overlaps] + assert got == pytest.approx(expected, rel=5e-3) + + @pytest.mark.parametrize("n, k, overlap", [(1383365, 1, 1e-5), (239448133, 2, 1e-7), + (3367460167391, 6, 1e-11)]) + def test_large_n_matches_high_precision(self, n, k, overlap): + expected = float(mpmath_failure(DPSSWindowState, 10, 6, n, k, overlap)) + assert multi_run_failure(DPSSWindowState, 10, 6, n, k, overlap) == pytest.approx( + expected, rel=1e-9) + + +class TestSizeTaperedQPE: + @pytest.mark.parametrize("overlap, runs, quantile, failure", [ + (1.0, 1, 1, 2.735e-4), + (0.9, 5, 2, 4.644e-4), + (0.5, 14, 2, 9.193e-4), + (0.25, 33, 2, 9.361e-4), + (0.1, 92, 2, 9.850e-4), + ]) + def test_reference_values(self, overlap, runs, quantile, failure): + plan = size_tapered_qpe(DPSSWindowState, PHASE_ERROR, 1e-3, overlap) + assert (plan.num_precision_qubits, plan.num_confidence_qubits) == (3, 4) + assert (plan.num_runs, plan.quantile) == (runs, quantile) + assert plan.multi_run_failure == pytest.approx(failure, rel=1e-3) + assert plan.single_run_failure == pytest.approx( + 1 - overlap * (1 - per_run_qpe_failure(DPSSWindowState, 3, 4)), rel=1e-12) + assert plan.overlap == overlap + assert plan.force_lowest_estimate is False + + def test_reference_preset_best_and_lowest_estimate(self): + # Energy error 2^-10 on a unit spectrum padded by the error on both sides; 95% confidence. + delta = 2**-10 + phase_error = delta / (1 + 2 * delta) + best = size_tapered_qpe(DPSSWindowState, phase_error, 0.05, 0.5) + assert (best.num_runs, best.quantile) == (8, 2) + assert (best.num_precision_qubits, best.num_confidence_qubits) == (10, 4) + assert 1 - best.multi_run_failure == pytest.approx(0.954553, abs=1e-6) + lowest = size_tapered_qpe(DPSSWindowState, phase_error, 0.05, 0.5, + force_lowest_estimate=True) + assert (lowest.num_runs, lowest.quantile) == (5, 1) + assert (lowest.num_precision_qubits, lowest.num_confidence_qubits) == (10, 5) + assert lowest.force_lowest_estimate is True + + @pytest.mark.parametrize("window_cls", TAPERED, ids=TAPERED_IDS) + @pytest.mark.parametrize("phase_error, probability_of_failure", + [(PHASE_ERROR, 1e-3), (0.03, 1e-8), (FINE_PHASE_ERROR, 0.2)]) + def test_overlap_one_matches_single_run_rule(self, window_cls, phase_error, + probability_of_failure): + plan = size_tapered_qpe(window_cls, phase_error, probability_of_failure, 1.0) + window = window_cls.from_requirements(phase_error, probability_of_failure) + assert (plan.num_runs, plan.quantile) == (1, 1) + assert plan.num_precision_qubits == window.num_precision_qubits + assert plan.num_confidence_qubits == window.num_confidence_qubits + + def test_band_limited_needs_at_least_as_much(self): + dpss = size_tapered_qpe(DPSSWindowState, PHASE_ERROR, 1e-6, 0.3) + for window_cls in TAPERED[1:]: + other = size_tapered_qpe(window_cls, PHASE_ERROR, 1e-6, 0.3) + assert (other.num_runs * 2**other.num_phase_qubits + >= dpss.num_runs * 2**dpss.num_phase_qubits) + + @pytest.mark.parametrize("overlap, probability_of_failure", [(0.5, 1e-3), (0.07, 1e-5)]) + @pytest.mark.parametrize("force", [False, True]) + def test_cheapest_over_confidence_qubits(self, overlap, probability_of_failure, force): + plan = size_tapered_qpe(TruncatedDPSSWindowState, PHASE_ERROR, probability_of_failure, + overlap, force) + cost = plan.num_runs * 2**plan.num_phase_qubits + for num_confidence in range(MIN_CONFIDENCE_QUBITS, MAX_CONFIDENCE_QUBITS + 1): + other = fewest_runs(TruncatedDPSSWindowState, 3, num_confidence, + probability_of_failure, overlap, force) + if other is not None: + assert other.num_runs * 2**other.num_phase_qubits >= cost + + @pytest.mark.parametrize("overlap, probability_of_failure, runs, quantile", [ + (1e-5, 1e-6, 1383365, 1), + (1e-5, 1e-9, 2393962, 2), + (1e-7, 1e-9, 239448133, 2), + (1e-11, 1e-9, 3367460167391, 6), + ]) + def test_large_n(self, overlap, probability_of_failure, runs, quantile): + plan = size_tapered_qpe(DPSSWindowState, FINE_PHASE_ERROR, probability_of_failure, + overlap) + assert (plan.num_precision_qubits, plan.num_confidence_qubits) == (10, 6) + assert (plan.num_runs, plan.quantile) == (runs, quantile) + assert mpmath_failure(DPSSWindowState, 10, 6, runs, quantile, + overlap) <= probability_of_failure + assert mpmath_failure(DPSSWindowState, 10, 6, runs - 1, quantile, + overlap) > probability_of_failure + + def test_infeasible_within_cost_ceiling(self): + # Even K = 1 needs about 2e15 runs of at least 2^12 queries each. + assert size_tapered_qpe(DPSSWindowState, FINE_PHASE_ERROR, 1e-9, 1e-14) is None + assert size_tapered_qpe(DPSSWindowState, FINE_PHASE_ERROR, 1e-9, 1e-14, True) is None + + @pytest.mark.parametrize("overlap", [0.0, -0.1, 1.5]) + def test_rejects_overlap(self, overlap): + with pytest.raises(ValueError, match="overlap"): + size_tapered_qpe(DPSSWindowState, PHASE_ERROR, 1e-3, overlap) + + @pytest.mark.parametrize("probability_of_failure", [0.0, 1.0]) + def test_rejects_probability_of_failure(self, probability_of_failure): + with pytest.raises(ValueError, match="probability of failure"): + size_tapered_qpe(DPSSWindowState, PHASE_ERROR, probability_of_failure, 0.5) + + +class TestFewestRuns: + @pytest.mark.parametrize("window_cls", TAPERED, ids=TAPERED_IDS) + # Answers range from N = 6 to 57, K = 1 to 17; forcing K = 1 is impossible in some, and with 2 + # confidence qubits nothing works. + @pytest.mark.parametrize("num_precision, num_confidence, overlap, probability_of_failure", [ + (1, 3, 0.6, 1e-2), + (2, 4, 0.6, 1e-2), + (0, 4, 0.8, 1e-4), + (1, 4, 0.3, 1e-3), + (2, 3, 0.5, 0.05), + (3, 2, 0.9, 0.3), + ]) + @pytest.mark.parametrize("force", [False, True]) + def test_matches_exhaustive_search(self, window_cls, num_precision, num_confidence, overlap, + probability_of_failure, force): + def failure(n, k): + return multi_run_failure(window_cls, num_precision, num_confidence, n, k, overlap) + + expected = None + for n in range(1, 120): + quantiles = [1] if force else range(1, n + 1) + if any(failure(n, k) <= probability_of_failure for k in quantiles): + expected = (n, min(quantiles, key=lambda k: failure(n, k))) + break + plan = fewest_runs(window_cls, num_precision, num_confidence, probability_of_failure, + overlap, force) + if expected is None: + assert plan is None + return + assert plan.num_runs == expected[0] + assert failure(plan.num_runs, plan.quantile) == pytest.approx( + failure(*expected), rel=1e-9) + assert plan.multi_run_failure <= probability_of_failure + + def test_no_runs_suffice(self): + # One confidence qubit: the per-run failure bound is 1, so no run is ever trusted. + assert fewest_runs(DPSSWindowState, 3, 1, 1e-3, 0.5) is None + + def test_respects_cost_ceiling(self): + num_precision, num_confidence = 10, 6 + plan = fewest_runs(DPSSWindowState, num_precision, num_confidence, 1e-9, 1e-11) + assert plan.num_runs * 2**(num_precision + num_confidence) <= COST_CEILING + assert fewest_runs(DPSSWindowState, num_precision, num_confidence, 1e-9, 1e-14) is None + + +def outcome_probs(window, phase): + """QPE outcome probabilities for a window state and a phase (turns): |QFT^-1 (w e^{2 pi i n + phase})|^2.""" + n = np.arange(len(window)) + return np.abs(np.fft.fft(window * np.exp(2j * np.pi * n * phase)))**2 / len(window) + + +def simulated_multi_run_failure(window, phases, weights, phase_error, num_runs, quantile): + """Exact failure of reporting the `quantile`-th smallest of `num_runs` runs, for an initial + state with weights `weights[j]` on eigenphases `phases[j]`; phases[0] is the target (lowest). + """ + num_outcomes = len(window) + probs = sum(c * outcome_probs(window, phase) for c, phase in zip(weights, phases)) + estimates = np.arange(num_outcomes) / num_outcomes + distance = np.abs(estimates - phases[0]) + accepted = np.minimum(distance, 1 - distance) <= phase_error + under = ~accepted & (estimates < phases[0]) + p_under, p_over = probs[under].sum(), probs[~accepted & ~under].sum() + # At least `quantile` underestimates, or at least num_runs - quantile + 1 overestimates. + return (binom.sf(quantile - 1, num_runs, p_under) + + binom.sf(num_runs - quantile, num_runs, p_over)) + + +class TestModelAgainstQPE: + """The multi-run bound against exact QPE outcome distributions for random spectra and initial + states, with every eigenphase at least the phase error from the window edges.""" + + @pytest.mark.parametrize("window_cls", TAPERED, ids=TAPERED_IDS) + def test_bound_holds_for_random_spectra(self, window_cls): + rng = np.random.default_rng(7) + for _ in range(40): + num_precision = int(rng.integers(1, 4)) + num_confidence = int(rng.integers(2, 5)) + overlap = float(rng.uniform(0.05, 1.0)) + probability_of_failure = 10**rng.uniform(-6, -1) + plan = fewest_runs(window_cls, num_precision, num_confidence, probability_of_failure, + overlap) + if plan is None: + continue + phase_error = 2.0**-(num_precision + 1) + num_excited = int(rng.integers(0, 6)) + phases = np.sort(rng.uniform(phase_error, 1 - phase_error, num_excited + 1)) + weights = np.concatenate( + [[overlap], rng.dirichlet(np.ones(num_excited)) * (1 - overlap) + if num_excited else []]) + if num_excited == 0: + weights = np.array([1.0]) + window = window_cls(num_precision, num_confidence).taper() + failure = simulated_multi_run_failure(window, phases, weights, phase_error, + plan.num_runs, plan.quantile) + assert failure <= plan.multi_run_failure * (1 + 1e-9) + + def test_bound_needs_margin_below_phase_one(self): + # An excited state within the phase error of phase 1: its correct estimates wrap to + # near 0 and count as underestimates, so the bound fails by a factor of about 27. This + # is the window-edge margin the driver does not yet guarantee. + num_precision, num_confidence, overlap = 3, 4, 0.5 + phase_error = 2.0**-(num_precision + 1) + plan = fewest_runs(DPSSWindowState, num_precision, num_confidence, 1e-3, overlap) + window = DPSSWindowState(num_precision, num_confidence).taper() + weights = [overlap, 1 - overlap] + inside = simulated_multi_run_failure(window, [0.3, 0.6], weights, phase_error, + plan.num_runs, plan.quantile) + near_edge = simulated_multi_run_failure(window, [0.3, 1 - 0.3 * phase_error], weights, + phase_error, plan.num_runs, plan.quantile) + assert inside <= plan.multi_run_failure + assert near_edge > 20 * plan.multi_run_failure