tapering QPE parameters

analysis/README.md

@@ -320,10 +320,44 @@ 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.
 
 ### Analyzing an Algorithm
 

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",):

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
 
 # -------------------------------------------------------------------------------------------------

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