open-standards/wai/quantum-cross-simulator/tools/tolerance.py

"""Every tolerance the corpus holds an engine to, each derived from that engine's arithmetic.

None is fitted to a result. Each is a first-order worst-case bound, and so is loose: an engine
that passes is within it, and the run also records how far inside it each engine is. Where an
assumption enters (a library's backward-error constant, which the library does not state), the
function says so and states the number it takes.

Floating point (unit roundoff u: 2^-53 in double, 2^-24 in single). One gate acting on w qubits
at once, applied to a normalised state, errs in the Euclidean norm by at most C(w)·u:

  C(w) = (2√2 + (2^w − 1) + 2 + π) · 2^(w/2)

- 2√2 u: a complex product of two rounded numbers (Higham, Accuracy and Stability of Numerical
  Algorithms, 2nd ed., §3.6);
- 2^w − 1: the additions of a row of 2^w terms, summed in any order;
- 2 + π: the matrix entries, each within 2u from its sine and cosine, and moved by at most π·u
  by the rounding of an angle of at most 2π;
- 2^(w/2): the norm of the matrix of absolute values, at most √(2^w) for a unitary.

Errors of successive gates add (each gate is unitary, so it carries the error before it unchanged
in norm). A fused gate of width w costs C(w) for each gate fused into it, since its matrix is a
product of theirs computed in the same arithmetic.
"""
import math

U_DOUBLE = 2.0**-53
U_SINGLE = 2.0**-24


def C(w):
    return (2 * math.sqrt(2) + (2**w - 1) + 2 + math.pi) * 2 ** (w / 2)


def gate_count(steps):
    n = 0
    for s in steps:
        for x in s.get("then", [s]):
            if "g" in x:
                n += 1
    return n


def _first_order(total):
    # The bound's own second-order term: (1 + x)·x for x below 1.
    return total * (1 + total) if total < 1 else float("inf")


def double_state(gates, n, width):
    return _first_order(gates * C(width) * U_DOUBLE)


def single_state(gates, n, width):
    """Single precision (complex64 amplitudes) with the same model and u = 2^-24."""
    return _first_order(gates * C(width) * U_SINGLE)


def contraction(steps, width):
    """A tensor network contracted pairwise: a contraction with an index of 2^width values is a
    row sum of that many terms, so each step costs C(width). The model assumes, as for a gate, that
    the norm of an intermediate's matrix of absolute values is at most √(2^width) times its own;
    for a network of unitaries contracted in a sensible order that holds, and the measured
    errors (recorded beside it) are the check on the assumption."""
    return _first_order(steps * C(width) * U_DOUBLE)


DD_TOLERANCE = 2.0**-52 * 1024  # mqt.core's dd::RealNumber::eps, read from its header


def decision_diagram(gates, n):
    """A decision-diagram engine merges edge weights that lie within its tolerance τ of an entry
    of its table, so each weight moves by at most τ; an amplitude is a product of n + 1 weights,
    so each gate moves it by at most (n + 1)·τ relatively, and the double arithmetic adds the
    state vector's bound."""
    return _first_order(gates * (n + 1) * DD_TOLERANCE) + double_state(gates, n, 2)


LAPACK_CONSTANT = 16  # p(m, n) in LAPACK's SVD backward error p(m, n)·u·||A||, taken as 16·m


def svd_steps(steps):
    """The two-site updates a matrix product state makes for a circuit: a gate on qubits d apart is
    routed by d − 1 swaps there and d − 1 back, each an update, and the gate itself one more."""
    total = 0
    for s in steps:
        for x in s.get("then", [s]):
            if "g" in x and len(x["q"]) >= 2:
                total += 2 * (max(x["q"]) - min(x["q"])) - 1
    return total


def mps(updates, chi, truncation):
    """A matrix product state: each two-site update is a contraction and an SVD of a 2χ by 2χ
    matrix, whose backward error LAPACK states as p(m, n)·u with p "a modestly growing function"
    of the shape; we take p = 16·2χ. `truncation` bounds what the engine's truncations moved the
    state: a truncation that drops weight w moves a normalised state by at most 2√w (√w dropped,
    at most w more by renormalising), and the unitary steps after it carry that unchanged, so
    these add. The weights are the ones the engine reports (truncation_from)."""
    svd = updates * (C(2) + LAPACK_CONSTANT * 2 * chi) * U_DOUBLE
    return _first_order(svd) + truncation


def truncation_from(weights=None, total=None, count=0):
    """2·Σ√wᵢ from each truncation's weight when the engine reports them one by one; from a total
    W over `count` truncations, 2·√(count·W), which bounds that sum (Cauchy–Schwarz)."""
    if weights is not None:
        return 2 * sum(math.sqrt(max(w, 0.0)) for w in weights)
    return 2 * math.sqrt(max(count, 0) * max(total or 0.0, 0.0))


def density(gates, n):
    """A density matrix: each gate is applied on both sides (2·C(2)), each depolarizing channel
    is a sum of up to 16 Pauli conjugations (16·C(2)); an expectation value of a Pauli string
    errs by at most √(2^n) times the error of ρ in the Frobenius norm."""
    return math.sqrt(2**n) * gates * (2 + 16) * C(2) * U_DOUBLE


def dem(noise_sites):
    """A detector error model's probabilities. A channel's per-Pauli rate, (1 − √(1 − 4p/3))/2 or
    its two-qubit form, is within 3u absolutely; merging two mechanisms, p + q − 2pq, adds at most
    3u more; a mechanism merges at most every noise channel of the circuit. Two engines differ by at
    most twice that: 2 · noise_sites · 6u, written 16u per site."""
    return noise_sites * 16 * U_DOUBLE


# Fixed point: the extension's own bound (quantum-sim.md §5). Per inexact gate, every component
# it writes is within 1 + √2/2 units of 2^-30 of the ideal gate applied to the same values: one
# unit from two products rounded half a unit each, √2/2 from table constants within half a unit.
PER_COMPONENT = 1 + math.sqrt(2) / 2
UNIT = 2.0**-30


def fixed_point_op(n, base, controls, param):
    """The Euclidean bound, in units, for one op of the bit-exact op-log: base gate codes are
    quantum-sim.md's (§3.0): I X Y Z H S S† T T† P = 0..9."""
    c = len(controls)
    if base in (0, 1, 2, 3, 5, 6):
        return 0.0  # entries 0, ±1, ±i: exact
    if base == 9 and param <= 2:
        return 0.0  # P(1) = -1 and P(2) = i: exact
    if base == 4:
        # H writes both amplitudes of each pair: four components per pair.
        return PER_COMPONENT * math.sqrt(4 * 2 ** (n - 1 - c))
    # T, T† and P(k) write one amplitude of each pair: two components.
    return PER_COMPONENT * math.sqrt(2 * 2 ** (n - 1 - c))


def fixed_point_bound(n, ops):
    """The bit-exact method's distance from the ideal state, bounded by its published per-gate
    bound summed over the op-log it ran, each gate carrying the error before it: a fixed-point
    gate's matrix is within 2·2^-30 of a unitary, so it can grow that error by at most that."""
    total = 0.0
    inexact = 0
    for base, controls, target, param in ops:
        b = fixed_point_op(n, base, controls, param)
        if b:
            inexact += 1
        total += b
    return total * UNIT * (1 + 2 * UNIT) ** inexact


def l1_distance(shots, outcomes, alpha):
    """ε with P(Σ|X_k/shots − p_k| ≥ ε) ≤ alpha for counts X of `shots` draws over `outcomes`
    categories: Weissman, Ordentlich, Seroussi, Verdú and Weinberger (2003), (2^K − 2)·exp(−shots·ε²/2).
    The L1 distance is twice the total variation. No normal approximation. It holds the whole
    distribution at once, where the per-outcome test below holds each count alone."""
    return math.sqrt(2 * (outcomes * math.log(2) + math.log(1 / alpha)) / shots)


def bernstein(shots, p, alpha):
    """t with P(|X − shots·p| ≥ t) ≤ alpha for X binomial(shots, p), from Bernstein's
    inequality 2·exp(−t² / (2(V + t/3))), V = shots·p·(1 − p). No normal approximation."""
    b = math.log(2 / alpha)
    v = shots * p * (1 - p)
    return b / 3 + math.sqrt(b * b / 9 + 2 * b * v)