wai-quantum 0.6.9: cut-rank amplitudes with the circuit's own global phase #350

Merged
dcharlot merged 1 commit from wai/quantum-exact-phase into main 2026-10-20 04:49:01 -04:00
18 changed files with 561 additions and 51 deletions

View file

@ -6,6 +6,36 @@ wire format may change and such changes are recorded here.
## v0.4.0 — current
### wai-quantum 0.6.9: cut-rank amplitudes with the circuit's own global phase
Additive. Every pinned digest is unchanged.
- **`PhasedTableau`** (`quantum_tableau`) is a stabilizer tableau with its global phase.
- It tracks one basis state `r` of the support and its exact amplitude
`⟨r|ψ⟩ = 2^{-k/2} e^{iπm/4}` through every gate.
- Paulis, phase gates, CX, CY, CZ and SWAP move `r` and the phase directly.
- For a Hadamard on qubit `a`, the amplitude of `r ⊕ e_a` comes from the stabilizer
whose X part is `e_a`, found by reducing the stabilizers (about `n³/64` word
operations). `√X` is `H S H`, exactly.
- Tested against the dense state vector after every gate of 40 random circuits
(2 to 8 qubits, over 2,000 gates). A mutant of the Hadamard update is caught.
- **`ProductMeasured::with_last`** keeps the graph form, order and steps with other closing
gates. **`ProductMeasured::scale`** multiplies every amplitude by a factor, folded into
one qubit's weights. Tested against fresh builds.
- **The service's `cutrank` method** returns amplitudes with the circuit's own global
phase.
- The factor is `e^{iΘ} a/b`. Here `a` is the tracked amplitude of the Clifford part,
`b` is the graph form's amplitude at the same basis state (closing gates set to the
identity), and `e^{iΘ}` is the phase the lowering to Clifford gates drops. That
phase is read gate by gate from each gate's matrix against its lowering's.
- Past about a minute of tracking (counted from the Hadamards), the result says
`global_phase` is not fixed, as before.
- `auto` now tries `cutrank` for amplitudes too, before `tensor_network`.
- Tests compare amplitudes with the state vector's directly (no phase aligned). The
circuits include phase-carrying lowerings: `rz`, `rx`, `ry`, `u`, `p`, `cp`, `rzz`,
`rxx`, `iswap`, `sx`, `sxdg`, `cy`, `swap`, `crz` and `cu`. Two mutants are caught:
the lowering phase dropped, and the tracked amplitude ignored.
### wai-quantum 0.6.8: threads only where they pay, and a cut-rank search that recounts only what moved
Additive. No result changes a bit: every kernel path gives the same bits at any thread count,

View file

@ -1017,7 +1017,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
[[package]]
name = "wai-quantum"
version = "0.6.8"
version = "0.6.9"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -599,7 +599,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
[[package]]
name = "wai-quantum"
version = "0.6.8"
version = "0.6.9"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -713,7 +713,7 @@ dependencies = [
[[package]]
name = "wai-quantum"
version = "0.6.8"
version = "0.6.9"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -422,7 +422,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
[[package]]
name = "wai-quantum"
version = "0.6.8"
version = "0.6.9"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -599,7 +599,7 @@ dependencies = [
[[package]]
name = "wai-quantum"
version = "0.6.8"
version = "0.6.9"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -857,7 +857,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
[[package]]
name = "wai-quantum"
version = "0.6.8"
version = "0.6.9"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -41,7 +41,7 @@ The program may use the standard gates, user `gate` definitions, and the control
| `pauli` | any circuit, up to 1024 qubits | expectations | bounded |
| `mps` | any circuit, up to 4096 qubits, at bond dimension up to 256; best when entanglement stays low | counts, expectations, amplitudes | exact when nothing is truncated, otherwise heuristic |
| `tensor_network` | any circuit | amplitudes | exact |
| `cutrank` | a Clifford circuit, then one layer of one-qubit gates at any angles, up to 4096 qubits, as wide as the state's cut-rank allows | counts, amplitudes | exact (amplitudes up to one shared phase); single when only a single-precision table fits |
| `cutrank` | a Clifford circuit, then one layer of one-qubit gates at any angles, up to 4096 qubits, as wide as the state's cut-rank allows | counts, amplitudes | exact; single when only a single-precision table fits |
| `density_matrix` | noisy circuits, up to 16 qubits within a job's memory | counts, expectations | exact |
| `trajectories` | noisy circuits, up to 32 qubits | counts, expectations | estimate |
@ -52,7 +52,7 @@ How the methods run:
- It writes the Clifford part's state as a graph state with local Cliffords, and sums the amplitude one qubit at a time along an order of the qubits, in a table of `2^cut` entries. Here `cut` is the state's cut-rank, its entanglement in bits, across the cut the sum has reached.
- The planner makes one short search for an order with narrower cuts, enough to bound the job's memory and work. The run searches longer, from several seeds on threads of their own, and keeps the cheapest order that is no wider than the planner's. The cost is set by the widest cut, `max_cut` in the result, not by the number of qubits.
- Counts are exact samples, drawn qubit by qubit from the table's conditional probabilities.
- Amplitudes are exact up to one phase shared by all of the job's amplitudes, because a stabilizer tableau does not fix the state's global phase; the result's `global_phase` says so.
- Amplitudes carry the circuit's own global phase. A stabilizer tableau does not fix it, so one amplitude of the Clifford part is tracked exactly, gate by gate, and the summed amplitudes are scaled to agree with it, together with the phase the lowering to Clifford gates drops. A circuit with so many Hadamards that the tracking would take more than about a minute reports `global_phase` as not fixed instead.
- The table is held in double precision when it fits the job's memory, and in single precision otherwise (accuracy class `single`).
- **pauli** reports each value with its truncation `bound`. It holds at most as many Pauli terms as the job's memory allows; a sum that would outgrow it raises its threshold until it fits (`threshold_reached`), and the bound covers everything that dropped. Unless `ladder` is `false`, it also reports a `ladder`: the value and bound at up to three thresholds a decade above the one asked for, so you can see the value settle as the threshold drops.
- **stabilizer**, with noise, samples by Pauli frames against one noiseless outcome of the tableau, 64 shots to a machine word. Its counts have the law of independent noisy shots. Its expectation values are estimates with standard errors, from at least 10,000 shots unless `trajectories` sets the number. Amplitude damping is not a Pauli channel and is refused.
@ -78,7 +78,7 @@ curl -s https://wai.transaction.science/quantum-jobs/v1/run -H 'content-type: ap
"observables": ["Z2", "X2"], "shots": 1000 }'
```
`auto` picks the first method that takes the job. With noise, it tries `density_matrix`, then `stabilizer`, then `trajectories`. For counts, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `cutrank`, then `mps`. For expectations alone, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `pauli`, then `mps`. For amplitudes alone, it tries `exact`, `statevector`, `statevector_single`, `tensor_network`, then `mps`. The job's `why` names the method chosen and the ones that could not take it.
`auto` picks the first method that takes the job. With noise, it tries `density_matrix`, then `stabilizer`, then `trajectories`. For counts, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `cutrank`, then `mps`. For expectations alone, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `pauli`, then `mps`. For amplitudes alone, it tries `exact`, `statevector`, `statevector_single`, `cutrank`, `tensor_network`, then `mps`. The job's `why` names the method chosen and the ones that could not take it.
The accuracy classes:

View file

@ -14,7 +14,7 @@
//! | `pauli` | any, ≤ 1024 qubits | expectations | bounded |
//! | `mps` | any, ≤ 4096 qubits, best on low-entanglement circuits | counts, expectations, amplitudes | exact, or heuristic with its fidelity estimate |
//! | `tensor_network` | any | amplitudes | exact |
//! | `cutrank` | Clifford, then one layer of one-qubit gates at any angles, ≤ 4096 qubits, as wide as its cut-rank allows | counts, amplitudes | exact (amplitudes up to one shared phase), or single |
//! | `cutrank` | Clifford, then one layer of one-qubit gates at any angles, ≤ 4096 qubits, as wide as its cut-rank allows | counts, amplitudes | exact, or single |
//! | `density_matrix` | any, noisy, ≤ 16 qubits in a job's memory | counts, expectations | exact |
//! | `trajectories` | any, noisy, ≤ 32 qubits | counts, expectations | estimate |
//!
@ -22,7 +22,7 @@
//! for counts, `stabilizer`, `exact`, `statevector`, `statevector_single`,
//! `cutrank`, `mps`; for expectations only, `stabilizer`, `exact`, `statevector`,
//! `statevector_single`, `pauli`, `mps`; for amplitudes only, `exact`,
//! `statevector`, `statevector_single`, `tensor_network`, `mps`. A method that cannot take the circuit is
//! `statevector`, `statevector_single`, `cutrank`, `tensor_network`, `mps`. A method that cannot take the circuit is
//! skipped, and the plan says why the one chosen was chosen.
use crate::qasm::{self, Program};
@ -107,7 +107,7 @@ pub fn describe(job_bytes: u64) -> Value {
{ "method": "pauli", "accuracy": "bounded", "circuits": "any", "max_qubits": PAULI_MAX_QUBITS, "answers": ["expectations"], "options": { "threshold": "coefficients below it are dropped (default 1e-5); the bound covers what was dropped", "max_weight": "Pauli strings heavier than this are dropped", "ladder": "also report the value at up to three thresholds a decade apart above it, to show it settling (default true)", "max_terms": "the most Pauli terms to hold (default: what the job's memory holds); past it the threshold rises, and the bound covers what drops" } },
{ "method": "mps", "accuracy": "exact when nothing is truncated, otherwise heuristic", "circuits": "any; best on low-entanglement circuits", "max_qubits": MPS_MAX_QUBITS, "answers": ["counts", "expectations", "amplitudes"], "options": { "max_bond": "the largest bond dimension kept (default 64, at most 256)", "cutoff": "singular values whose squared share is below it are dropped (default 1e-12)" } },
{ "method": "tensor_network", "accuracy": "exact", "circuits": "any", "max_qubits": qasm::MAX_QUBITS, "answers": ["amplitudes"] },
{ "method": "cutrank", "accuracy": "exact (amplitudes up to one phase shared by the job's amplitudes); single when only a single-precision table fits", "circuits": "a Clifford circuit, then one layer of one-qubit gates at any angles: each qubit's gates from its first non-Clifford one on act on that qubit alone; the cost is set by the state's cut-rank along an order of the qubits, not by their number", "max_qubits": CUTRANK_MAX_QUBITS, "answers": ["counts", "amplitudes"] },
{ "method": "cutrank", "accuracy": "exact; single when only a single-precision table fits", "circuits": "a Clifford circuit, then one layer of one-qubit gates at any angles: each qubit's gates from its first non-Clifford one on act on that qubit alone; the cost is set by the state's cut-rank along an order of the qubits, not by their number", "max_qubits": CUTRANK_MAX_QUBITS, "answers": ["counts", "amplitudes"] },
{ "method": "density_matrix", "accuracy": "exact", "circuits": "any, with noise", "max_qubits": density_max_qubits(job_bytes), "answers": ["counts", "expectations"] },
{ "method": "trajectories", "accuracy": "estimate", "circuits": "any, with noise", "max_qubits": TRAJECTORY_MAX_QUBITS.min(statevector_max_qubits(job_bytes)), "answers": ["counts", "expectations"], "options": { "trajectories": "how many (default: the shots, or 1000 for expectations alone); shots that take the same errors share one state" } },
])
@ -291,7 +291,7 @@ pub fn product_measured(r: &Request, thorough: bool) -> Result<wai_quantum::quan
if n > CUTRANK_MAX_QUBITS {
return Err(format!("cutrank holds at most {CUTRANK_MAX_QUBITS} qubits"));
}
let (bulk, last) = r.program.product_measured()?;
let (bulk, last, _) = r.program.product_measured()?;
let t = Tableau::run(n, &bulk).map_err(|e| e.to_string())?;
let pm = ProductMeasured::new(&t, &last).map_err(|e| e.to_string())?;
// Recounting every cut of an order takes about k · cut_k · words
@ -331,6 +331,49 @@ pub fn product_measured(r: &Request, thorough: bool) -> Result<wai_quantum::quan
Ok(best)
}
/// The word operations a `cutrank` job may spend fixing its global phase: an
/// exact amplitude of the Clifford part tracked gate by gate, each Hadamard
/// a reduction of the stabilizers (about `n³/64`). About a minute.
pub const CUTRANK_PHASE_OPS: u64 = 50_000_000_000;
/// The factor that gives `pm`'s amplitudes the circuit's own global phase, or
/// `None` when finding it would cost more than [`CUTRANK_PHASE_OPS`].
///
/// The Clifford part's state, tracked with its phase, has the exact amplitude
/// `a` on a basis state `r`; the graph form `pm` sums has `b` there, with the
/// identity in place of the closing gates. The lowering of the gates to
/// Clifford ones drops a phase `e^{iΘ}` ([`qasm::Program::product_measured`]).
/// The factor is `e^{iΘ} a / b`, of modulus 1.
pub fn cutrank_phase(r: &Request, pm: &wai_quantum::quantum_cutrank::ProductMeasured, cut: usize, precision: wai_quantum::quantum_cutrank::Precision) -> Result<Option<wai_quantum::gates::C>, String> {
use wai_quantum::gates::complex;
use wai_quantum::quantum_tableau::{Clifford, PhasedTableau};
let n = r.program.n_qubits as u64;
let (bulk, _, theta) = r.program.product_measured()?;
let hadamards = bulk.iter().map(|g| match g {
Clifford::H(_) => 1,
Clifford::SX(_) | Clifford::SXdg(_) => 2,
_ => 0,
}).sum::<u64>();
if hadamards.saturating_mul(n * n * n.div_ceil(64)) > CUTRANK_PHASE_OPS {
return Ok(None);
}
let tracked = PhasedTableau::run(n as u32, &bulk).map_err(|e| e.to_string())?;
let (bits, (are, aim)) = tracked.reference();
let (o, one) = (complex(0.0, 0.0), complex(1.0, 0.0));
let bulk_only = pm.with_last(&vec![[one, o, o, one]; n as usize]).map_err(|e| e.to_string())?;
let b = bulk_only.amplitude_in(&bits, cut, precision).map_err(|e| e.to_string())?;
let bb = b.re * b.re + b.im * b.im;
// a / b, then times e^{iΘ}.
let ratio = complex((are * b.re + aim * b.im) / bb, (aim * b.re - are * b.im) / bb);
let k = complex(theta.re * ratio.re - theta.im * ratio.im, theta.re * ratio.im + theta.im * ratio.re);
let off = (k.re.hypot(k.im) - 1.0).abs();
let tolerance = if precision == wai_quantum::quantum_cutrank::Precision::Double { 1e-9 } else { 1e-5 };
if off > tolerance {
return Err(format!("internal: the tracked and summed amplitudes differ in modulus by {off:e}"));
}
Ok(Some(k))
}
/// The memory a `cutrank` job takes: its largest table with a margin, the
/// graph form, and room to run.
pub fn cutrank_bytes(pm: &wai_quantum::quantum_cutrank::ProductMeasured, precision: wai_quantum::quantum_cutrank::Precision) -> u64 {
@ -602,7 +645,7 @@ pub fn plan(v: &Value, job_bytes: u64) -> Result<Plan, String> {
} else if !r.observables.is_empty() && r.amplitudes.is_empty() {
&["stabilizer", "exact", "statevector", "statevector_single", "pauli", "mps"]
} else {
&["exact", "statevector", "statevector_single", "tensor_network", "mps"]
&["exact", "statevector", "statevector_single", "cutrank", "tensor_network", "mps"]
};
let mut refusals = Vec::new();
for m in order {

View file

@ -1244,6 +1244,73 @@ fn one_qubit(kind: Kind, p: &[f64]) -> [C; 4] {
}
}
/// A program split for the `cutrank` method: its Clifford gates, each qubit's
/// closing one-qubit matrix, and the global phase the lowering drops.
pub type ProductMeasuredParts = (Vec<Clifford>, Vec<[C; 4]>, C);
fn cmul(a: C, b: C) -> C {
complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re)
}
/// A Clifford gate as a dense matrix on its qubits (the first the most
/// significant), as the stabilizer tableau reads it.
fn clifford_matrix(c: Clifford) -> Matrix {
let (z, o, i) = (complex(0.0, 0.0), complex(1.0, 0.0), complex(0.0, 1.0));
match c {
Clifford::H(q) => Matrix::h(q),
Clifford::S(q) => Matrix::s(q),
Clifford::Sdg(q) => Matrix::sdg(q),
Clifford::X(q) => Matrix::x(q),
Clifford::Y(q) => Matrix::y(q),
Clifford::Z(q) => Matrix::z(q),
Clifford::SX(q) => Matrix::sqrt_x(q),
Clifford::SXdg(q) => {
let (a, b) = (complex(0.5, -0.5), complex(0.5, 0.5));
Matrix { qubits: vec![q], matrix: vec![a, b, b, a] }
}
Clifford::CX(c, t) => Matrix::cx(c, t),
Clifford::CY(c, t) => Matrix { qubits: vec![c, t], matrix: vec![o, z, z, z, z, o, z, z, z, z, z, complex(0.0, -1.0), z, z, i, z] },
Clifford::CZ(a, b) => Matrix::cz(a, b),
Clifford::Swap(a, b) => Matrix::swap(a, b),
}
}
/// `e^{iθ}` with `op`'s matrix equal to `e^{iθ}` times the product of the
/// Clifford gates `gates` (in order) on the same qubits: the product is built
/// column by column on the op's own qubits, and the phase read off its
/// largest entry.
fn lowering_phase(op: &Op, gates: &[Clifford]) -> C {
let k = op.qubits.len() as u32;
let local = |q: u32| op.qubits.iter().position(|&p| p == q).expect("a lowered gate stays on its op's qubits") as u32;
let relabel = |c: Clifford| -> Matrix {
let mut m = clifford_matrix(c);
m.qubits = m.qubits.iter().map(|&q| local(q)).collect();
m
};
let local_gates: Vec<Matrix> = gates.iter().map(|&c| relabel(c)).collect();
let d = 1usize << k;
let want = op.matrix().matrix;
let mut best = (0.0f64, complex(1.0, 0.0));
for col in 0..d {
let mut sv = wai_quantum::quantum_sv::StateVector::zero(k).expect("at most three qubits");
sv.re[0] = 0.0;
sv.re[col] = 1.0;
for g in &local_gates {
sv.apply(g).expect("a lowered gate on the op's qubits");
}
for row in 0..d {
let (lr, li) = (sv.re[row], sv.im[row]);
let mag = lr * lr + li * li;
if mag > best.0 + 1e-12 {
// want / lowered
let w = want[row * d + col];
best = (mag, complex((w.re * lr + w.im * li) / mag, (w.im * lr - w.re * li) / mag));
}
}
}
best.1
}
impl Op {
/// The gate as a matrix on its qubits (the first the most significant).
pub fn matrix(&self) -> Matrix {
@ -1737,7 +1804,11 @@ impl Program {
/// into one matrix (the identity for a qubit with none). The input of the
/// `cutrank` method. Refused when a gate on two or more qubits comes after
/// a qubit's first non-Clifford gate, or is not Clifford itself.
pub fn product_measured(&self) -> Result<(Vec<Clifford>, Vec<[C; 4]>), String> {
///
/// The third member is the global phase the Clifford lowering drops: each
/// gate's matrix is `e^{iθ}` times the product of the Clifford gates it is
/// written as, and this is the product of the `e^{iθ}`.
pub fn product_measured(&self) -> Result<ProductMeasuredParts, String> {
let one = complex(1.0, 0.0);
let zero = complex(0.0, 0.0);
let mul = |a: &[C], b: &[C; 4]| -> [C; 4] {
@ -1746,6 +1817,7 @@ impl Program {
[add(m(a[0], b[0]), m(a[1], b[2])), add(m(a[0], b[1]), m(a[1], b[3])), add(m(a[2], b[0]), m(a[3], b[2])), add(m(a[2], b[1]), m(a[3], b[3]))]
};
let mut bulk = Vec::new();
let mut phase = one;
let mut tail: Vec<Option<[C; 4]>> = vec![None; self.n_qubits as usize];
for op in &self.ops {
if let [q] = op.qubits[..] {
@ -1754,7 +1826,10 @@ impl Program {
tail[q] = Some(mul(&op.matrix().matrix, &t));
} else {
match Self::lower_clifford(std::slice::from_ref(op)) {
Ok((gates, _)) => bulk.extend(gates),
Ok((gates, _)) => {
phase = cmul(phase, lowering_phase(op, &gates));
bulk.extend(gates);
}
Err(_) => tail[q] = Some(mul(&op.matrix().matrix, &[one, zero, zero, one])),
}
}
@ -1764,9 +1839,12 @@ impl Program {
return Err(format!("{:?} on qubits {:?} comes after a gate on qubit {q} that is not Clifford; cutrank takes a Clifford circuit, then one layer of one-qubit gates", op.kind, op.qubits));
}
let (gates, _) = Self::lower_clifford(std::slice::from_ref(op)).map_err(|e| format!("{e}; cutrank takes a Clifford circuit, then one layer of one-qubit gates"))?;
phase = cmul(phase, lowering_phase(op, &gates));
bulk.extend(gates);
}
Ok((bulk, tail.into_iter().map(|t| t.unwrap_or([one, zero, zero, one])).collect()))
// Held to modulus 1 against the rounding of many products.
let m = phase.re.hypot(phase.im);
Ok((bulk, tail.into_iter().map(|t| t.unwrap_or([one, zero, zero, one])).collect(), complex(phase.re / m, phase.im / m)))
}
/// The Clifford lowering, and where each op's gates end in it.

View file

@ -372,7 +372,7 @@ pub fn run(request: &Value, plan: &Plan) -> Result<Value, String> {
"class": if precision == Precision::Double { "exact" } else { "single" },
"max_cut": cut,
"order_work": pm.work(),
"note": "the state's sum taken qubit by qubit along an order of the qubits, in a table of 2^max_cut entries (single class: the table stored in single precision, its arithmetic in double); counts are exact samples, seeded; amplitudes are exact up to one phase shared by all of the job's amplitudes, since a stabilizer tableau does not fix the state's global phase"
"note": "the state's sum taken qubit by qubit along an order of the qubits, in a table of 2^max_cut entries (single class: the table stored in single precision, its arithmetic in double); counts are exact samples, seeded; amplitudes carry the circuit's own global phase, found by tracking one amplitude of the Clifford part exactly through the circuit"
});
if r.shots > 0 {
let mut map: HashMap<String, u64> = HashMap::new();
@ -382,6 +382,14 @@ pub fn run(request: &Value, plan: &Plan) -> Result<Value, String> {
merge(&mut out, counts(map));
}
if !r.amplitudes.is_empty() {
// The circuit's global phase: the lowering's own, times the
// ratio of one exactly tracked amplitude of the Clifford part
// to the graph form's at the same basis state.
let mut pm = pm;
let fixed = catalog::cutrank_phase(&r, &pm, cut, precision)?;
if let Some(kappa) = fixed {
pm.scale(kappa);
}
let mut entries = Vec::new();
for b in &r.amplitudes {
let x: Vec<bool> = (0..n).map(|q| b.as_bytes()[(n - 1 - q) as usize] == b'1').collect();
@ -389,7 +397,9 @@ pub fn run(request: &Value, plan: &Plan) -> Result<Value, String> {
entries.push(json!({ "bits": b, "re": a.re, "im": a.im }));
}
out["amplitudes"] = Value::Array(entries);
out["global_phase"] = json!("not fixed: the amplitudes share one unknown phase; their magnitudes and relative phases are exact");
if fixed.is_none() {
out["global_phase"] = json!("not fixed: tracking it exactly would take more than one job may; the amplitudes share one unknown phase, and their magnitudes and relative phases are exact");
}
}
}
"tensor_network" => {
@ -943,35 +953,34 @@ mod tests {
q
}
#[test]
fn cutrank_amplitudes_are_the_state_vectors_up_to_one_phase() {
for (n, seed) in [(6u32, 1u64), (9, 2), (11, 3)] {
let qasm = clifford_then_layer(n, 3 * n as usize, seed);
let bits: Vec<String> = (0..24u64).map(|k| (0..n).map(|q| if (k.wrapping_mul(0x9e37_79b9) >> q) & 1 == 1 { '1' } else { '0' }).collect()).collect();
let cr = go(json!({ "qasm": qasm, "method": "cutrank", "amplitudes": bits }));
let sv = go(json!({ "qasm": qasm, "method": "statevector", "amplitudes": bits }));
assert_eq!(cr["method"], "cutrank", "{cr}");
assert_eq!(cr["accuracy"]["class"], "exact", "{cr}");
assert!(cr["global_phase"].as_str().is_some_and(|g| g.starts_with("not fixed")), "{cr}");
let amp = |v: &Value, i: usize| (v["amplitudes"][i]["re"].as_f64().unwrap(), v["amplitudes"][i]["im"].as_f64().unwrap());
// The shared phase, from the largest amplitude.
let k = (0..bits.len()).max_by(|&a, &b| {
let (x, y) = (amp(&sv, a), amp(&sv, b));
(x.0 * x.0 + x.1 * x.1).total_cmp(&(y.0 * y.0 + y.1 * y.1))
}).unwrap();
let ((cr_k_re, cr_k_im), (sv_k_re, sv_k_im)) = (amp(&cr, k), amp(&sv, k));
let norm = sv_k_re * sv_k_re + sv_k_im * sv_k_im;
// phase = cr_k / sv_k
let (pr, pi) = ((cr_k_re * sv_k_re + cr_k_im * sv_k_im) / norm, (cr_k_im * sv_k_re - cr_k_re * sv_k_im) / norm);
assert!(((pr * pr + pi * pi) - 1.0).abs() < 1e-9, "n {n}: |phase| {}", (pr * pr + pi * pi).sqrt());
for (i, b) in bits.iter().enumerate() {
let ((cr_re, cr_im), (sv_re, sv_im)) = (amp(&cr, i), amp(&sv, i));
let (want_re, want_im) = (pr * sv_re - pi * sv_im, pr * sv_im + pi * sv_re);
assert!((cr_re - want_re).abs() < 1e-10 && (cr_im - want_im).abs() < 1e-10, "n {n}, {b}: cutrank {cr_re} {cr_im}, statevector × phase {want_re} {want_im}");
}
/// cutrank's amplitudes, global phase and all, are the state vector's.
fn assert_amplitudes_match(qasm: &str, n: u32, label: &str) {
let bits: Vec<String> = (0..24u64).map(|k| (0..n).map(|q| if (k.wrapping_mul(0x9e37_79b9) >> q) & 1 == 1 { '1' } else { '0' }).collect()).collect();
let cr = go(json!({ "qasm": qasm, "method": "cutrank", "amplitudes": bits }));
let sv = go(json!({ "qasm": qasm, "method": "statevector", "amplitudes": bits }));
assert_eq!(cr["method"], "cutrank", "{label}: {cr}");
assert_eq!(cr["accuracy"]["class"], "exact", "{label}: {cr}");
assert!(cr.get("global_phase").is_none(), "{label}: the phase is fixed: {cr}");
for (i, b) in bits.iter().enumerate() {
let get = |v: &Value| (v["amplitudes"][i]["re"].as_f64().unwrap(), v["amplitudes"][i]["im"].as_f64().unwrap());
let ((cre, cim), (sre, sim)) = (get(&cr), get(&sv));
assert!((cre - sre).abs() < 1e-10 && (cim - sim).abs() < 1e-10, "{label}, {b}: cutrank {cre} {cim}, statevector {sre} {sim}");
}
}
#[test]
fn cutrank_amplitudes_are_the_state_vectors_phase_and_all() {
for (n, seed) in [(6u32, 1u64), (9, 2), (11, 3)] {
assert_amplitudes_match(&clifford_then_layer(n, 3 * n as usize, seed), n, &format!("n {n}"));
}
// Gates whose lowering to Clifford ones drops a phase, in the Clifford part.
let phased = "OPENQASM 3.0; qubit[5] q; h q; rz(pi/2) q[0]; rx(pi) q[1]; ry(pi/2) q[2]; u(pi/2, 0, pi) q[3]; p(pi/2) q[4]; \
cp(pi) q[0], q[1]; rzz(pi/2) q[1], q[2]; iswap q[2], q[3]; sx q[4]; sxdg q[0]; cy q[3], q[4]; swap q[0], q[2]; \
rxx(pi/2) q[3], q[1]; crz(pi) q[4], q[0]; cu(pi, 0, pi, pi/2) q[2], q[1]; \
t q[0]; ry(0.4) q[1]; rz(1.3) q[2]; h q[3]; rx(0.7) q[3]; tdg q[4];";
assert_amplitudes_match(phased, 5, "phase-carrying lowerings");
}
#[test]
fn cutrank_counts_follow_the_exact_distribution() {
let n = 5u32;

View file

@ -1079,7 +1079,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
[[package]]
name = "wai-quantum"
version = "0.6.8"
version = "0.6.9"
dependencies = [
"blake3",
"block2",

View file

@ -1,6 +1,6 @@
[package]
name = "wai-quantum"
version = "0.6.8"
version = "0.6.9"
edition = "2024"
rust-version = "1.97"
license = "Apache-2.0"

View file

@ -467,6 +467,16 @@ the search takes one qubit to another place, and only the cuts between the two
places are recounted, so 300 000 moves on the 70-qubit circuit below take a few
seconds. `sample` prints each exact sample with `2^n · p`.
A stabilizer tableau fixes a state only up to a global phase, and the graph
form inherits that. `PhasedTableau` carries the phase: alongside the tableau it
tracks one basis state of the support and its exact amplitude, a power of
`2^{-1/2}` times an eighth root of unity, gate by gate. A Hadamard reads the
stabilizer whose X part is that one qubit, which is a reduction of the
stabilizers. `ProductMeasured::scale` then brings every amplitude to the
circuit's own phase: by the ratio of the tracked amplitude to the graph form's
at the same basis state (`with_last` with identities), times whatever phase
the caller's lowering to Clifford gates dropped.
A 70-qubit graph-state sampling circuit has a nearest-neighbour Clifford bulk,
then a non-Clifford basis on every qubit. With an annealed order, every cut is
29 bits or fewer. On one workstation, one amplitude takes about half a second

View file

@ -67,6 +67,12 @@ fn mul(a: &Unitary, b: &Unitary) -> Unitary {
[add(m(a[0], b[0]), m(a[1], b[2])), add(m(a[0], b[1]), m(a[1], b[3])), add(m(a[2], b[0]), m(a[3], b[2])), add(m(a[2], b[1]), m(a[3], b[3]))]
}
/// `W_q / √2 = V_q L_q / √2` for each qubit.
fn weights(local: &[Unitary], last: &[Unitary]) -> Vec<Unitary> {
let r = std::f64::consts::FRAC_1_SQRT_2;
local.iter().zip(last).map(|(l, v)| mul(v, l).map(|c| complex(c.re * r, c.im * r))).collect()
}
fn get(v: &[u64], i: usize) -> bool {
v[i / 64] >> (i % 64) & 1 == 1
}
@ -78,6 +84,8 @@ pub struct ProductMeasured {
adjacency: Vec<Vec<u64>>,
/// `W_q / √2`.
w: Vec<Unitary>,
/// The graph form's local Cliffords `L_q`, so the closing gates can change.
local: Vec<Unitary>,
order: Vec<usize>,
/// The cut-rank after each qubit of the order.
cuts: Vec<usize>,
@ -105,7 +113,7 @@ impl ProductMeasured {
let h = [complex(r, 0.0), complex(r, 0.0), complex(r, 0.0), complex(-r, 0.0)];
let s = [one, o, o, complex(0.0, 1.0)];
let z = [one, o, o, complex(-1.0, 0.0)];
let w = (0..n)
let local: Vec<Unitary> = (0..n)
.map(|q| {
let l = form.local[q];
let mut m = [one, o, o, one];
@ -118,16 +126,37 @@ impl ProductMeasured {
if l.z {
m = mul(&m, &z);
}
let wq = mul(&last[q], &m);
wq.map(|c| complex(c.re * r, c.im * r))
m
})
.collect();
let mut pm = ProductMeasured { n, adjacency: form.adjacency.clone(), w, order: (0..n).collect(), cuts: Vec::new(), steps: Vec::new() };
let w = weights(&local, last);
let mut pm = ProductMeasured { n, adjacency: form.adjacency.clone(), w, local, order: (0..n).collect(), cuts: Vec::new(), steps: Vec::new() };
pm.cuts = pm.cut_ranks_of(&pm.order);
pm.steps = plan(&pm.adjacency, n, &pm.order);
pm
}
/// The same graph form, order and steps with other closing gates:
/// `last[q]` on qubit q.
pub fn with_last(&self, last: &[Unitary]) -> Result<ProductMeasured, CutRankError> {
if last.len() != self.n {
return Err(CutRankError::Gates { qubits: self.n, gates: last.len() });
}
let mut pm = self.clone();
pm.w = weights(&self.local, last);
Ok(pm)
}
/// Every amplitude times `factor`. Each term of the sum holds exactly one
/// entry of each qubit's weights, so it is folded into one qubit's; a
/// factor of modulus 1 (a phase) leaves every probability as it was.
pub fn scale(&mut self, factor: C) {
let q = self.order.first().copied().unwrap_or(0);
if let Some(wq) = self.w.get_mut(q) {
*wq = wq.map(|c| complex(c.re * factor.re - c.im * factor.im, c.re * factor.im + c.im * factor.re));
}
}
/// The same state, summed in `order`.
pub fn with_order(mut self, order: &[usize]) -> Result<ProductMeasured, CutRankError> {
let mut seen = vec![false; self.n];
@ -1164,6 +1193,34 @@ mod tests {
assert!((a.re - b.re).hypot(a.im - b.im) < 1e-9 * a.re.hypot(a.im).max(1e-300), "{a:?} vs {b:?}");
}
/// Other closing gates on the same graph form and order give the same
/// amplitudes as a fresh build; a phase scales every amplitude by itself.
#[test]
fn closing_gates_change_and_a_phase_scales() {
let n = 10usize;
let mut s = 11u64;
let ops = random_clifford(n as u32, 6 * n, 23);
let t = Tableau::run(n as u32, &ops).unwrap();
let first: Vec<Unitary> = (0..n).map(|_| random_unitary(&mut s)).collect();
let second: Vec<Unitary> = (0..n).map(|_| random_unitary(&mut s)).collect();
let pm = ProductMeasured::new(&t, &first).unwrap();
let order = pm.anneal_order(200, 1);
let pm = pm.with_order(&order).unwrap();
let swapped = pm.with_last(&second).unwrap();
let fresh = ProductMeasured::new(&t, &second).unwrap().with_order(&order).unwrap();
let phase = complex(0.6, 0.8);
let mut scaled = fresh.clone();
scaled.scale(phase);
for k in 0..20 {
let x: Vec<bool> = (0..n).map(|q| (k * 7 + q * 3) % 5 < 2).collect();
let (a, b, c) = (swapped.amplitude(&x, n).unwrap(), fresh.amplitude(&x, n).unwrap(), scaled.amplitude(&x, n).unwrap());
assert!((a.re - b.re).abs() < 1e-14 && (a.im - b.im).abs() < 1e-14, "{a:?} vs {b:?}");
let want = complex(b.re * phase.re - b.im * phase.im, b.re * phase.im + b.im * phase.re);
assert!((c.re - want.re).abs() < 1e-14 && (c.im - want.im).abs() < 1e-14, "{c:?} vs {want:?}");
}
assert!(matches!(pm.with_last(&second[..3]), Err(CutRankError::Gates { .. })));
}
#[test]
fn cut_ranks_are_the_entanglement() {
// A GHZ chain: one bit across every cut. A product state: none.

View file

@ -691,6 +691,211 @@ impl Tableau {
}
}
impl Tableau {
/// The stabilizer whose X part is exactly `e_a`, if the group holds one:
/// its Z part and its sign, as a row of letters (`Y` where both are set).
/// Found by reducing the stabilizers' X parts: in reduced form a vector
/// lies in their span only through the row whose pivot it sets.
fn stabilizer_with_x(&self, a: usize) -> Option<(Vec<u64>, bool)> {
let n = self.n;
let mut rows = self.rows(n, 2 * n);
let mut k = 0;
let mut pivot_row = None;
for q in 0..n {
let Some(p) = (k..n).find(|&i| get(rows.xr(i), q)) else { continue };
rows.swap(k, p);
for i in 0..n {
if i != k && get(rows.xr(i), q) {
rows.mul_into(i, k);
}
}
if q == a {
pivot_row = Some(k);
}
k += 1;
}
let j = pivot_row?;
let x = rows.xr(j);
let only_a = x.iter().enumerate().all(|(w, &bits)| bits == if w == a / 64 { 1u64 << (a % 64) } else { 0 });
only_a.then(|| (rows.zr(j).to_vec(), rows.r[j]))
}
}
/// A stabilizer state with its global phase: the [`Tableau`], which fixes the
/// state only up to a phase, and one basis state `r` of its support with the
/// exact amplitude `⟨r|ψ⟩ = 2^{−k/2} e^{iπm/4}`. Every gate but `H` moves `r`
/// and the phase directly. For `H` on qubit `a` the amplitude of `r ⊕ e_a` is
/// needed: it is `c · ⟨r|ψ⟩` with `c` read from the stabilizer whose X part is
/// `e_a` (`c = 0` when there is none), which costs a reduction of the
/// stabilizers, about `n³/64` word operations.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct PhasedTableau {
t: Tableau,
r: Vec<u64>,
/// `|⟨r|ψ⟩| = 2^{−k/2}`.
k: u32,
/// The phase, in eighths of a turn.
m: u8,
}
impl PhasedTableau {
/// `|0…0⟩`, with amplitude 1 on `|0…0⟩`.
pub fn new(n: u32) -> Result<PhasedTableau, TableauError> {
Ok(PhasedTableau { t: Tableau::new(n)?, r: vec![0; words(n as usize)], k: 0, m: 0 })
}
/// The circuit `ops` on `n` qubits from `|0…0⟩`.
pub fn run(n: u32, ops: &[Clifford]) -> Result<PhasedTableau, TableauError> {
let mut t = PhasedTableau::new(n)?;
for &op in ops {
t.apply(op)?;
}
Ok(t)
}
/// The tableau, without the phase.
pub fn tableau(&self) -> &Tableau {
&self.t
}
/// The reference basis state (qubit `q` at index `q`) and `⟨r|ψ⟩`.
pub fn reference(&self) -> (Vec<bool>, (f64, f64)) {
let bits = (0..self.t.n).map(|q| get(&self.r, q)).collect();
let mag = (-(self.k as f64) / 2.0).exp2();
// Exact for the eight angles, so the entries are the same bits everywhere.
let (c, s) = match self.m % 8 {
0 => (1.0, 0.0),
2 => (0.0, 1.0),
4 => (-1.0, 0.0),
6 => (0.0, -1.0),
1 => (std::f64::consts::FRAC_1_SQRT_2, std::f64::consts::FRAC_1_SQRT_2),
3 => (-std::f64::consts::FRAC_1_SQRT_2, std::f64::consts::FRAC_1_SQRT_2),
5 => (-std::f64::consts::FRAC_1_SQRT_2, -std::f64::consts::FRAC_1_SQRT_2),
_ => (std::f64::consts::FRAC_1_SQRT_2, -std::f64::consts::FRAC_1_SQRT_2),
};
(bits, (mag * c, mag * s))
}
fn bit(&self, q: usize) -> bool {
get(&self.r, q)
}
fn flip(&mut self, q: usize) {
self.r[q / 64] ^= 1 << (q % 64);
}
fn turn(&mut self, eighths: u8) {
self.m = (self.m + eighths) % 8;
}
/// `H` on qubit `a`, the reference first (it needs the state before).
fn h(&mut self, a: usize) {
// c = ⟨r ⊕ e_a|ψ⟩ / ⟨r|ψ⟩ in quarter turns, or None when r ⊕ e_a is
// outside the support. With P = (−1)^s i^y X_a Z^z a stabilizer,
// P|r⟩ = (−1)^{s + z·r} i^y |r ⊕ e_a⟩.
let c = self.t.stabilizer_with_x(a).map(|(z, s)| {
let y = get(&z, a);
let parity = z.iter().zip(&self.r).map(|(a, b)| (a & b).count_ones()).sum::<u32>() % 2 == 1;
(u8::from(s ^ parity) * 2 + u8::from(y)) % 4
});
let b = self.bit(a);
// ⟨r|Hψ⟩ / ⟨r|ψ⟩ is (1 + c)/√2 for b = 0 and (c − 1)/√2 for b = 1; when
// that is zero, r ⊕ e_a takes over with (c + 1)/√2 or (1 − c)/√2.
match (c, b) {
(None, false) => self.k += 1,
(None, true) => {
self.k += 1;
self.turn(4);
}
(Some(0), false) => self.k -= 1,
(Some(0), true) => {
self.k -= 1;
self.flip(a);
}
(Some(2), false) => {
self.k -= 1;
self.flip(a);
}
(Some(2), true) => {
self.k -= 1;
self.turn(4);
}
(Some(1), false) => self.turn(1),
(Some(1), true) => self.turn(3),
(Some(_), false) => self.turn(7),
(Some(_), true) => self.turn(5),
}
self.t.h(a);
}
/// Applies one gate.
pub fn apply(&mut self, op: Clifford) -> Result<(), TableauError> {
match op {
Clifford::H(a) => {
let a = self.t.check(a)?;
self.h(a);
return Ok(());
}
Clifford::SX(a) | Clifford::SXdg(a) => {
// √X = H S H and √X† = H S† H, exactly.
let a = self.t.check(a)?;
self.h(a);
self.apply(if matches!(op, Clifford::SX(_)) { Clifford::S(a as u32) } else { Clifford::Sdg(a as u32) })?;
self.h(a);
return Ok(());
}
Clifford::X(a) => {
let a = self.t.check(a)?;
self.flip(a);
}
Clifford::Y(a) => {
// Y|b⟩ = i(−1)^b |1 − b⟩.
let a = self.t.check(a)?;
self.turn(2 + 4 * u8::from(self.bit(a)));
self.flip(a);
}
Clifford::Z(a) => {
let a = self.t.check(a)?;
self.turn(4 * u8::from(self.bit(a)));
}
Clifford::S(a) => {
let a = self.t.check(a)?;
self.turn(2 * u8::from(self.bit(a)));
}
Clifford::Sdg(a) => {
let a = self.t.check(a)?;
self.turn(6 * u8::from(self.bit(a)));
}
Clifford::CX(a, b) => {
let (a, b) = self.t.pair(a, b)?;
if self.bit(a) {
self.flip(b);
}
}
Clifford::CY(a, b) => {
let (a, b) = self.t.pair(a, b)?;
if self.bit(a) {
self.turn(2 + 4 * u8::from(self.bit(b)));
self.flip(b);
}
}
Clifford::CZ(a, b) => {
let (a, b) = self.t.pair(a, b)?;
self.turn(4 * u8::from(self.bit(a) && self.bit(b)));
}
Clifford::Swap(a, b) => {
let (a, b) = self.t.pair(a, b)?;
let (x, y) = (self.bit(a), self.bit(b));
if x != y {
self.flip(a);
self.flip(b);
}
}
}
self.t.apply(op)
}
}
/// The computational-basis outcomes of a stabilizer state: uniform over
/// `x₀ + span(basis)`. Bit `q` of an outcome is qubit `q`.
#[derive(Clone, Debug, PartialEq, Eq)]
@ -1593,3 +1798,81 @@ mod frame_tests {
assert!(e[0].0 > 0.8 && e[0].0 < 1.0, "{e:?}");
}
}
#[cfg(all(test, feature = "quantum_sv"))]
mod phase_tests {
use super::*;
use crate::gates::{complex, Gate};
use crate::quantum_sv::StateVector;
fn gate(c: Clifford) -> Gate {
let (z, o, i) = (complex(0.0, 0.0), complex(1.0, 0.0), complex(0.0, 1.0));
match c {
Clifford::H(q) => Gate::h(q),
Clifford::S(q) => Gate::s(q),
Clifford::Sdg(q) => Gate::sdg(q),
Clifford::X(q) => Gate::x(q),
Clifford::Y(q) => Gate::y(q),
Clifford::Z(q) => Gate::z(q),
Clifford::SX(q) => Gate::sqrt_x(q),
Clifford::SXdg(q) => {
let (a, b) = (complex(0.5, -0.5), complex(0.5, 0.5));
Gate { qubits: vec![q], matrix: vec![a, b, b, a] }
}
Clifford::CX(c, t) => Gate::cx(c, t),
Clifford::CY(c, t) => Gate { qubits: vec![c, t], matrix: vec![o, z, z, z, z, o, z, z, z, z, z, complex(0.0, -1.0), z, z, i, z] },
Clifford::CZ(a, b) => Gate::cz(a, b),
Clifford::Swap(a, b) => Gate::swap(a, b),
}
}
fn circuit(n: u32, len: usize, seed: u64) -> Vec<Clifford> {
let mut s = seed;
(0..len)
.map(|_| {
let a = (splitmix(&mut s) % u64::from(n)) as u32;
let b = (a + 1 + (splitmix(&mut s) % u64::from(n - 1)) as u32) % n;
// Hadamards often, so every case of the reference's update is met.
match splitmix(&mut s) % 14 {
0..=2 => Clifford::H(a),
3 => Clifford::S(a),
4 => Clifford::Sdg(a),
5 => Clifford::SX(a),
6 => Clifford::SXdg(a),
7 => Clifford::X(a),
8 => Clifford::Y(a),
9 => Clifford::Z(a),
10 => Clifford::CX(a, b),
11 => Clifford::CY(a, b),
12 => Clifford::CZ(a, b),
_ => Clifford::Swap(a, b),
}
})
.collect()
}
/// The reference amplitude is the state vector's, phase and all, after
/// every gate of random Clifford circuits.
#[test]
fn the_reference_amplitude_is_the_state_vectors() {
let mut checked = 0;
for seed in 0..40u64 {
let n = 2 + (seed % 7) as u32;
let ops = circuit(n, 12 * n as usize, seed);
let mut t = PhasedTableau::new(n).unwrap();
let mut sv = StateVector::zero(n).unwrap();
for (g, &op) in ops.iter().enumerate() {
t.apply(op).unwrap();
sv.apply(&gate(op)).unwrap();
let (bits, (re, im)) = t.reference();
let index = bits.iter().enumerate().fold(0usize, |x, (q, &b)| x | (usize::from(b) << (n as usize - 1 - q)));
let a = sv.amplitude(index);
assert!((a.re - re).abs() < 1e-12 && (a.im - im).abs() < 1e-12, "seed {seed}, gate {g} ({op:?}): tableau {re} {im}, state vector {} {}", a.re, a.im);
assert!(re.hypot(im) > 1e-9, "the reference stays in the support");
checked += 1;
}
assert_eq!(t.tableau(), &Tableau::run(n, &ops).unwrap(), "the tableau is the plain one");
}
assert!(checked > 2000);
}
}

2
wai/wai-rs/Cargo.lock generated
View file

@ -4310,7 +4310,7 @@ version = "0.1.0"
[[package]]
name = "wai-quantum"
version = "0.6.8"
version = "0.6.9"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -1241,7 +1241,7 @@ version = "0.1.0"
[[package]]
name = "wai-quantum"
version = "0.6.8"
version = "0.6.9"
dependencies = [
"blake3",
"ed25519-dalek",