# Variational algorithms and the algorithm library

`engine::variational` and `engine::algorithms` answer five calls, the same on
every surface: the command line (`wai-quantum gradient`, `optimize`, `vqe`,
`qaoa`, `build`), both MCP servers (`gradient_qasm`, `optimize_qasm`,
`vqe_ground_state`, `qaoa_maxcut`, `build_algorithm`), the job service
(`POST /v1/gradient`, `/v1/optimize`, `/v1/vqe`, `/v1/qaoa`, `/v1/build`) and
the browser engine (`QuantumEngine.variational(op, args)`). Each answer is
`{"request", "result"}`, the result sealed in a job's receipt and following
`wai.quantum.algorithms/1` ([the schema](../schemas/algorithms.v1.json)), so
`wai-quantum verify` checks it, and runs it again: every call is
deterministic and carries no clock, so the same request gives the same bytes.

## Parameters

A program's parameters are its real-valued inputs, `input float theta;` and
`input angle phi;`. Any job binds them with `"params"`: each value is written
into the program as a `const` of its type at exactly that double, before the
program is read, and every line keeps its number. A parameter left without a
value is refused (`Q114`, `J050`), and so is a name the program does not
declare. A sweep (`engine::sets`) may fix some parameters in `"params"` while
it sweeps the others.

```sh runnable
cat > ansatz.qasm <<'EOF'
OPENQASM 3.0;
input float t;
qubit[2] q;
h q[0];
crx(t) q[0], q[1];
EOF
wai-quantum run ansatz.qasm --param t=0.4 -o Z1       # ⟨Z1⟩ = (1 + cos 0.4)/2
wai-quantum run ansatz.qasm --sweep t=0,0.4,0.8 -o Z1 # one set, each point bound
wai-quantum run ansatz.qasm -o Z1                     # exits 1
```

## Gradients

`gradient` gives `E = ⟨ψ|H|ψ⟩` for a Hamiltonian `H`, a sum of Pauli terms in
the form of a job's observables (`"Z0 Z1"`, `"0.5 X2"`, a constant `"-1.5"`, or
`{"pauli", "coeff"}`), and its derivative by each parameter, on the
double-precision state vector (at most 20 qubits, no noise, measurements only
at the end).

The reader carries every gate's derivative by every parameter
(`qasm::tangents`): through sums, products, quotients, constant powers, the
built-in functions (`sin`, `cos`, `tan`, `exp`, `ln`, `sqrt`, `arcsin`,
`arccos`, `arctan`), loops, broadcast and gate definitions, into each gate's
arguments. What it cannot follow is refused on its line with `Q129`: a
parameter under a modifier (write `crx(θ)` or `rx(-θ)` instead of
`ctrl @ rx(θ)` or `inv @ rx(θ)`), in a classical variable, a condition, an
index, a size, a loop's range or a subroutine, or a type that rounds it
(`angle[n]`, `float[32]`).

Two methods, both exact up to rounding:

- `adjoint` (the default): the circuit run forward once, then back once, each
  gate's derivative applied once against the Hamiltonian applied to the state.
- `parameter_shift`: each gate a parameter moves run again at shifted angles.
  Where its generator has two eigenvalues (`rx`, `ry`, `rz`, `p`, `cp`, `rxx`,
  `ryy`, `rzz`, each angle of `u`, and `cu`'s phases) the rule has two terms,
  `[E(θ + π/2) − E(θ − π/2)]/2`. Where it has three (`crx`, `cry`, `crz`, and
  `cu`'s `θ`: eigenvalues `0, ±½`) it has four, with coefficients
  `(√2 ± 1)/(4√2)` at `±π/2` and `±3π/2`. Every standard gate with a parameter
  has one of the two; a gate with no rule would be refused.

The arithmetic is IEEE-754 double precision, so every number is approximate,
and the answer says so: `precision` gives `"arithmetic": "f64"`,
`"exact": false`, and `rounding_bound`, an a priori bound on the rounding of
the value and of each component. It counts each gate's entries (within four
units of their value) and each amplitude's sum of products, the expectation's
sums over the state in 64 chunks and its terms, the backward pass, and the
reader's own arithmetic for each gate's derivative, to first order in the unit
roundoff `2⁻⁵³`, and doubles the total. It bounds the reported number's
distance from the value the same circuit has in exact arithmetic.

```sh runnable
wai-quantum gradient ansatz.qasm --param t=0.4 --term Z1 --term "0.5 X0"
wai-quantum gradient ansatz.qasm --param t=0.4 --term Z1 --gradient parameter_shift --json > grad.json
wai-quantum verify grad.json
```

The tests hold each method to gradients worked out by hand for every gate
with a parameter within that bound, to central finite differences on a
circuit with every such gate, and to each other.

## Optimisers

`optimize` minimises `E` from the request's `params`, deterministically:

| `optimizer` | step | settings (defaults) |
|---|---|---|
| `gradient_descent` | `θ ← θ − step_size · ∇E` | `step_size` 0.1 |
| `adam` (the default) | the moment estimates with their bias corrections | `step_size` 0.05, `beta1` 0.9, `beta2` 0.999, `epsilon` 1e-8 |
| `spsa` | `ĝ = (E(θ + c_kΔ) − E(θ − c_kΔ))/(2c_k) · Δ`, `θ ← θ − a_k ĝ`, `Δ ∈ {±1}ⁿ` from a seeded stream | `step_size` (`a`) 0.2, `c` 0.1, `A` steps/10, `alpha` 0.602, `gamma` 0.101, `seed` 1 |

with `a_k = a/(k + 1 + A)^α` and `c_k = c/(k + 1)^γ`, at most 2000 `steps`.
The optimisers use only `+ − × ÷` and `√`, the crate's own `exp` and `ln` for
SPSA's gains, and the powers of `β` by repeated product: the same request
gives the same bytes on every platform (pinned in the tests, and held in
wasm32 by the browser engine's).

Every step is listed with its parameters and value (and gradient, or SPSA's
perturbation, estimate and gains). Each step is a link of a digest chain:
`start` is BLAKE3 of `wai.quantum.trace/1\0` and the request's canonical
text, each step's `digest` is BLAKE3 of the previous link and the step's
canonical text without its digest, and `chain` is the last. The receipt seals
the request (the program, the parameters it started from, the method and the
seed) and every step; `wai-quantum verify` checks the chain and runs the
optimisation again. `best` is the lowest value along the path, and `final`
the last point: a minimum along this path, not a proof of the global one.

```sh runnable
wai-quantum optimize ansatz.qasm --param t=0.4 --term Z1 --optimizer adam --steps 50
wai-quantum optimize ansatz.qasm --param t=0.4 --term Z1 --optimizer spsa --steps 50 --seed 7 --json > opt.json
wai-quantum verify opt.json
```

An optimisation can stop between any two steps and go on later, to the same
bytes. Its checkpoint (`wai.quantum.optimize-checkpoint/1`) holds:

- the request it is of, by the chain's `start` and the request's digest;
- the steps taken, each with its `digest`, and the last link;
- the optimiser's whole state: the point, Adam's moments and the powers of its
  betas, SPSA's stream, the best point and the counts. Every number is the
  hex of its f64's bits, so it reads back exactly.

A checkpoint goes on only from its own request, and only if each step is still
the link its digest says. Any other is refused (`J052`).

On the command line, `--checkpoint FILE` keeps the state in FILE as the run goes,
at most once a second and whenever it stops, and goes on from FILE when it is
there. A run killed and started again with the same arguments loses at most a
second of steps, and ends with the same bytes as a run that never stopped. The
file is removed when the run ends. `--slice N` takes at most N steps, then stops
with the state kept (exit 2), so a long optimisation can be run a slice at a
time.

```sh runnable
wai-quantum optimize ansatz.qasm --param t=0.4 --term Z1 --steps 200 --json > whole.json
wai-quantum optimize ansatz.qasm --param t=0.4 --term Z1 --steps 200 --json --checkpoint run.json --slice 100 # exits 2
wai-quantum optimize ansatz.qasm --param t=0.4 --term Z1 --steps 200 --json --checkpoint run.json > resumed.json
cmp whole.json resumed.json
```

## The variational eigensolver

`vqe` optimises an ansatz on a Hamiltonian: the request's program with its
`params`, or the built-in hardware-efficient ansatz (`--layers L`, or
`"ansatz": {"layers": L}`): layers of `ry` then `rz` on every qubit and a ring
of `cx`, the ansatz the variational machine-learning module trains
(`quantum_vml`), its parameters `t0`, `t1`, … starting at 0.1. Up to 8 qubits
the answer sets the best energy beside the exact ground energy by
diagonalisation (`reference.ground_energy`, by the cyclic Jacobi method on the
Hamiltonian's real symmetric form) and states the `gap`. The energy found is
never below the ground energy; an ansatz that cannot reach the ground state
stops above it.

```sh runnable
wai-quantum vqe --term "X0 X1" --term "Y0 Y1" --term "Z0 Z1" --layers 2 --steps 300
```

## QAOA for Max-Cut

`qaoa` builds the program and its cost from a graph (`{"nodes": n, "edges":
[[i, j] | [i, j, w], …]}`, at most 20 vertices): `|+⟩ⁿ`, then for each layer
`l`, `rzz(gamma_l · w)` on every edge and `rx(2 · beta_l)` on every qubit,
with the Hamiltonian `Σ (w/2) Z_i Z_j − (Σ w)/2`, whose expectation is minus
the expected cut. It optimises from the request's `params` or a linear ramp
and reads the best point's distribution against the maximum cut found by
enumerating every assignment: the expected cut and its ratio to the maximum,
the likeliest assignment and whether it is optimal, and the probability of an
optimal one. Assignments are written vertex `n − 1` first, as a job's counts
write a register's bits.

```sh runnable
printf '{"nodes": 4, "edges": [[0, 1], [1, 2], [2, 3], [3, 0]]}\n' > ring.json
wai-quantum qaoa --graph ring.json --layers 2 --steps 100
```

## The algorithm builders

`build` writes an OpenQASM 3 program for a standard algorithm, says how to
read its result and what it gives, and checks it with the equivalence ladder
against a reference built from the algorithm's definition, never from the
program's gates: every basis state at most 8 qubits (`exact`), random product
states at most 10 (`random`, evidence), the adder's truth table at most 16;
wider programs say `unchecked`. Without `--json` the program goes to standard
output and the check to standard error.

| `algorithm` | arguments | the program | read | its reference |
|---|---|---|---|---|
| `qpe` | `counting` (4), `phase` φ in turns | `p(2πφ)` on its eigenstate `\|1⟩`, controlled powers from `t` counting qubits, the inverse Fourier transform | the counting register's integer `k`: `φ ≈ k/2ᵗ` | the Fourier transform by its formula and the powers by their phases |
| `amplitude_estimation` | `counting` (4), `amplitude` a | `A = ry(2·arcsin √a)`, then phase estimation of the Grover operator `Q = A·Z·A†·Z` | `a ≈ sin²(πk/2ᵗ)` | `Q` as `A`'s matrices multiplied, its powers by repeated product |
| `adder` | `bits`, optionally `a`, `b`, `carry_in` | the ripple-carry adder in place: `b ← a + b + c_in mod 2ⁿ`, the carry into `c_out` | with `a` and `b`, the measured sum and carry | the truth table in integer arithmetic |
| `qsvt` | `phases` φ₀…φ_d, and `x` or `matrix` | signal processing `e^{iφ₀Z} Π W(x) e^{iφₖZ}`, `W(x) = e^{i arccos(x) X}`; or the singular value transformation of a real symmetric 2 × 2 `A` through `[[A, √(I−A²)], [√(I−A²), −A]]` | `⟨0\|U\|0⟩ = P(x)`, or the signal qubit's `⟨0\|·\|0⟩` block `P(A)` | the 2 × 2 products, and the block encoding with `√(I − A²)` by its closed form |

Phase estimation's outcome `k` has the Fejér distribution
`P(k) = |2⁻ᵗ Σⱼ e^{2πi j(φ − k/2ᵗ)}|²`, exactly `k = 2ᵗφ` when that is whole;
amplitude estimation's is `½F(k; θ/π) + ½F(k; −θ/π)` with `θ = arcsin √a`.
For `qsvt` the phases are in the signal's convention (`W(x)`); with a matrix
they are turned to the reflection's (`φ₀ + (2d − 1)π/4`, the inner ones
`− π/2`, `φ_d − π/4`), so the block is the same polynomial of `A`. All phases
zero give the Chebyshev polynomial `T_d`. The tests hold phase and amplitude
estimation to those distributions, the adder to a truth table written apart
from the builder's and run as a job on the bit-exact method, and `qsvt` to
`T_d` of a signal and of a matrix by its recursion.

```sh runnable
wai-quantum build qpe --counting 5 --phase 0.3 > qpe.qasm
wai-quantum run qpe.qasm --shots 1000
wai-quantum build adder --bits 4 --a 9 --b 7 > add.qasm
wai-quantum run add.qasm --shots 10
wai-quantum build qsvt --phases 0,0,0,0 --matrix 0.3,0.2,-0.5 --json > qsvt.json
wai-quantum verify qsvt.json
```

## Refusals

| Code | When |
|---|---|
| `J050` | `params` that do not read, a name the program does not declare as a parameter, or a parameter left without a value (with `Q114`) |
| `J051` | a gradient that cannot be taken: a differentiation method that does not exist, a parameter the reader cannot follow (with `Q129`), a program that measures in the middle, noise, no Hamiltonian, or more qubits or work than it takes |
| `J052` | an optimisation, a VQE or a QAOA run that does not read: an unknown optimiser, settings or steps past their limits, a Hamiltonian or a graph that does not read |
| `J053` | a builder's request that does not read: an unknown algorithm, sizes past its limits, or arguments it does not take |
| `J031` | the job service's own bounds (16 qubits, 500 steps, about 4·10⁹ gate amplitudes a call): the message names the command that runs it on your machine |