From 9bb3da572e57376aa03428b3a4366552a2ad0698 Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Thu, 1 Oct 2026 20:40:14 -0400 Subject: [PATCH] wai-quantum 0.3.23: belief-propagation tensor networks, a second method for the same circuits MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit quantum_bptn holds a tensor-network state on the hardware's coupling graph: one tensor per qubit, one bond per edge, in Vidal gauge. - A two-qubit gate on an edge is a simple update: QR-reduce both sides, apply the gate, SVD, keep at most χ singular values. - Expectations take each tensor's environment from its bond weights, or from belief propagation run to a fixed point on the norm network. regauge() brings the network into the belief-propagation gauge, and run_gauged() does so after every layer. - Its errors are unrelated to those of sparse Pauli dynamics, so the two can check each other on the same circuit. Measured: - On trees with no truncation it is exact: 3.7e-15 against exact sparse Pauli dynamics, and the same with belief-propagation environments. - 127-qubit heavy-hex, 5 steps, χ = 32: the exactly computed magnetization to within 3.5e-7 at θ_h = kπ/32 for k = 2, 4, 6, 8. - 127-qubit heavy-hex, 20 steps, ⟨Z62⟩: 0.9786 at θ_h = 0.2 and 0.8498 at θ_h = 0.4, the published belief-propagation values (arXiv:2306.14887) at the same angles. - At θ_h = 0.6 it approaches the published value as χ grows (0.696 at χ = 16, 0.676 at χ = 32, against 0.669). edge_layers() partitions a graph's edges into layers of vertex-disjoint edges: Δ layers on a bipartite graph (König), three on heavy-hex. kicked_ising() now applies its RZZ gates a layer at a time, as the device does. For sparse Pauli dynamics this changes no bit, because the gates are commuting exact Cliffords. A layer's gates touch disjoint tensors and bonds, so the tensor network computes them in parallel and commits them in order, bit-identical to applying them one at a time. src/linalg.rs holds the complex type, the Jacobi SVD and a Householder QR, shared with the MPS backend, whose arithmetic is unchanged. The SVD the new backend uses runs to convergence. The original's stopping rule leaves the singular vectors orthonormal only to about 1e-7. That is harmless where a backend keeps U·S·V† whole, as MPS does, and is wrong to 1e-8 where expectations assume orthonormal tensors. Its |γ| is computed from / and √: libm's hypot made the native and WASI results differ in the last bit until it went. tests/bit_identity.rs pins a tensor-network run, confirmed natively and under wasm32-wasip2. tests/utility_reference.rs holds the 5-step values by both methods. `wai-quantum utility --method bptn --chi N` runs it from the CLI, with output byte-identical natively and as a WASI component. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum-serve/Cargo.lock | 2 +- wai/wai-quantum/Cargo.lock | 2 +- wai/wai-quantum/Cargo.toml | 13 +- wai/wai-quantum/README.md | 5 +- wai/wai-quantum/examples/utility_bptn.rs | 52 + wai/wai-quantum/src/lib.rs | 6 + wai/wai-quantum/src/linalg.rs | 416 ++++++++ wai/wai-quantum/src/quantum_bptn.rs | 1100 ++++++++++++++++++++ wai/wai-quantum/src/quantum_mps.rs | 136 +-- wai/wai-quantum/src/quantum_spd.rs | 114 +- wai/wai-quantum/tests/bit_identity.rs | 16 +- wai/wai-quantum/tests/utility_reference.rs | 16 + wai/wai-rs/Cargo.lock | 2 +- wai/wai-rs/Cargo.toml | 3 + wai/wai-rs/src/bin/wai_quantum.rs | 33 +- wai/wai-rs/src/lib.rs | 3 + wai/wai-web/Cargo.lock | 2 +- 17 files changed, 1773 insertions(+), 148 deletions(-) create mode 100644 wai/wai-quantum/examples/utility_bptn.rs create mode 100644 wai/wai-quantum/src/linalg.rs create mode 100644 wai/wai-quantum/src/quantum_bptn.rs diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index 5d62fd3c3..5bdd91421 100644 --- a/wai/wai-quantum-serve/Cargo.lock +++ b/wai/wai-quantum-serve/Cargo.lock @@ -589,7 +589,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.22" +version = "0.3.23" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index dfcb512b1..8016907e4 100644 --- a/wai/wai-quantum/Cargo.lock +++ b/wai/wai-quantum/Cargo.lock @@ -415,7 +415,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.3.22" +version = "0.3.23" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index 783681c9a..17792bacf 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -1,10 +1,10 @@ [package] name = "wai-quantum" -version = "0.3.22" +version = "0.3.23" edition = "2024" rust-version = "1.97" license = "Apache-2.0" -description = "A deterministic quantum stack in pure Rust: byte-exact circuit simulation (statevector / stabilizer / tensor-network MPS / sparse-Pauli backends), sparse Pauli dynamics at utility scale (arbitrary angles, 1024 qubits), error mitigation, qLDPC decoding, noise learning, circuit-equivalence proofs, a phasor interference-ML layer, information-theoretic limits, noisy channels and state tomography, and signed energy-accounted receipts. No QPU, no cloud, no system libraries — identical results native, in the browser, and as a WASI component at the edge." +description = "A deterministic quantum stack in pure Rust: byte-exact circuit simulation (statevector / stabilizer / tensor-network MPS / sparse-Pauli backends), sparse Pauli dynamics at utility scale (arbitrary angles, 1024 qubits), belief-propagation tensor networks on the hardware graph, error mitigation, qLDPC decoding, noise learning, circuit-equivalence proofs, a phasor interference-ML layer, information-theoretic limits, noisy channels and state tomography, and signed energy-accounted receipts. No QPU, no cloud, no system libraries — identical results native, in the browser, and as a WASI component at the edge." repository = "https://git.transaction.science/Transaction-Science/open-standards" homepage = "https://wai.transaction.science/quantum-ide" documentation = "https://docs.rs/wai-quantum" @@ -45,6 +45,9 @@ quantum_pauli = ["quantum"] quantum_spd = [] # Signed, energy-accounted receipts for it, re-checkable bit for bit anywhere. quantum_spd_receipt = ["quantum_spd", "_signing"] +# Tensor-network states on the hardware graph: simple-update evolution, belief- +# propagation environments; exact on trees, an independent check on loopy graphs. +quantum_bptn = ["quantum_spd"] # Mixed states and the observables that distinguish quantum theory from a classical # account: density matrices, the partial trace, exact Pauli expectations, CHSH. quantum_info = ["quantum"] @@ -84,7 +87,7 @@ quantum_phasor_meter = ["_signing"] # Everything. full = [ - "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel", + "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel", "quantum_receipt", "quantum_ops", "quantum_cal", "quantum_mitigate", "quantum_qec", "quantum_control", "quantum_noise", "quantum_compile", "quantum_atom", "quantum_qbom", "quantum_vml", "quantum_phasor", "quantum_qfhrr", "quantum_kernel", "quantum_qdata", @@ -99,6 +102,10 @@ required-features = ["quantum_spd"] name = "utility_mz" required-features = ["quantum_spd"] +[[example]] +name = "utility_bptn" +required-features = ["quantum_bptn"] + [[example]] name = "bell" required-features = ["quantum", "quantum_toolchain"] diff --git a/wai/wai-quantum/README.md b/wai/wai-quantum/README.md index 282470b08..9359ff085 100644 --- a/wai/wai-quantum/README.md +++ b/wai/wai-quantum/README.md @@ -61,6 +61,7 @@ never creates information. | `quantum_pauli` | sparse-Pauli / Heisenberg observable propagation over the dyadic gate set (up to 64 qubits) | | `quantum_spd` | **sparse Pauli dynamics at utility scale**: arbitrary-angle Pauli rotations on up to 1024 qubits, exact Clifford and quarter-turn rotations, a rigorous bound on what truncation can move a result, and the 127-qubit heavy-hex kicked-Ising circuits — whose exactly computed 5-step magnetization it reproduces on all 127 qubits | | `quantum_spd_receipt` | a signed, energy-accounted receipt over a sparse-Pauli-dynamics run, re-checkable bit for bit on any machine | +| `quantum_bptn` | **tensor-network states on the hardware graph**: simple-update evolution with belief-propagation environments and re-gauging, the device's couplings applied a layer at a time in parallel (bit-identical to one at a time). Exact on trees; on the 127-qubit lattice it gives the exactly computed 5-step magnetization to 4e-7 at bond dimension 32 — an independent check on sparse Pauli dynamics, whose errors it does not share | | `quantum_receipt` | signed, energy-accounted receipt over a reconstruction | | `quantum_ops` | the operations/attestation layer | | `quantum_cal`, `quantum_control`, `quantum_noise` | calibration, filter-function robust control, DD noise spectroscopy, Cycle-Benchmarking noise learning | @@ -108,7 +109,9 @@ computes ⟨Z₆₂⟩ after 20 Trotter steps of the 127-qubit heavy-hex kicked- circuit in a fraction of a second, with the bound on how far truncation can have moved it. `--example utility_mz` checks the 5-step magnetization against the exactly computed values for every angle, and `tests/utility_reference.rs` holds the -published references. +published references. `quantum_bptn` computes the same circuits as tensor +networks (`--example utility_bptn`), so a result can be checked by a second +method that does not share the first one's approximation. ## Determinism diff --git a/wai/wai-quantum/examples/utility_bptn.rs b/wai/wai-quantum/examples/utility_bptn.rs new file mode 100644 index 000000000..2e7cdfc9a --- /dev/null +++ b/wai/wai-quantum/examples/utility_bptn.rs @@ -0,0 +1,52 @@ +//! The 127-qubit heavy-hex kicked-Ising experiment by belief-propagation +//! tensor networks, against published values. +//! +//! `cargo run --release --example utility_bptn --features quantum_bptn -- [chi ...]` +//! +//! For each bond dimension χ: the 5-step magnetization against its exactly +//! computed values (arXiv:2308.05077), and ⟨Z₆₂⟩ after 20 steps against the +//! published belief-propagation tensor-network values at the same angles +//! (arXiv:2306.14887, extrapolated in χ). + +use std::time::Instant; +use wai_quantum::quantum_bptn::Tns; +use wai_quantum::quantum_spd::{Pauli, heavy_hex_127, kicked_ising}; + +/// (k, exact 5-step Mz at θ_h = kπ/32) +const MZ5: [(u32, f64); 4] = [(2, 0.979836), (4, 0.904843), (6, 0.731464), (8, 0.452804)]; +/// (θ_h, published 20-step ⟨Z₆₂⟩, extrapolated in χ) +const Z62: [(f64, f64); 4] = [(0.2, 0.9786), (0.4, 0.8498), (0.6, 0.6686), (0.8, 0.2267)]; + +fn main() { + let chis: Vec = std::env::args().skip(1).map(|a| a.parse().expect("a bond dimension")).collect(); + let chis = if chis.is_empty() { vec![8, 16, 32] } else { chis }; + let edges = heavy_hex_127(); + for chi in chis { + println!("== χ = {chi}"); + for (k, exact) in MZ5 { + let started = Instant::now(); + let theta = f64::from(k) * std::f64::consts::PI / 32.0; + let mut t = Tns::new(127, &edges, chi).unwrap(); + t.run(&kicked_ising(127, &edges, theta, 5)).unwrap(); + let mz: f64 = (0..127).map(|q| t.expect1(q, Pauli::Z).unwrap()).sum::() / 127.0; + println!( + " 5-step Mz θ_h=kπ/32 k={k} {mz:.6} exact {exact:.6} |err| {:.1e} discarded {:.1e} {:.1}s", + (mz - exact).abs(), + t.discarded(), + started.elapsed().as_secs_f64() + ); + } + for (theta, published) in Z62 { + let started = Instant::now(); + let mut t = Tns::new(127, &edges, chi).unwrap(); + t.run(&kicked_ising(127, &edges, theta, 20)).unwrap(); + let z = t.expect1(62, Pauli::Z).unwrap(); + println!( + " 20-step Z62 θ_h={theta} {z:.4} published {published:.4} |diff| {:.1e} discarded {:.1e} {:.1}s", + (z - published).abs(), + t.discarded(), + started.elapsed().as_secs_f64() + ); + } + } +} diff --git a/wai/wai-quantum/src/lib.rs b/wai/wai-quantum/src/lib.rs index 3b77cb02c..59ef74f18 100644 --- a/wai/wai-quantum/src/lib.rs +++ b/wai/wai-quantum/src/lib.rs @@ -17,6 +17,9 @@ #[cfg(feature = "_signing")] mod merkle; +#[cfg(any(feature = "quantum_mps", feature = "quantum_bptn"))] +mod linalg; + #[cfg(any( feature = "quantum_spd", feature = "quantum_pauli", @@ -78,6 +81,9 @@ pub mod quantum_spd; #[cfg(feature = "quantum_spd_receipt")] pub mod quantum_spd_receipt; +#[cfg(feature = "quantum_bptn")] +pub mod quantum_bptn; + #[cfg(feature = "quantum_info")] pub mod quantum_info; diff --git a/wai/wai-quantum/src/linalg.rs b/wai/wai-quantum/src/linalg.rs new file mode 100644 index 000000000..48e2f95a7 --- /dev/null +++ b/wai/wai-quantum/src/linalg.rs @@ -0,0 +1,416 @@ +//! Dense complex linear algebra shared by the tensor-network backends: a +//! complex number, a one-sided Jacobi SVD and a Householder QR. Only +//! `+ − × ÷ √`, in a fixed order, so every result is the same on every machine. + +#![allow(clippy::should_implement_trait, clippy::needless_range_loop)] + +// =========================================================================== +// Complex f64 (deterministic: only + - * / sqrt, no FMA/transcendentals) +// =========================================================================== + +#[derive(Clone, Copy, Debug)] +pub struct C { + pub re: f64, + pub im: f64, +} +impl C { + pub(crate) const ZERO: C = C { re: 0.0, im: 0.0 }; + pub(crate) const ONE: C = C { re: 1.0, im: 0.0 }; + #[inline] + pub(crate) fn new(re: f64, im: f64) -> C { + C { re, im } + } + #[inline] + pub(crate) fn add(self, o: C) -> C { + C { re: self.re + o.re, im: self.im + o.im } + } + #[inline] + pub(crate) fn sub(self, o: C) -> C { + C { re: self.re - o.re, im: self.im - o.im } + } + #[inline] + pub(crate) fn mul(self, o: C) -> C { + C { re: self.re * o.re - self.im * o.im, im: self.re * o.im + self.im * o.re } + } + #[inline] + pub(crate) fn conj(self) -> C { + C { re: self.re, im: -self.im } + } + #[inline] + pub(crate) fn scale(self, s: f64) -> C { + C { re: self.re * s, im: self.im * s } + } + #[inline] + pub(crate) fn norm2(self) -> f64 { + self.re * self.re + self.im * self.im + } +} + +// =========================================================================== +// Complex one-sided Jacobi SVD (M = U · diag(S) · V†) +// =========================================================================== + +/// SVD of an `rows × cols` matrix given by its **columns** (`m[c]` is column `c`, a +/// length-`rows` vector). Returns `(u_cols, s, v_cols)` where each `u_cols[k]` is a +/// length-`rows` left singular vector, `s[k]` its singular value, `v_cols[k]` a +/// length-`cols` right singular vector — sorted by `s` descending. Deterministic: +/// fixed sweep order, fixed sweep cap, only `+ − × ÷ √`. +#[cfg(feature = "quantum_mps")] +pub(crate) fn jacobi_svd(mut m: Vec>, rows: usize, cols: usize) -> (Vec>, Vec, Vec>) { + // V starts as identity (cols × cols), stored by columns. + let mut v: Vec> = (0..cols) + .map(|j| (0..cols).map(|i| if i == j { C::ONE } else { C::ZERO }).collect()) + .collect(); + let eps = 1e-14; + for _sweep in 0..60 { + let mut off = 0.0f64; + for p in 0..cols { + for q in (p + 1)..cols { + // gram entries of columns p,q + let mut alpha = 0.0; // + let mut beta = 0.0; // + let mut gamma = C::ZERO; // = sum conj(c_p)·c_q + for k in 0..rows { + alpha += m[p][k].norm2(); + beta += m[q][k].norm2(); + gamma = gamma.add(m[p][k].conj().mul(m[q][k])); + } + let g2 = gamma.norm2(); + off += g2; + if g2 <= eps * alpha * beta || g2 == 0.0 { + continue; + } + let gabs = g2.sqrt(); + // phase-rotate column q so becomes real +gabs + let ph = C::new(gamma.re / gabs, gamma.im / gabs); // unit + let cph = ph.conj(); + for k in 0..rows { + m[q][k] = m[q][k].mul(cph); + } + for k in 0..cols { + v[q][k] = v[q][k].mul(cph); + } + // real Jacobi rotation with a=alpha,b=beta,g=gabs + let tau = (beta - alpha) / (2.0 * gabs); + let t = if tau >= 0.0 { + 1.0 / (tau + (1.0 + tau * tau).sqrt()) + } else { + -1.0 / (-tau + (1.0 + tau * tau).sqrt()) + }; + let cs = 1.0 / (1.0 + t * t).sqrt(); + let sn = t * cs; + for k in 0..rows { + let x = m[p][k]; + let y = m[q][k]; + m[p][k] = x.scale(cs).sub(y.scale(sn)); + m[q][k] = x.scale(sn).add(y.scale(cs)); + } + for k in 0..cols { + let x = v[p][k]; + let y = v[q][k]; + v[p][k] = x.scale(cs).sub(y.scale(sn)); + v[q][k] = x.scale(sn).add(y.scale(cs)); + } + } + } + if off <= eps { + break; + } + } + // singular values = column norms; left vectors = normalized columns + let mut s = vec![0.0; cols]; + for j in 0..cols { + s[j] = (0..rows).map(|k| m[j][k].norm2()).sum::().sqrt(); + } + let mut u: Vec> = m; + for j in 0..cols { + if s[j] > 1e-300 { + let inv = 1.0 / s[j]; + for k in 0..rows { + u[j][k] = u[j][k].scale(inv); + } + } + } + // sort by singular value descending (deterministic) + let mut order: Vec = (0..cols).collect(); + order.sort_by(|&a, &b| s[b].partial_cmp(&s[a]).unwrap_or(std::cmp::Ordering::Equal)); + let s2: Vec = order.iter().map(|&i| s[i]).collect(); + let u2: Vec> = order.iter().map(|&i| u[i].clone()).collect(); + let v2: Vec> = order.iter().map(|&i| v[i].clone()).collect(); + (u2, s2, v2) +} + +/// [`jacobi_svd`] run to convergence: a pair of columns is rotated until +/// their correlation is at rounding level, and the sweeps stop when one +/// rotates nothing. The original keeps its looser stopping rule (columns +/// orthogonal to about 1e-7 relative) because the MPS backend's results are +/// pinned to it. +#[cfg(feature = "quantum_bptn")] +pub(crate) fn jacobi_svd_strict(mut m: Vec>, rows: usize, cols: usize) -> (Vec>, Vec, Vec>) { + // V starts as identity (cols × cols), stored by columns. + let mut v: Vec> = (0..cols) + .map(|j| (0..cols).map(|i| if i == j { C::ONE } else { C::ZERO }).collect()) + .collect(); + // |⟨c_p,c_q⟩|² ≤ eps·‖c_p‖²‖c_q‖²: columns orthogonal to rounding. + let eps = 1e-30; + // A column this small against the whole matrix is numerically null: it + // holds no singular value anyone keeps, and rotating it only feeds + // rounding noise into V. + let null = 1e-30 * m.iter().flatten().map(|x| x.norm2()).sum::(); + for _sweep in 0..100 { + let mut rotated = false; + for p in 0..cols { + for q in (p + 1)..cols { + // gram entries of columns p,q + let mut alpha = 0.0; // + let mut beta = 0.0; // + let mut gamma = C::ZERO; // = sum conj(c_p)·c_q + for k in 0..rows { + alpha += m[p][k].norm2(); + beta += m[q][k].norm2(); + gamma = gamma.add(m[p][k].conj().mul(m[q][k])); + } + if alpha <= null || beta <= null { + continue; + } + let g2 = gamma.norm2(); + if g2 <= eps * alpha * beta || g2 == 0.0 { + continue; + } + rotated = true; + // |γ| scaled by its larger component, so a tiny γ keeps its + // phase at unit modulus without underflowing. Not `hypot`: + // libm's last bit differs between platforms; `/` and `√` do not. + let big = gamma.re.abs().max(gamma.im.abs()); + let (r, i) = (gamma.re / big, gamma.im / big); + let gabs = big * (r * r + i * i).sqrt(); + // phase-rotate column q so becomes real +gabs + let ph = C::new(gamma.re / gabs, gamma.im / gabs); // unit + let cph = ph.conj(); + for k in 0..rows { + m[q][k] = m[q][k].mul(cph); + } + for k in 0..cols { + v[q][k] = v[q][k].mul(cph); + } + // real Jacobi rotation with a=alpha,b=beta,g=gabs + let tau = (beta - alpha) / (2.0 * gabs); + let t = if tau >= 0.0 { + 1.0 / (tau + (1.0 + tau * tau).sqrt()) + } else { + -1.0 / (-tau + (1.0 + tau * tau).sqrt()) + }; + let cs = 1.0 / (1.0 + t * t).sqrt(); + let sn = t * cs; + for k in 0..rows { + let x = m[p][k]; + let y = m[q][k]; + m[p][k] = x.scale(cs).sub(y.scale(sn)); + m[q][k] = x.scale(sn).add(y.scale(cs)); + } + for k in 0..cols { + let x = v[p][k]; + let y = v[q][k]; + v[p][k] = x.scale(cs).sub(y.scale(sn)); + v[q][k] = x.scale(sn).add(y.scale(cs)); + } + } + } + if !rotated { + break; + } + } + // singular values = column norms; left vectors = normalized columns + let mut s = vec![0.0; cols]; + for j in 0..cols { + s[j] = (0..rows).map(|k| m[j][k].norm2()).sum::().sqrt(); + } + let mut u: Vec> = m; + for j in 0..cols { + if s[j] > 1e-300 { + let inv = 1.0 / s[j]; + for k in 0..rows { + u[j][k] = u[j][k].scale(inv); + } + } + } + // sort by singular value descending (deterministic) + let mut order: Vec = (0..cols).collect(); + order.sort_by(|&a, &b| s[b].partial_cmp(&s[a]).unwrap_or(std::cmp::Ordering::Equal)); + let s2: Vec = order.iter().map(|&i| s[i]).collect(); + let u2: Vec> = order.iter().map(|&i| u[i].clone()).collect(); + let v2: Vec> = order.iter().map(|&i| v[i].clone()).collect(); + (u2, s2, v2) +} + + +// =========================================================================== +// Householder QR (A = Q · R) +// =========================================================================== + +/// Thin QR of the row-major `m × n` matrix `a`: `Q` (`m × k`, orthonormal +/// columns) and `R` (`k × n`, upper triangular), both row-major, with +/// `k = min(m, n)`. Householder reflections in a fixed order; deterministic. +#[cfg(feature = "quantum_bptn")] +pub(crate) fn qr(a: &[C], m: usize, n: usize) -> (Vec, Vec, usize) { + let k = m.min(n); + let mut r = a.to_vec(); + let mut vs: Vec> = Vec::with_capacity(k); + for j in 0..k { + let norm = (j..m).map(|i| r[i * n + j].norm2()).sum::().sqrt(); + let mut v: Vec = (j..m).map(|i| r[i * n + j]).collect(); + if norm == 0.0 { + vs.push(vec![C::ZERO; m - j]); + continue; + } + // alpha = −(x₀/|x₀|)·‖x‖, so the reflection maps x to alpha·e₁ + let x0 = v[0]; + let a0 = x0.norm2().sqrt(); + let phase = if a0 == 0.0 { C::ONE } else { C::new(x0.re / a0, x0.im / a0) }; + let alpha = phase.scale(-norm); + v[0] = v[0].sub(alpha); + let vn = v.iter().map(|x| x.norm2()).sum::().sqrt(); + if vn == 0.0 { + vs.push(vec![C::ZERO; m - j]); + continue; + } + for x in v.iter_mut() { + *x = x.scale(1.0 / vn); + } + // r[j.., j..] -= 2 v (v† r[j.., j..]) + for c in j..n { + let mut dot = C::ZERO; + for (t, vi) in v.iter().enumerate() { + dot = dot.add(vi.conj().mul(r[(j + t) * n + c])); + } + let d2 = dot.scale(2.0); + for (t, vi) in v.iter().enumerate() { + r[(j + t) * n + c] = r[(j + t) * n + c].sub(vi.mul(d2)); + } + } + vs.push(v); + } + // Q = H₀ H₁ … H_{k−1} applied to the first k columns of the identity. + let mut q = vec![C::ZERO; m * k]; + for i in 0..k { + q[i * k + i] = C::ONE; + } + for j in (0..k).rev() { + let v = &vs[j]; + for c in 0..k { + let mut dot = C::ZERO; + for (t, vi) in v.iter().enumerate() { + dot = dot.add(vi.conj().mul(q[(j + t) * k + c])); + } + let d2 = dot.scale(2.0); + for (t, vi) in v.iter().enumerate() { + q[(j + t) * k + c] = q[(j + t) * k + c].sub(vi.mul(d2)); + } + } + } + let mut rr = vec![C::ZERO; k * n]; + for i in 0..k { + for c in i..n { + rr[i * n + c] = r[i * n + c]; + } + } + (q, rr, k) +} + +#[cfg(all(test, feature = "quantum_bptn"))] +mod tests { + use super::*; + + #[test] + fn qr_reconstructs_and_q_is_orthonormal() { + let mut s = 7u64; + let mut next = || { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + (s >> 11) as f64 / (1u64 << 53) as f64 - 0.5 + }; + for (m, n) in [(5, 3), (3, 5), (8, 8), (1, 4), (6, 1)] { + let a: Vec = (0..m * n).map(|_| C::new(next(), next())).collect(); + let (q, r, k) = qr(&a, m, n); + for i in 0..m { + for c in 0..n { + let mut x = C::ZERO; + for t in 0..k { + x = x.add(q[i * k + t].mul(r[t * n + c])); + } + let d = x.sub(a[i * n + c]); + assert!(d.norm2() < 1e-26, "({m}x{n}) entry {i},{c}"); + } + } + for a1 in 0..k { + for b1 in 0..k { + let mut x = C::ZERO; + for i in 0..m { + x = x.add(q[i * k + a1].conj().mul(q[i * k + b1])); + } + let want = if a1 == b1 { 1.0 } else { 0.0 }; + assert!((x.re - want).abs() < 1e-13 && x.im.abs() < 1e-13); + } + } + for i in 0..k { + for c in 0..i.min(n) { + assert_eq!(r[i * n + c].norm2(), 0.0, "R is upper triangular"); + } + } + } + } +} + +#[cfg(all(test, feature = "quantum_bptn"))] +mod svd_tests { + use super::*; + + /// The converged SVD reconstructs `A = U·S·V†` and returns orthonormal + /// singular vectors, both to rounding, on full-rank and rank-deficient + /// matrices. (The original's looser stopping rule leaves `U` orthonormal + /// only to about 1e-7.) + #[test] + fn the_converged_svd_reconstructs_with_orthonormal_vectors() { + let mut s = 99u64; + let mut next = || { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + (s >> 11) as f64 / (1u64 << 53) as f64 - 0.5 + }; + for trial in 0..200 { + let rows = 2 + trial % 9; + let cols = 2 + (trial / 9) % 7; + let rank = 1 + trial % rows.min(cols); + let b: Vec = (0..rows * rank).map(|_| C::new(next(), next())).collect(); + let c: Vec = (0..cols * rank).map(|_| C::new(next(), next())).collect(); + let a: Vec = (0..rows * cols) + .map(|i| { + let (r, q) = (i / cols, i % cols); + (0..rank).fold(C::ZERO, |acc, k| acc.add(b[r * rank + k].mul(c[q * rank + k]))) + }) + .collect(); + let columns: Vec> = (0..cols).map(|q| (0..rows).map(|r| a[r * cols + q]).collect()).collect(); + let (u, sv, v) = jacobi_svd_strict(columns, rows, cols); + for r in 0..rows { + for q in 0..cols { + let x = (0..cols).fold(C::ZERO, |acc, k| acc.add(u[k][r].scale(sv[k]).mul(v[k][q].conj()))); + assert!(x.sub(a[r * cols + q]).norm2() < 1e-26, "trial {trial}: reconstruction"); + } + } + for k1 in 0..cols { + for k2 in 0..cols { + for (left, vecs, len) in [(true, &u, rows), (false, &v, cols)] { + if left && (sv[k1] < 1e-8 * sv[0] || sv[k2] < 1e-8 * sv[0]) { + continue; // a null column of U is not normalised + } + let d = (0..len).fold(C::ZERO, |acc, r| acc.add(vecs[k1][r].conj().mul(vecs[k2][r]))); + let want = if k1 == k2 { 1.0 } else { 0.0 }; + assert!((d.re - want).abs() < 1e-13 && d.im.abs() < 1e-13, "trial {trial}: orthonormality"); + } + } + } + } + } +} diff --git a/wai/wai-quantum/src/quantum_bptn.rs b/wai/wai-quantum/src/quantum_bptn.rs new file mode 100644 index 000000000..6155937ac --- /dev/null +++ b/wai/wai-quantum/src/quantum_bptn.rs @@ -0,0 +1,1100 @@ +//! Tensor-network states on the hardware graph, evolved by simple update and +//! read through belief-propagation environments (`wai.quantum.bptn`). +//! +//! The state is a tensor network laid on the device's coupling graph: one +//! tensor per qubit, with a physical index and one virtual index per edge. It +//! is kept in Vidal gauge, with a vector of bond weights `λ_e` on every edge. A +//! single-qubit gate acts on one tensor. A two-qubit gate on an edge acts on +//! the two tensors and the bond between them: +//! +//! 1. Each tensor, with the weights of its other bonds absorbed, is QR-reduced. +//! 2. The gate is applied to the two small factors and the bond weight between +//! them. +//! 3. An SVD splits the result again, keeping at most `χ` singular values. +//! These become the new bond weights. +//! +//! An expectation value takes each tensor's environment to be the product of +//! its bond weights. That is the rank-one environment belief propagation +//! converges to, and on a tree it is exact. +//! +//! What to expect from it: +//! +//! - **On a tree** (a chain, or any loop-free coupling graph), with `χ` large +//! enough to truncate nothing, the evolution and the expectation values are +//! exact up to rounding. The tests hold it to sparse Pauli dynamics there. +//! - **On a graph with loops**, such as heavy-hex with its 12-qubit rings, the +//! environments neglect correlations around the loops. Truncation discards +//! weight as well, and [`Tns::discarded`] reports how much. +//! +//! Its errors are unrelated to those of sparse Pauli dynamics, which is what +//! makes it worth having beside that method: where the two agree, neither +//! method's approximation is likely to be what they agree on. +//! +//! Arithmetic: `+ − × ÷ √` in a fixed order (Householder QR, one-sided Jacobi +//! SVD), with `sin` and `cos` computed in-crate. Results are bit-identical on +//! every platform. + +#![allow(clippy::needless_range_loop)] + +use crate::linalg::{C, jacobi_svd_strict, qr}; +use crate::quantum_spd::{Op, Pauli, RotCircuit}; +use crate::repro::sin_cos; +use std::fmt; + +/// Why a circuit could not be run on a graph. +#[derive(Clone, Debug, PartialEq)] +pub enum TnsError { + /// A qubit outside the graph. + QubitOutOfRange(u32), + /// A two-qubit gate on qubits the graph does not connect. + NotAnEdge(u32, u32), + /// A rotation about a string of more than two qubits. + TooWide(usize), + /// A circuit `quantum_spd` would refuse. + Invalid(String), +} + +impl fmt::Display for TnsError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + TnsError::QubitOutOfRange(q) => write!(f, "qubit {q} is not on the graph"), + TnsError::NotAnEdge(a, b) => write!(f, "no edge between {a} and {b}"), + TnsError::TooWide(w) => write!(f, "a rotation on {w} qubits; at most 2 here"), + TnsError::Invalid(e) => write!(f, "{e}"), + } + } +} + +impl std::error::Error for TnsError {} + +/// The state and how it was truncated. +pub struct Tns { + n: usize, + edges: Vec<(u32, u32)>, + /// Per vertex: incident edge ids, ascending. A tensor's virtual indices + /// follow this order. + nbrs: Vec>, + /// Per vertex: physical index outermost, then the virtual indices in + /// `nbrs` order, row-major. + t: Vec>, + /// Per edge: the bond weights, descending, with unit 2-norm. + lam: Vec>, + chi: usize, + discarded: f64, +} + +/// A singular value is dropped below this fraction of the largest. Dividing +/// a tensor by a smaller weight would amplify rounding past use. +const CUTOFF: f64 = 1e-10; + +impl Tns { + /// `|0…0⟩` on the graph with these edges, keeping at most `chi` singular + /// values per bond. + pub fn new(n: u32, edges: &[(u32, u32)], chi: usize) -> Result { + let n = n as usize; + let mut nbrs = vec![Vec::new(); n]; + for (e, &(a, b)) in edges.iter().enumerate() { + for q in [a, b] { + if q as usize >= n { + return Err(TnsError::QubitOutOfRange(q)); + } + } + if a == b { + return Err(TnsError::NotAnEdge(a, b)); + } + nbrs[a as usize].push(e); + nbrs[b as usize].push(e); + } + // Every bond starts at dimension 1, so a tensor holds two numbers: + // amplitude 1 for |0⟩, 0 for |1⟩. + let t = (0..n).map(|_| vec![C::ONE, C::ZERO]).collect(); + Ok(Tns { n, edges: edges.to_vec(), nbrs, t, lam: vec![vec![1.0]; edges.len()], chi: chi.max(1), discarded: 0.0 }) + } + + /// The summed fraction of squared singular-value weight the truncations + /// have discarded: 0 means no truncation has happened yet. + pub fn discarded(&self) -> f64 { + self.discarded + } + + /// The largest bond dimension in the network. + pub fn max_bond(&self) -> usize { + self.lam.iter().map(Vec::len).max().unwrap_or(1) + } + + fn bond_dims(&self, v: usize) -> Vec { + self.nbrs[v].iter().map(|&e| self.lam[e].len()).collect() + } + + fn edge_of(&self, a: u32, b: u32) -> Option { + self.nbrs.get(a as usize)?.iter().copied().find(|&e| { + let (x, y) = self.edges[e]; + (x, y) == (a, b) || (x, y) == (b, a) + }) + } + + /// Apply every op of `circuit` in order. + /// + /// Consecutive two-qubit gates on vertex-disjoint edges, such as a layer of + /// the device's couplings, touch disjoint tensors and bonds and read no + /// bond another of them writes. They are computed in parallel and + /// committed in order, which gives the same bits as applying them one by + /// one (and is what happens on targets without threads). + pub fn run(&mut self, circuit: &RotCircuit) -> Result<(), TnsError> { + self.run_inner(circuit, None) + } + + /// [`Self::run`], re-gauging the network into the belief-propagation gauge + /// (at most `sweeps` sweeps to tolerance `tol`) after every layer of + /// two-qubit gates, so each truncation is made against the environment + /// belief propagation sees rather than the bond weights alone. + pub fn run_gauged(&mut self, circuit: &RotCircuit, sweeps: usize, tol: f64) -> Result<(), TnsError> { + self.run_inner(circuit, Some((sweeps, tol))) + } + + fn run_inner(&mut self, circuit: &RotCircuit, gauge: Option<(usize, f64)>) -> Result<(), TnsError> { + circuit.validate().map_err(|e| TnsError::Invalid(e.to_string()))?; + if circuit.n as usize != self.n { + return Err(TnsError::Invalid(format!("a {}-qubit circuit on a {}-qubit graph", circuit.n, self.n))); + } + let mut batch: Vec<(u32, u32, [[C; 4]; 4])> = Vec::new(); + let mut busy = vec![false; self.n]; + for op in &circuit.ops { + match gate_of(op)? { + Gate::One(q, g) => { + if busy[q as usize] { + self.flush(&mut batch, &mut busy, gauge)?; + } + self.gate1(q, g)?; + } + Gate::Two(a, b, g) => { + for q in [a, b] { + if q as usize >= self.n { + return Err(TnsError::QubitOutOfRange(q)); + } + } + if busy[a as usize] || busy[b as usize] { + self.flush(&mut batch, &mut busy, gauge)?; + } + busy[a as usize] = true; + busy[b as usize] = true; + batch.push((a, b, g)); + } + } + } + self.flush(&mut batch, &mut busy, gauge) + } + + fn flush( + &mut self, + batch: &mut Vec<(u32, u32, [[C; 4]; 4])>, + busy: &mut [bool], + gauge: Option<(usize, f64)>, + ) -> Result<(), TnsError> { + if batch.is_empty() { + return Ok(()); + } + let threads = if cfg!(target_arch = "wasm32") { + 1 + } else { + std::thread::available_parallelism().map_or(1, |n| n.get()).min(batch.len()) + }; + let updates: Vec> = if threads <= 1 { + batch.iter().map(|&(a, b, g)| self.compute2(a, b, g)).collect() + } else { + let this = &*self; + let chunk = batch.len().div_ceil(threads); + std::thread::scope(|scope| { + let handles: Vec<_> = batch + .chunks(chunk) + .map(|part| scope.spawn(move || part.iter().map(|&(a, b, g)| this.compute2(a, b, g)).collect::>())) + .collect(); + handles.into_iter().flat_map(|h| h.join().expect("a gate's computation panicked")).collect() + }) + }; + for up in updates { + self.commit(up?); + } + for &(a, b, _) in batch.iter() { + busy[a as usize] = false; + busy[b as usize] = false; + } + batch.clear(); + if let Some((sweeps, tol)) = gauge { + let env = self.belief_propagation(sweeps, tol); + self.regauge(&env); + } + Ok(()) + } + + /// Apply one op. + pub fn apply(&mut self, op: &Op) -> Result<(), TnsError> { + match gate_of(op)? { + Gate::One(q, g) => self.gate1(q, g), + Gate::Two(a, b, g) => { + let up = self.compute2(a, b, g)?; + self.commit(up); + Ok(()) + } + } + } + + fn gate1(&mut self, q: u32, g: [[C; 2]; 2]) -> Result<(), TnsError> { + let v = q as usize; + if v >= self.n { + return Err(TnsError::QubitOutOfRange(q)); + } + let half = self.t[v].len() / 2; + let t = &mut self.t[v]; + for i in 0..half { + let (a, b) = (t[i], t[half + i]); + t[i] = g[0][0].mul(a).add(g[0][1].mul(b)); + t[half + i] = g[1][0].mul(a).add(g[1][1].mul(b)); + } + Ok(()) + } + + /// What a two-qubit gate `g` (row and column index `2·s_a + s_b`) on edge + /// `(a, b)` does, computed without changing the state: [`Self::commit`] + /// applies it. + fn compute2(&self, a: u32, b: u32, g: [[C; 4]; 4]) -> Result { + for q in [a, b] { + if q as usize >= self.n { + return Err(TnsError::QubitOutOfRange(q)); + } + } + let e = self.edge_of(a, b).ok_or(TnsError::NotAnEdge(a, b))?; + let (u, v) = (a as usize, b as usize); + let d = self.lam[e].len(); + + // 1. Each side as a matrix: rows over its other bonds (their weights + // absorbed), columns over (physical, this bond); then QR. + let (mu, ru, rest_u) = self.side_matrix(u, e); + let (mv, rv, rest_v) = self.side_matrix(v, e); + let (qu, rru, ku) = qr(&mu, ru, 2 * d); + let (qv, rrv, kv) = qr(&mv, rv, 2 * d); + + // 2. θ[(x,s),(t,y)] = Σ_i Ru[x,(s,i)] λ[i] Rv[y,(t,i)], then the gate. + let lam = &self.lam[e]; + let rows = 2 * ku; + let cols = 2 * kv; + let mut theta = vec![C::ZERO; rows * cols]; + for x in 0..ku { + for s in 0..2 { + for y in 0..kv { + for t in 0..2 { + let mut acc = C::ZERO; + for i in 0..d { + acc = acc.add(rru[x * 2 * d + s * d + i].scale(lam[i]).mul(rrv[y * 2 * d + t * d + i])); + } + theta[(x * 2 + s) * cols + t * kv + y] = acc; + } + } + } + } + let mut gated = vec![C::ZERO; rows * cols]; + for x in 0..ku { + for y in 0..kv { + for s2 in 0..2 { + for t2 in 0..2 { + let mut acc = C::ZERO; + for s in 0..2 { + for t in 0..2 { + acc = acc.add(g[2 * s2 + t2][2 * s + t].mul(theta[(x * 2 + s) * cols + t * kv + y])); + } + } + gated[(x * 2 + s2) * cols + t2 * kv + y] = acc; + } + } + } + } + + // 3. SVD, truncated. + let (left, sv, right) = svd_rows(&gated, rows, cols); + let total: f64 = sv.iter().map(|x| x * x).sum(); + let top = sv.first().copied().unwrap_or(0.0); + let keep = sv.iter().take(self.chi).take_while(|&&x| x > CUTOFF * top).count().max(1); + let kept: f64 = sv[..keep].iter().map(|x| x * x).sum(); + let discarded = if total > 0.0 { (total - kept) / total } else { 0.0 }; + let norm = kept.sqrt(); + let new_lam: Vec = sv[..keep].iter().map(|x| x / norm).collect(); + + // 4. New tensors: Q times the singular vectors, the other bonds' + // weights divided back out. + let mut new_u = vec![C::ZERO; rest_u * 2 * keep]; + for o in 0..rest_u { + for s in 0..2 { + for k in 0..keep { + let mut acc = C::ZERO; + for x in 0..ku { + acc = acc.add(qu[o * ku + x].mul(left[k][x * 2 + s])); + } + new_u[(o * 2 + s) * keep + k] = acc; + } + } + } + let mut new_v = vec![C::ZERO; rest_v * 2 * keep]; + for o in 0..rest_v { + for t in 0..2 { + for k in 0..keep { + let mut acc = C::ZERO; + for y in 0..kv { + acc = acc.add(qv[o * kv + y].mul(right[k][t * kv + y].conj())); + } + new_v[(o * 2 + t) * keep + k] = acc; + } + } + } + Ok(Update { + e, + u, + v, + tu: self.rebuilt(u, e, &new_u, keep), + tv: self.rebuilt(v, e, &new_v, keep), + lam: new_lam, + discarded, + }) + } + + fn commit(&mut self, up: Update) { + self.discarded += up.discarded; + self.lam[up.e] = up.lam; + self.t[up.u] = up.tu; + self.t[up.v] = up.tv; + } + + + /// Vertex `v`'s tensor as a matrix: rows over its bonds other than `e` + /// (row-major in `nbrs` order, each scaled by its bond weight), columns + /// over (physical, bond `e`). Returns the matrix, its row count, and that + /// count again for the caller's reshaping. + fn side_matrix(&self, v: usize, e: usize) -> (Vec, usize, usize) { + let dims = self.bond_dims(v); + let pos = self.nbrs[v].iter().position(|&x| x == e).expect("the edge is incident"); + let d = dims[pos]; + let bonds: usize = dims.iter().product(); + let rest = bonds / d; + let mut m = vec![C::ZERO; rest * 2 * d]; + let mut idx = vec![0usize; dims.len()]; + for flat in 0..bonds { + // unflatten + let mut r = flat; + for j in (0..dims.len()).rev() { + idx[j] = r % dims[j]; + r /= dims[j]; + } + let mut row = 0usize; + let mut w = 1.0; + for j in 0..dims.len() { + if j != pos { + row = row * dims[j] + idx[j]; + w *= self.lam[self.nbrs[v][j]][idx[j]]; + } + } + for s in 0..2 { + m[row * 2 * d + s * d + idx[pos]] = self.t[v][s * bonds + flat].scale(w); + } + } + (m, rest, rest) + } + + /// Vertex `v`'s new tensor from `m` (rows over the other bonds, as in + /// `side_matrix`, then physical, then the new bond `e` of size `keep`), + /// with the other bonds' weights divided back out. + fn rebuilt(&self, v: usize, e: usize, m: &[C], keep: usize) -> Vec { + let pos = self.nbrs[v].iter().position(|&x| x == e).expect("the edge is incident"); + let mut dims = self.bond_dims(v); + dims[pos] = keep; + let bonds: usize = dims.iter().product(); + let mut t = vec![C::ZERO; 2 * bonds]; + let mut idx = vec![0usize; dims.len()]; + for flat in 0..bonds { + let mut r = flat; + for j in (0..dims.len()).rev() { + idx[j] = r % dims[j]; + r /= dims[j]; + } + let mut row = 0usize; + let mut w = 1.0; + for j in 0..dims.len() { + if j != pos { + row = row * dims[j] + idx[j]; + w *= self.lam[self.nbrs[v][j]][idx[j]]; + } + } + let inv = 1.0 / w; + for s in 0..2 { + t[s * bonds + flat] = m[(row * 2 + s) * keep + idx[pos]].scale(inv); + } + } + t + } + + /// `⟨P⟩` for a Pauli on one qubit, from the tensor and the squared weights + /// of its bonds: the rank-one environment belief propagation converges to. + pub fn expect1(&self, q: u32, p: Pauli) -> Result { + let v = q as usize; + if v >= self.n { + return Err(TnsError::QubitOutOfRange(q)); + } + let dims = self.bond_dims(v); + let bonds: usize = dims.iter().product(); + let mut rho = [[C::ZERO; 2]; 2]; + let mut idx = vec![0usize; dims.len()]; + for flat in 0..bonds { + let mut r = flat; + for j in (0..dims.len()).rev() { + idx[j] = r % dims[j]; + r /= dims[j]; + } + let w: f64 = (0..dims.len()).map(|j| self.lam[self.nbrs[v][j]][idx[j]]).map(|l| l * l).product(); + for s in 0..2 { + for s2 in 0..2 { + let x = self.t[v][s * bonds + flat].mul(self.t[v][s2 * bonds + flat].conj()).scale(w); + rho[s][s2] = rho[s][s2].add(x); + } + } + } + let pm = pauli(p); + let mut num = C::ZERO; + for s in 0..2 { + for s2 in 0..2 { + num = num.add(rho[s][s2].mul(pm[s2][s])); + } + } + let den = rho[0][0].re + rho[1][1].re; + Ok(num.re / den) + } +} + +/// Belief-propagation messages on the norm network `⟨ψ|ψ⟩`. For each edge +/// `e = (a, b)` there is a message from `a` to `b` and one from `b` to `a`: +/// `χ_e × χ_e` matrices, the row on the ket's bond index and the column on the +/// bra's, normalised to unit trace. +pub struct BpEnv { + msgs: Vec<[Vec; 2]>, + /// Sweeps run. + pub sweeps: usize, + /// The largest change in any message entry during the last sweep. + pub residual: f64, +} + +impl Tns { + /// Vertex `v`'s tensor with the square root of each of its bond weights + /// absorbed, so that the norm network is the tensors and their conjugates + /// alone. + fn absorbed(&self, v: usize) -> Vec { + let dims = self.bond_dims(v); + let bonds: usize = dims.iter().product(); + let mut out = self.t[v].clone(); + let mut idx = vec![0usize; dims.len()]; + for flat in 0..bonds { + let mut r = flat; + for j in (0..dims.len()).rev() { + idx[j] = r % dims[j]; + r /= dims[j]; + } + let w: f64 = (0..dims.len()).map(|j| self.lam[self.nbrs[v][j]][idx[j]].sqrt()).product(); + for s in 0..2 { + out[s * bonds + flat] = out[s * bonds + flat].scale(w); + } + } + out + } + + /// The message arriving at `v` across its bond at position `j`. + fn incoming<'a>(&self, msgs: &'a [[Vec; 2]], v: usize, j: usize) -> &'a [C] { + let e = self.nbrs[v][j]; + if self.edges[e].0 as usize == v { &msgs[e][1] } else { &msgs[e][0] } + } + + /// `a` (vertex `v`'s absorbed tensor) with the incoming message contracted + /// into every bond but `skip`: each such index now runs over the bra side. + fn with_messages(&self, a: &[C], v: usize, msgs: &[[Vec; 2]], skip: Option) -> Vec { + let dims = self.bond_dims(v); + let bonds: usize = dims.iter().product(); + let mut cur = a.to_vec(); + for j in 0..dims.len() { + if skip == Some(j) { + continue; + } + let m = self.incoming(msgs, v, j); + let d = dims[j]; + let stride: usize = dims[j + 1..].iter().product(); + let mut next = vec![C::ZERO; cur.len()]; + for s in 0..2 { + for flat in 0..bonds { + let x = cur[s * bonds + flat]; + if x.re == 0.0 && x.im == 0.0 { + continue; + } + let i = (flat / stride) % d; + let base = flat - i * stride; + for ip in 0..d { + let k = s * bonds + base + ip * stride; + next[k] = next[k].add(x.mul(m[i * d + ip])); + } + } + } + cur = next; + } + cur + } + + /// The message from `u` across its bond at position `p`, given `msgs`. + fn message(&self, abs: &[Vec], msgs: &[[Vec; 2]], u: usize, p: usize) -> Vec { + let a = &abs[u]; + let c = self.with_messages(a, u, msgs, Some(p)); + let dims = self.bond_dims(u); + let bonds: usize = dims.iter().product(); + let d = dims[p]; + let stride: usize = dims[p + 1..].iter().product(); + let mut m = vec![C::ZERO; d * d]; + for s in 0..2 { + for base in 0..bonds { + if (base / stride) % d != 0 { + continue; + } + for i in 0..d { + let ci = c[s * bonds + base + i * stride]; + for ip in 0..d { + m[i * d + ip] = m[i * d + ip].add(ci.mul(a[s * bonds + base + ip * stride].conj())); + } + } + } + } + let tr: f64 = (0..d).map(|i| m[i * d + i].re).sum(); + if tr > 0.0 { + for x in m.iter_mut() { + *x = x.scale(1.0 / tr); + } + } + m + } + + /// Run belief propagation on the norm network to a fixed point: at most + /// `max_sweeps` sweeps, stopping once no message entry changes by more + /// than `tol`. Each sweep computes every message from the previous + /// sweep's (in parallel where there are threads), so the result is the + /// same everywhere. + pub fn belief_propagation(&self, max_sweeps: usize, tol: f64) -> BpEnv { + let abs: Vec> = (0..self.n).map(|v| self.absorbed(v)).collect(); + let mut msgs: Vec<[Vec; 2]> = self + .lam + .iter() + .map(|l| { + // In Vidal form the consistent message is diag(λ): start there. + let d = l.len(); + let sum: f64 = l.iter().sum(); + let diag: Vec = (0..d * d).map(|k| if k % (d + 1) == 0 { C::new(l[k / (d + 1)] / sum, 0.0) } else { C::ZERO }).collect(); + [diag.clone(), diag] + }) + .collect(); + // Directed edges, in order: (edge, direction, sender, sender's bond position). + let directed: Vec<(usize, usize, usize, usize)> = (0..self.edges.len()) + .flat_map(|e| { + let (a, b) = (self.edges[e].0 as usize, self.edges[e].1 as usize); + let pa = self.nbrs[a].iter().position(|&x| x == e).expect("incident"); + let pb = self.nbrs[b].iter().position(|&x| x == e).expect("incident"); + [(e, 0, a, pa), (e, 1, b, pb)] + }) + .collect(); + let mut sweeps = 0; + let mut residual = f64::INFINITY; + while sweeps < max_sweeps && residual > tol { + let threads = if cfg!(target_arch = "wasm32") { + 1 + } else { + std::thread::available_parallelism().map_or(1, |n| n.get()).min(directed.len()) + }; + let fresh: Vec> = if threads <= 1 { + directed.iter().map(|&(_, _, u, p)| self.message(&abs, &msgs, u, p)).collect() + } else { + let (this, abs_ref, msgs_ref) = (&*self, &abs, &msgs); + let chunk = directed.len().div_ceil(threads); + std::thread::scope(|scope| { + let handles: Vec<_> = directed + .chunks(chunk) + .map(|part| scope.spawn(move || part.iter().map(|&(_, _, u, p)| this.message(abs_ref, msgs_ref, u, p)).collect::>())) + .collect(); + handles.into_iter().flat_map(|h| h.join().expect("a message's computation panicked")).collect() + }) + }; + residual = 0.0; + for (&(e, dir, _, _), m) in directed.iter().zip(fresh) { + for (old, new) in msgs[e][dir].iter().zip(&m) { + residual = residual.max(old.sub(*new).norm2().sqrt()); + } + msgs[e][dir] = m; + } + sweeps += 1; + } + BpEnv { msgs, sweeps, residual } + } + + /// Bring the network into the gauge belief propagation's messages define: + /// on each edge, factor the two messages as `M = X X†`, take the SVD + /// `Xₐᵀ X_b = U Λ V†`, and put `(Xₐ⁺)ᵀ U` on one side and `(V† X_b⁺)ᵀ` on the + /// other. That inserts the identity on the bond (on the support the + /// messages see), so the state does not change, and `Λ` becomes the bond's + /// weights: afterwards [`Self::expect1`]'s environment is the BP one. + pub fn regauge(&mut self, env: &BpEnv) { + let abs: Vec> = (0..self.n).map(|v| self.absorbed(v)).collect(); + let mut transforms: Vec<[Vec; 2]> = Vec::with_capacity(self.edges.len()); + let mut new_lam: Vec> = Vec::with_capacity(self.edges.len()); + for e in 0..self.edges.len() { + let d = self.lam[e].len(); + let (xa, xap, ra) = psd_sqrt(&env.msgs[e][0], d); + let (xb, xbp, rb) = psd_sqrt(&env.msgs[e][1], d); + let mut bm = vec![C::ZERO; ra * rb]; + for p in 0..ra { + for q in 0..rb { + let mut acc = C::ZERO; + for i in 0..d { + acc = acc.add(xa[i * ra + p].mul(xb[i * rb + q])); + } + bm[p * rb + q] = acc; + } + } + let (u, sv, v) = svd_rows(&bm, ra, rb); + let top = sv.first().copied().unwrap_or(0.0); + let k = sv.iter().take_while(|&&x| x > CUTOFF * top).count().max(1); + let norm = sv[..k].iter().map(|x| x * x).sum::().sqrt(); + new_lam.push(sv[..k].iter().map(|x| x / norm).collect()); + let mut ta = vec![C::ZERO; d * k]; + let mut tb = vec![C::ZERO; d * k]; + for i in 0..d { + for n in 0..k { + let mut a_acc = C::ZERO; + for p in 0..ra { + a_acc = a_acc.add(xap[p * d + i].mul(u[n][p])); + } + ta[i * k + n] = a_acc; + let mut b_acc = C::ZERO; + for q in 0..rb { + b_acc = b_acc.add(v[n][q].conj().mul(xbp[q * d + i])); + } + tb[i * k + n] = b_acc; + } + } + transforms.push([ta, tb]); + } + for v in 0..self.n { + let mut dims = self.bond_dims(v); + let mut cur = abs[v].clone(); + for j in 0..dims.len() { + let e = self.nbrs[v][j]; + let side = if self.edges[e].0 as usize == v { 0 } else { 1 }; + let k = new_lam[e].len(); + cur = transform_bond(&cur, &dims, j, &transforms[e][side], k); + dims[j] = k; + } + self.t[v] = cur; + } + self.lam = new_lam; + } + + /// `⟨P⟩` for a Pauli on one qubit, with the environment `env` belief + /// propagation converged to. + pub fn expect1_bp(&self, env: &BpEnv, q: u32, p: Pauli) -> Result { + let v = q as usize; + if v >= self.n { + return Err(TnsError::QubitOutOfRange(q)); + } + let a = self.absorbed(v); + let c = self.with_messages(&a, v, &env.msgs, None); + let bonds = a.len() / 2; + let mut rho = [[C::ZERO; 2]; 2]; + for s in 0..2 { + for s2 in 0..2 { + for k in 0..bonds { + rho[s][s2] = rho[s][s2].add(c[s * bonds + k].mul(a[s2 * bonds + k].conj())); + } + } + } + let pm = pauli(p); + let mut num = C::ZERO; + for s in 0..2 { + for s2 in 0..2 { + num = num.add(rho[s][s2].mul(pm[s2][s])); + } + } + Ok(num.re / (rho[0][0].re + rho[1][1].re)) + } +} + +/// A gate as matrices: one qubit, or two (row and column index `2·s_a + s_b`). +enum Gate { + One(u32, [[C; 2]; 2]), + Two(u32, u32, [[C; 4]; 4]), +} + +/// What a two-qubit gate does to the state, ready to commit. +struct Update { + e: usize, + u: usize, + v: usize, + tu: Vec, + tv: Vec, + lam: Vec, + discarded: f64, +} + +fn gate_of(op: &Op) -> Result { + Ok(match op { + Op::H(q) | Op::S(q) | Op::Sdg(q) | Op::SX(q) | Op::SXdg(q) | Op::X(q) | Op::Y(q) | Op::Z(q) => Gate::One(*q, fixed1(op)), + Op::CX(a, b) => Gate::Two(*a, *b, cx()), + Op::CZ(a, b) => Gate::Two(*a, *b, cz()), + Op::Swap(a, b) => Gate::Two(*a, *b, swap()), + Op::Rot { axis, theta } => { + let (s, c) = sin_cos(0.5 * theta); + rotation(axis, c, s)? + } + Op::QuarterRot { axis, k } => { + // exp(−i(kπ/2)P/2): cos and sin of kπ/4, exactly. + let h = std::f64::consts::FRAC_1_SQRT_2; + let (c, s) = match k { + 1 => (h, h), + 2 => (0.0, 1.0), + _ => (-h, h), + }; + rotation(axis, c, s)? + } + }) +} + +/// `exp(−iθP/2) = cos(θ/2)·I − i·sin(θ/2)·P`, given `cos(θ/2), sin(θ/2)`. +fn rotation(axis: &[(u32, Pauli)], c: f64, s: f64) -> Result { + let minus_is = C::new(0.0, -s); + match axis { + [(q, p)] => { + let pm = pauli(*p); + let mut g = [[C::ZERO; 2]; 2]; + for r in 0..2 { + for col in 0..2 { + let id = if r == col { C::new(c, 0.0) } else { C::ZERO }; + g[r][col] = id.add(minus_is.mul(pm[r][col])); + } + } + Ok(Gate::One(*q, g)) + } + [(a, pa), (b, pb)] => { + let (ma, mb) = (pauli(*pa), pauli(*pb)); + let mut g = [[C::ZERO; 4]; 4]; + for r in 0..4 { + for col in 0..4 { + let kron = ma[r >> 1][col >> 1].mul(mb[r & 1][col & 1]); + let id = if r == col { C::new(c, 0.0) } else { C::ZERO }; + g[r][col] = id.add(minus_is.mul(kron)); + } + } + Ok(Gate::Two(*a, *b, g)) + } + _ => Err(TnsError::TooWide(axis.len())), + } +} + +/// A Hermitian positive semi-definite `d × d` matrix (row-major) as `X X†`: +/// `X` (`d × r`, row-major), its pseudo-inverse `X⁺` (`r × d`), and `r`, the +/// number of eigenvalues above the cutoff. +fn psd_sqrt(m: &[C], d: usize) -> (Vec, Vec, usize) { + let (u, s, _) = svd_rows(m, d, d); + let top = s.first().copied().unwrap_or(0.0); + let r = s.iter().take_while(|&&x| x > CUTOFF * top).count().max(1); + let mut x = vec![C::ZERO; d * r]; + let mut xp = vec![C::ZERO; r * d]; + for k in 0..r { + let root = s[k].sqrt(); + for i in 0..d { + x[i * r + k] = u[k][i].scale(root); + xp[k * d + i] = u[k][i].conj().scale(1.0 / root); + } + } + (x, xp, r) +} + +/// Contract bond `j` of a tensor (physical index outermost, then the bonds +/// row-major with dimensions `dims`) with the `dims[j] × k` matrix `m`. +fn transform_bond(t: &[C], dims: &[usize], j: usize, m: &[C], k: usize) -> Vec { + let d = dims[j]; + let bonds: usize = dims.iter().product(); + let stride: usize = dims[j + 1..].iter().product(); + let new_bonds = bonds / d * k; + let mut out = vec![C::ZERO; 2 * new_bonds]; + for s in 0..2 { + for flat in 0..bonds { + let x = t[s * bonds + flat]; + if x.re == 0.0 && x.im == 0.0 { + continue; + } + let pre = flat / (d * stride); + let i = (flat / stride) % d; + let post = flat % stride; + for n in 0..k { + let at = s * new_bonds + pre * k * stride + n * stride + post; + out[at] = out[at].add(x.mul(m[i * k + n])); + } + } + } + out +} + +/// SVD of a row-major `rows × cols` matrix: left singular vectors (each of +/// length `rows`), singular values (descending), right singular vectors (each +/// of length `cols`), with `A = Σ_k s_k · u_k · v_k†`. The one-sided Jacobi +/// sweeps run over the shorter dimension, to convergence: the singular vectors +/// must be orthonormal to rounding, or every later expectation inherits the error. +fn svd_rows(a: &[C], rows: usize, cols: usize) -> (Vec>, Vec, Vec>) { + if cols <= rows { + let columns: Vec> = (0..cols).map(|c| (0..rows).map(|r| a[r * cols + c]).collect()).collect(); + jacobi_svd_strict(columns, rows, cols) + } else { + // A† = V S U†: decompose the conjugate transpose, whose columns are + // the conjugated rows of A, and swap the roles. + let columns: Vec> = (0..rows).map(|r| (0..cols).map(|c| a[r * cols + c].conj()).collect()).collect(); + let (v, s, u) = jacobi_svd_strict(columns, cols, rows); + (u, s, v) + } +} + +fn pauli(p: Pauli) -> [[C; 2]; 2] { + let (z, o) = (C::ZERO, C::ONE); + match p { + Pauli::X => [[z, o], [o, z]], + Pauli::Y => [[z, C::new(0.0, -1.0)], [C::new(0.0, 1.0), z]], + Pauli::Z => [[o, z], [z, C::new(-1.0, 0.0)]], + } +} + +fn fixed1(op: &Op) -> [[C; 2]; 2] { + let (z, o) = (C::ZERO, C::ONE); + let (i, ni) = (C::new(0.0, 1.0), C::new(0.0, -1.0)); + let h = std::f64::consts::FRAC_1_SQRT_2; + let (hp, hm) = (C::new(h, 0.0), C::new(-h, 0.0)); + let (pp, pm) = (C::new(0.5, 0.5), C::new(0.5, -0.5)); + match op { + Op::H(_) => [[hp, hp], [hp, hm]], + Op::S(_) => [[o, z], [z, i]], + Op::Sdg(_) => [[o, z], [z, ni]], + Op::SX(_) => [[pp, pm], [pm, pp]], + Op::SXdg(_) => [[pm, pp], [pp, pm]], + Op::X(_) => pauli(Pauli::X), + Op::Y(_) => pauli(Pauli::Y), + _ => pauli(Pauli::Z), + } +} + +fn cx() -> [[C; 4]; 4] { + let (z, o) = (C::ZERO, C::ONE); + [[o, z, z, z], [z, o, z, z], [z, z, z, o], [z, z, o, z]] +} + +fn cz() -> [[C; 4]; 4] { + let (z, o) = (C::ZERO, C::ONE); + [[o, z, z, z], [z, o, z, z], [z, z, o, z], [z, z, z, C::new(-1.0, 0.0)]] +} + +fn swap() -> [[C; 4]; 4] { + let (z, o) = (C::ZERO, C::ONE); + [[o, z, z, z], [z, z, o, z], [z, o, z, z], [z, z, z, o]] +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::quantum_spd::{SpdConfig, expect_z, expectation, heavy_hex_127, kicked_ising}; + + /// Exact sparse Pauli dynamics: an independent oracle. + fn exact(c: &RotCircuit, q: u32, p: Pauli) -> f64 { + expectation(c, &[(1.0, vec![(q, p)])], &[], &SpdConfig { threshold: 0.0, max_weight: None }).unwrap().value + } + + struct Rng(u64); + impl Rng { + fn next(&mut self) -> u64 { + self.0 ^= self.0 << 13; + self.0 ^= self.0 >> 7; + self.0 ^= self.0 << 17; + self.0 + } + fn below(&mut self, n: u64) -> u64 { + self.next() % n + } + fn angle(&mut self) -> f64 { + (self.next() >> 11) as f64 / (1u64 << 53) as f64 * 6.0 - 3.0 + } + } + + /// A random circuit whose two-qubit gates all lie on `edges`. + fn random_on(rng: &mut Rng, n: u32, edges: &[(u32, u32)], len: usize) -> RotCircuit { + let p = |r: &mut Rng| [Pauli::X, Pauli::Y, Pauli::Z][r.below(3) as usize]; + let mut c = RotCircuit::new(n); + for _ in 0..len { + let q = rng.below(u64::from(n)) as u32; + let (a, b) = edges[rng.below(edges.len() as u64) as usize]; + let (a, b) = if rng.below(2) == 0 { (a, b) } else { (b, a) }; + let op = match rng.below(12) { + 0 => Op::H(q), + 1 => Op::S(q), + 2 => Op::SX(q), + 3 => Op::SXdg(q), + 4 => Op::Y(q), + 5 => Op::CX(a, b), + 6 => Op::CZ(a, b), + 7 => Op::Swap(a, b), + 8 => Op::QuarterRot { axis: vec![(a, p(rng)), (b, p(rng))], k: 1 + rng.below(3) as u8 }, + 9 => Op::Rot { axis: vec![(a, p(rng)), (b, p(rng))], theta: rng.angle() }, + _ => Op::Rot { axis: vec![(q, p(rng))], theta: rng.angle() }, + }; + c.push(op); + } + c + } + + #[test] + fn on_a_tree_with_no_truncation_it_is_exact() { + // A chain and a branching tree. With χ = 2^⌊n/2⌋ nothing is ever + // truncated, so the Vidal-gauge evolution and its environments are + // exact, and must agree with exact sparse Pauli dynamics to rounding. + let chain: Vec<(u32, u32)> = (0..7).map(|q| (q, q + 1)).collect(); + let tree = vec![(0, 1), (1, 2), (1, 3), (3, 4), (3, 5), (0, 6), (6, 7)]; + let mut rng = Rng(0x5eed); + for edges in [chain, tree] { + for trial in 0..25 { + let c = random_on(&mut rng, 8, &edges, 40); + let mut tns = Tns::new(8, &edges, 16).unwrap(); + tns.run(&c).unwrap(); + assert!(tns.discarded() < 1e-20, "nothing truncated: {}", tns.discarded()); + for q in 0..8 { + for p in [Pauli::X, Pauli::Y, Pauli::Z] { + let (got, want) = (tns.expect1(q, p).unwrap(), exact(&c, q, p)); + assert!((got - want).abs() < 1e-12, "trial {trial} q {q} {p:?}: {got} vs {want}"); + } + } + } + } + } + + #[test] + fn belief_propagation_is_exact_on_a_tree() { + let tree = vec![(0, 1), (1, 2), (1, 3), (3, 4), (3, 5), (0, 6), (6, 7)]; + let mut rng = Rng(0xbe1); + for trial in 0..15 { + let c = random_on(&mut rng, 8, &tree, 40); + let mut tns = Tns::new(8, &tree, 16).unwrap(); + tns.run(&c).unwrap(); + let env = tns.belief_propagation(200, 1e-14); + assert!(env.residual <= 1e-14, "trial {trial}: BP converged ({} sweeps, {:e})", env.sweeps, env.residual); + for q in 0..8 { + for p in [Pauli::X, Pauli::Y, Pauli::Z] { + let (got, want) = (tns.expect1_bp(&env, q, p).unwrap(), exact(&c, q, p)); + assert!((got - want).abs() < 1e-11, "trial {trial} q {q} {p:?}: {got} vs {want}"); + } + } + } + } + + #[test] + fn regauging_changes_no_expectation_and_makes_the_bond_weights_the_environment() { + let edges = heavy_hex_127(); + let mut tns = Tns::new(127, &edges, 8).unwrap(); + tns.run(&kicked_ising(127, &edges, 0.6, 4)).unwrap(); + let env = tns.belief_propagation(200, 1e-13); + let before: Vec = (0..127).map(|q| tns.expect1_bp(&env, q, Pauli::Z).unwrap()).collect(); + tns.regauge(&env); + let env2 = tns.belief_propagation(200, 1e-13); + for q in 0..127 { + let (after_bp, after_simple) = (tns.expect1_bp(&env2, q, Pauli::Z).unwrap(), tns.expect1(q, Pauli::Z).unwrap()); + assert!((after_bp - before[q as usize]).abs() < 1e-9, "q {q}: the state changed: {after_bp} vs {}", before[q as usize]); + assert!((after_simple - after_bp).abs() < 1e-9, "q {q}: bond weights {after_simple} vs BP {after_bp}"); + } + } + + #[test] + fn gauged_evolution_is_exact_on_a_tree() { + let tree = vec![(0, 1), (1, 2), (1, 3), (3, 4), (3, 5), (0, 6), (6, 7)]; + let mut rng = Rng(0x9a); + for trial in 0..10 { + let c = random_on(&mut rng, 8, &tree, 40); + let mut tns = Tns::new(8, &tree, 16).unwrap(); + tns.run_gauged(&c, 100, 1e-14).unwrap(); + for q in 0..8 { + let (got, want) = (tns.expect1(q, Pauli::Z).unwrap(), exact(&c, q, Pauli::Z)); + assert!((got - want).abs() < 1e-10, "trial {trial} q {q}: {got} vs {want}"); + } + } + } + + #[test] + fn belief_propagation_is_bit_identical_however_it_is_parallelised() { + let edges = heavy_hex_127(); + let mut tns = Tns::new(127, &edges, 4).unwrap(); + tns.run(&kicked_ising(127, &edges, 0.5, 3)).unwrap(); + let a = tns.belief_propagation(30, 1e-12); + let b = tns.belief_propagation(30, 1e-12); + assert_eq!((a.sweeps, a.residual.to_bits()), (b.sweeps, b.residual.to_bits())); + assert_eq!(tns.expect1_bp(&a, 62, Pauli::Z).unwrap().to_bits(), tns.expect1_bp(&b, 62, Pauli::Z).unwrap().to_bits()); + } + + #[test] + fn truncation_is_reported_and_bounded_by_chi() { + let chain: Vec<(u32, u32)> = (0..9).map(|q| (q, q + 1)).collect(); + let mut rng = Rng(77); + let c = random_on(&mut rng, 10, &chain, 200); + let mut tns = Tns::new(10, &chain, 4).unwrap(); + tns.run(&c).unwrap(); + assert!(tns.max_bond() <= 4); + assert!(tns.discarded() > 0.0, "a deep random circuit at χ = 4 must truncate"); + } + + #[test] + fn a_gate_off_the_graph_is_refused() { + let mut tns = Tns::new(3, &[(0, 1), (1, 2)], 4).unwrap(); + assert_eq!(tns.apply(&Op::CX(0, 2)), Err(TnsError::NotAnEdge(0, 2))); + assert_eq!( + tns.apply(&Op::Rot { axis: vec![(0, Pauli::Z), (1, Pauli::Z), (2, Pauli::Z)], theta: 0.1 }), + Err(TnsError::TooWide(3)) + ); + } + + #[test] + fn results_are_bit_identical_from_run_to_run() { + let c = kicked_ising(127, &heavy_hex_127(), 0.6, 3); + let run = || { + let mut t = Tns::new(127, &heavy_hex_127(), 8).unwrap(); + t.run(&c).unwrap(); + (t.expect1(62, Pauli::Z).unwrap().to_bits(), t.discarded().to_bits()) + }; + assert_eq!(run(), run()); + } + + #[test] + fn a_layer_computed_in_parallel_gives_the_bits_of_one_gate_at_a_time() { + let edges = heavy_hex_127(); + let c = kicked_ising(127, &edges, 0.7, 3); + let mut batched = Tns::new(127, &edges, 8).unwrap(); + batched.run(&c).unwrap(); + let mut serial = Tns::new(127, &edges, 8).unwrap(); + for op in &c.ops { + serial.apply(op).unwrap(); + } + assert_eq!(batched.discarded().to_bits(), serial.discarded().to_bits()); + for q in 0..127 { + assert_eq!(batched.expect1(q, Pauli::Z).unwrap().to_bits(), serial.expect1(q, Pauli::Z).unwrap().to_bits()); + } + } + + #[test] + fn the_heavy_hex_magnetization_tracks_sparse_pauli_dynamics() { + // On the 127-qubit lattice, with its loops, belief propagation is an + // approximation. Five steps at small angle: the two independent + // methods agree closely. + let edges = heavy_hex_127(); + let c = kicked_ising(127, &edges, std::f64::consts::PI / 16.0, 5); + let mut tns = Tns::new(127, &edges, 8).unwrap(); + tns.run(&c).unwrap(); + let spd = expect_z(&c, 62, &SpdConfig { threshold: 1e-7, max_weight: None }).unwrap().value; + let bp = tns.expect1(62, Pauli::Z).unwrap(); + assert!((bp - spd).abs() < 2e-3, "BP {bp} vs SPD {spd}"); + } +} diff --git a/wai/wai-quantum/src/quantum_mps.rs b/wai/wai-quantum/src/quantum_mps.rs index fcf815120..d32e7920b 100644 --- a/wai/wai-quantum/src/quantum_mps.rs +++ b/wai/wai-quantum/src/quantum_mps.rs @@ -23,140 +23,8 @@ use crate::quantum::{BaseGate, Circuit, Gate, StateVector}; -// =========================================================================== -// Complex f64 (deterministic: only + - * / sqrt, no FMA/transcendentals) -// =========================================================================== - -#[derive(Clone, Copy, Debug)] -pub struct C { - pub re: f64, - pub im: f64, -} -impl C { - const ZERO: C = C { re: 0.0, im: 0.0 }; - const ONE: C = C { re: 1.0, im: 0.0 }; - #[inline] - fn new(re: f64, im: f64) -> C { - C { re, im } - } - #[inline] - fn add(self, o: C) -> C { - C { re: self.re + o.re, im: self.im + o.im } - } - #[inline] - fn sub(self, o: C) -> C { - C { re: self.re - o.re, im: self.im - o.im } - } - #[inline] - fn mul(self, o: C) -> C { - C { re: self.re * o.re - self.im * o.im, im: self.re * o.im + self.im * o.re } - } - #[inline] - fn conj(self) -> C { - C { re: self.re, im: -self.im } - } - #[inline] - fn scale(self, s: f64) -> C { - C { re: self.re * s, im: self.im * s } - } - #[inline] - fn norm2(self) -> f64 { - self.re * self.re + self.im * self.im - } -} - -// =========================================================================== -// Complex one-sided Jacobi SVD (M = U · diag(S) · V†) -// =========================================================================== - -/// SVD of an `rows × cols` matrix given by its **columns** (`m[c]` is column `c`, a -/// length-`rows` vector). Returns `(u_cols, s, v_cols)` where each `u_cols[k]` is a -/// length-`rows` left singular vector, `s[k]` its singular value, `v_cols[k]` a -/// length-`cols` right singular vector — sorted by `s` descending. Deterministic: -/// fixed sweep order, fixed sweep cap, only `+ − × ÷ √`. -fn jacobi_svd(mut m: Vec>, rows: usize, cols: usize) -> (Vec>, Vec, Vec>) { - // V starts as identity (cols × cols), stored by columns. - let mut v: Vec> = (0..cols) - .map(|j| (0..cols).map(|i| if i == j { C::ONE } else { C::ZERO }).collect()) - .collect(); - let eps = 1e-14; - for _sweep in 0..60 { - let mut off = 0.0f64; - for p in 0..cols { - for q in (p + 1)..cols { - // gram entries of columns p,q - let mut alpha = 0.0; // - let mut beta = 0.0; // - let mut gamma = C::ZERO; // = sum conj(c_p)·c_q - for k in 0..rows { - alpha += m[p][k].norm2(); - beta += m[q][k].norm2(); - gamma = gamma.add(m[p][k].conj().mul(m[q][k])); - } - let g2 = gamma.norm2(); - off += g2; - if g2 <= eps * alpha * beta || g2 == 0.0 { - continue; - } - let gabs = g2.sqrt(); - // phase-rotate column q so becomes real +gabs - let ph = C::new(gamma.re / gabs, gamma.im / gabs); // unit - let cph = ph.conj(); - for k in 0..rows { - m[q][k] = m[q][k].mul(cph); - } - for k in 0..cols { - v[q][k] = v[q][k].mul(cph); - } - // real Jacobi rotation with a=alpha,b=beta,g=gabs - let tau = (beta - alpha) / (2.0 * gabs); - let t = if tau >= 0.0 { - 1.0 / (tau + (1.0 + tau * tau).sqrt()) - } else { - -1.0 / (-tau + (1.0 + tau * tau).sqrt()) - }; - let cs = 1.0 / (1.0 + t * t).sqrt(); - let sn = t * cs; - for k in 0..rows { - let x = m[p][k]; - let y = m[q][k]; - m[p][k] = x.scale(cs).sub(y.scale(sn)); - m[q][k] = x.scale(sn).add(y.scale(cs)); - } - for k in 0..cols { - let x = v[p][k]; - let y = v[q][k]; - v[p][k] = x.scale(cs).sub(y.scale(sn)); - v[q][k] = x.scale(sn).add(y.scale(cs)); - } - } - } - if off <= eps { - break; - } - } - // singular values = column norms; left vectors = normalized columns - let mut s = vec![0.0; cols]; - for j in 0..cols { - s[j] = (0..rows).map(|k| m[j][k].norm2()).sum::().sqrt(); - } - let mut u: Vec> = m; - for j in 0..cols { - if s[j] > 1e-300 { - let inv = 1.0 / s[j]; - for k in 0..rows { - u[j][k] = u[j][k].scale(inv); - } - } - } - // sort by singular value descending (deterministic) - let mut order: Vec = (0..cols).collect(); - order.sort_by(|&a, &b| s[b].partial_cmp(&s[a]).unwrap_or(std::cmp::Ordering::Equal)); - let s2: Vec = order.iter().map(|&i| s[i]).collect(); - let u2: Vec> = order.iter().map(|&i| u[i].clone()).collect(); - let v2: Vec> = order.iter().map(|&i| v[i].clone()).collect(); - (u2, s2, v2) -} +pub use crate::linalg::C; +use crate::linalg::jacobi_svd; // =========================================================================== // MPS diff --git a/wai/wai-quantum/src/quantum_spd.rs b/wai/wai-quantum/src/quantum_spd.rs index 49b844e79..695034b60 100644 --- a/wai/wai-quantum/src/quantum_spd.rs +++ b/wai/wai-quantum/src/quantum_spd.rs @@ -687,18 +687,100 @@ pub fn heavy_hex_127() -> Vec<(u32, u32)> { edges } +/// The edges as layers of vertex-disjoint edges: an edge colouring. A graph +/// whose maximum degree is `Δ` and which is bipartite, such as heavy-hex, gets +/// exactly `Δ` layers (König's theorem): each edge takes a colour free at both +/// ends, and when none is, an alternating path is recoloured to free one. A +/// graph with an odd cycle may need a layer more, which it is given. Within a +/// layer the edges are in their input order, and the layers in colour order. +pub fn edge_layers(n: u32, edges: &[(u32, u32)]) -> Vec> { + let n = n as usize; + let mut deg = vec![0usize; n]; + for &(a, b) in edges { + deg[a as usize] += 1; + deg[b as usize] += 1; + } + let mut ncol = deg.iter().copied().max().unwrap_or(0).max(1); + // at[v][c]: the edge of colour c at vertex v, if any + let mut at: Vec>> = vec![vec![None; ncol]; n]; + let mut colour = vec![usize::MAX; edges.len()]; + let other = |f: usize, x: usize| { + let (a, b) = edges[f]; + if a as usize == x { b as usize } else { a as usize } + }; + for (e, &(u, v)) in edges.iter().enumerate() { + let (u, v) = (u as usize, v as usize); + let free_u = (0..ncol).find(|&c| at[u][c].is_none()); + let free_v = (0..ncol).find(|&c| at[v][c].is_none()); + let chosen = match (free_u, free_v) { + (Some(a), Some(_)) if at[v][a].is_none() => Some(a), + (Some(a), Some(b)) => { + // Walk the a/b alternating path from v; flipping it frees a at v. + let mut path = Vec::new(); + let (mut x, mut c) = (v, a); + let mut reaches_u = false; + while let Some(f) = at[x][c] { + path.push(f); + x = other(f, x); + if x == u { + reaches_u = true; + break; + } + c = if c == a { b } else { a }; + } + if reaches_u { + None + } else { + for &f in &path { + let (p, q) = (edges[f].0 as usize, edges[f].1 as usize); + at[p][colour[f]] = None; + at[q][colour[f]] = None; + } + for &f in &path { + let (p, q) = (edges[f].0 as usize, edges[f].1 as usize); + colour[f] = if colour[f] == a { b } else { a }; + at[p][colour[f]] = Some(f); + at[q][colour[f]] = Some(f); + } + Some(a) + } + } + _ => None, + }; + let c = chosen.unwrap_or_else(|| { + ncol += 1; + for row in at.iter_mut() { + row.push(None); + } + ncol - 1 + }); + colour[e] = c; + at[u][c] = Some(e); + at[v][c] = Some(e); + } + (0..ncol) + .map(|c| edges.iter().zip(&colour).filter(|&(_, &k)| k == c).map(|(&e, _)| e).collect::>()) + .filter(|layer| !layer.is_empty()) + .collect() +} + /// The kicked-Ising circuit of the heavy-hex utility experiment: `steps` /// Trotter steps, each `RX(θ_h)` on every qubit and then `RZZ(−π/2)` on every -/// edge (applied exactly, as quarter turns). +/// edge, layer by layer as [`edge_layers`] partitions them: three layers of +/// vertex-disjoint edges on heavy-hex, as on the device. The `RZZ` gates are applied +/// exactly, as quarter turns. pub fn kicked_ising(n: u32, edges: &[(u32, u32)], theta_h: f64, steps: u32) -> RotCircuit { + let layers = edge_layers(n, edges); let mut c = RotCircuit::new(n); for _ in 0..steps { for q in 0..n { c.rx(q, theta_h); } - for &(a, b) in edges { - // RZZ(−π/2) = exp(+iπ/4·ZZ) = exp(−i(3π/2)ZZ/2) up to a global phase. - c.push(Op::QuarterRot { axis: vec![(a, Pauli::Z), (b, Pauli::Z)], k: 3 }); + for layer in &layers { + for &(a, b) in layer { + // RZZ(−π/2) = exp(+iπ/4·ZZ) = exp(−i(3π/2)ZZ/2) up to a global phase. + c.push(Op::QuarterRot { axis: vec![(a, Pauli::Z), (b, Pauli::Z)], k: 3 }); + } } } c @@ -1175,6 +1257,30 @@ mod tests { } } + #[test] + fn edge_layers_are_proper_colourings_with_as_few_layers_as_the_degree() { + let check = |n: u32, edges: &[(u32, u32)], want: usize| { + let layers = edge_layers(n, edges); + assert_eq!(layers.len(), want); + let mut all: Vec<(u32, u32)> = layers.iter().flatten().copied().collect(); + all.sort_unstable(); + let mut sorted = edges.to_vec(); + sorted.sort_unstable(); + assert_eq!(all, sorted, "every edge exactly once"); + for layer in &layers { + let mut seen = std::collections::BTreeSet::new(); + for &(a, b) in layer { + assert!(seen.insert(a) && seen.insert(b), "vertex-disjoint"); + } + } + }; + check(127, &heavy_hex_127(), 3); + let chain: Vec<(u32, u32)> = (0..9).map(|q| (q, q + 1)).collect(); + check(10, &chain, 2); + // A triangle is not bipartite: three edges, degree two, three layers. + check(3, &[(0, 1), (1, 2), (0, 2)], 3); + } + #[test] fn the_kicked_ising_clifford_points_are_exact() { // θ_h = 0: no rotation, every ⟨Z_q⟩ = 1. θ_h = π/2 written as quarter diff --git a/wai/wai-quantum/tests/bit_identity.rs b/wai/wai-quantum/tests/bit_identity.rs index b93dc549d..20d2f74c4 100644 --- a/wai/wai-quantum/tests/bit_identity.rs +++ b/wai/wai-quantum/tests/bit_identity.rs @@ -13,7 +13,8 @@ feature = "quantum_kernel", feature = "quantum_mps", feature = "quantum_pauli", - feature = "quantum_spd" + feature = "quantum_spd", + feature = "quantum_bptn" ))] use wai_quantum::*; @@ -107,3 +108,16 @@ fn sparse_pauli_dynamics() { let r = quantum_spd::expect_z(&circuit, 7, &cfg).unwrap(); assert_eq!(digest(&[r.value, r.truncation_bound]), 0xfc44_f4b6_0e70_68b2); } + +#[test] +fn tensor_network_evolution() { + // Every qubit's ⟨Z⟩ and the discarded weight after three steps on the + // 127-qubit lattice at χ = 8: QR, the converged SVD and the parallel + // layers all feed these bits. + let edges = quantum_spd::heavy_hex_127(); + let mut t = quantum_bptn::Tns::new(127, &edges, 8).unwrap(); + t.run(&quantum_spd::kicked_ising(127, &edges, 0.6, 3)).unwrap(); + let mut v: Vec = (0..127).map(|q| t.expect1(q, quantum_spd::Pauli::Z).unwrap()).collect(); + v.push(t.discarded()); + assert_eq!(digest(&v), 0x0e51_33ff_7987_aa26); +} diff --git a/wai/wai-quantum/tests/utility_reference.rs b/wai/wai-quantum/tests/utility_reference.rs index d5c5d87a2..69c593f77 100644 --- a/wai/wai-quantum/tests/utility_reference.rs +++ b/wai/wai-quantum/tests/utility_reference.rs @@ -12,6 +12,22 @@ use wai_quantum::quantum_spd::{Pauli, SpdConfig, expect_z, expectation, heavy_hex_127, kicked_ising}; +/// The same magnetization by belief-propagation tensor networks at χ = 16, an +/// independent method: the light cone of five steps fits, and the values are +/// exact to rounding of the published six decimals. +#[cfg(feature = "quantum_bptn")] +#[test] +fn five_step_magnetization_by_tensor_networks_matches_the_exact_values() { + use wai_quantum::quantum_bptn::Tns; + let edges = heavy_hex_127(); + for (k, exact) in [(2, 0.979836), (4, 0.904843)] { + let mut t = Tns::new(127, &edges, 16).unwrap(); + t.run(&kicked_ising(127, &edges, theta(k), 5)).unwrap(); + let mz: f64 = (0..127).map(|q| t.expect1(q, Pauli::Z).unwrap()).sum::() / 127.0; + assert!((mz - exact).abs() < 1e-6, "k = {k}: Mz = {mz} against {exact}"); + } +} + fn theta(k: u32) -> f64 { f64::from(k) * std::f64::consts::PI / 32.0 } diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index 8f783e1d6..f575aca24 100644 --- a/wai/wai-rs/Cargo.lock +++ b/wai/wai-rs/Cargo.lock @@ -3900,7 +3900,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.22" +version = "0.3.23" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-rs/Cargo.toml b/wai/wai-rs/Cargo.toml index 552ccadff..f7496f134 100644 --- a/wai/wai-rs/Cargo.toml +++ b/wai/wai-rs/Cargo.toml @@ -339,6 +339,8 @@ quantum_pauli = ["quantum", "wai-quantum/quantum_pauli"] # and its signed receipts. quantum_spd = ["quantum", "wai-quantum/quantum_spd"] quantum_spd_receipt = ["quantum_spd", "wai-quantum/quantum_spd_receipt"] +# Tensor-network states on the hardware graph (simple update, belief propagation). +quantum_bptn = ["quantum_spd", "wai-quantum/quantum_bptn"] quantum_qir = ["quantum_toolchain", "wai-quantum/quantum_qir"] # Variational quantum machine learning (extensions/quantum-sim § VML, # `wai.quantum.vml`): the interference-based ML the thermodynamic-computing pivot @@ -465,6 +467,7 @@ required-features = [ "quantum_comm", "quantum_qudit", "quantum_spd_receipt", + "quantum_bptn", ] [[bin]] diff --git a/wai/wai-rs/src/bin/wai_quantum.rs b/wai/wai-rs/src/bin/wai_quantum.rs index a48f7cb95..1036a78b9 100644 --- a/wai/wai-rs/src/bin/wai_quantum.rs +++ b/wai/wai-rs/src/bin/wai_quantum.rs @@ -32,6 +32,7 @@ //! //! ```text //! wai-quantum utility --theta 0.2 --steps 20 # the 127-qubit heavy-hex experiment +//! wai-quantum utility --method bptn --chi 32 # the same, by tensor networks //! wai-quantum utility --receipt | wai-quantum verify //! ``` //! @@ -125,7 +126,8 @@ fn usage() -> ! { wai-quantum chsh [--json]\n \ wai-quantum qutrit [-n QUTRITS] [--json]\n\n\ UTILITY SCALE:\n \ - wai-quantum utility [--theta T] [--steps N] [--threshold D] [-q QUBIT] [--json] [--receipt]\n\n\ + wai-quantum utility [--theta T] [--steps N] [--threshold D] [-q QUBIT] [--json] [--receipt]\n \ + wai-quantum utility --method bptn [--chi X] [--theta T] [--steps N] [-q QUBIT] [--json]\n\n\ Channels: identity, bitflip, dephasing, depolarizing, damping.\n\ Run `wai-quantum list` for the algorithm catalog." ); @@ -324,6 +326,35 @@ fn main() { let steps: u32 = num(args, "--steps", "", 20); let threshold: f64 = num(args, "--threshold", "", 1e-3); let qubit: u32 = num(args, "--qubit", "-q", 62); + // The same circuit by belief-propagation tensor networks: an + // independent method, whose errors are unrelated to truncating a + // Pauli sum. + if opt(args, "--method", "").as_deref() == Some("bptn") { + if flag(args, "--receipt") { + die("--receipt seals sparse Pauli dynamics runs; tensor-network runs have no receipt yet"); + } + let chi: usize = num(args, "--chi", "", 16); + let edges = heavy_hex_127(); + let mut tns = wai::quantum_bptn::Tns::new(127, &edges, chi).unwrap_or_else(|e| die(&e.to_string())); + tns.run(&kicked_ising(127, &edges, theta, steps)).unwrap_or_else(|e| die(&e.to_string())); + let value = tns.expect1(qubit, Pauli::Z).unwrap_or_else(|e| die(&e.to_string())); + if flag(args, "--json") { + println!( + "{{\"chi\":{chi},\"discarded\":{},\"max_bond\":{},\"method\":\"bptn\",\"qubit\":{qubit},\ + \"steps\":{steps},\"theta_h\":{theta},\"value\":{value}}}", + tns.discarded(), + tns.max_bond() + ); + } else { + println!(" {value:.6}"); + println!("discarded {:.3e} (squared singular-value weight truncated)", tns.discarded()); + println!("max bond {}", tns.max_bond()); + } + return; + } + if opt(args, "--method", "").as_deref().is_some_and(|m| m != "spd") { + die("--method is spd (the default) or bptn"); + } let job = SpdJob { circuit: kicked_ising(127, &heavy_hex_127(), theta, steps), observable: vec![(1.0, vec![(qubit, Pauli::Z)])], diff --git a/wai/wai-rs/src/lib.rs b/wai/wai-rs/src/lib.rs index ff3e83f15..d1e8661f7 100644 --- a/wai/wai-rs/src/lib.rs +++ b/wai/wai-rs/src/lib.rs @@ -587,6 +587,9 @@ pub use wai_quantum::quantum_spd; #[cfg(feature = "quantum_spd_receipt")] pub use wai_quantum::quantum_spd_receipt; + +#[cfg(feature = "quantum_bptn")] +pub use wai_quantum::quantum_bptn; #[cfg(feature = "quantum_qir")] pub use wai_quantum::quantum_qir; diff --git a/wai/wai-web/Cargo.lock b/wai/wai-web/Cargo.lock index c75264d66..4720f02a9 100644 --- a/wai/wai-web/Cargo.lock +++ b/wai/wai-web/Cargo.lock @@ -734,7 +734,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.22" +version = "0.3.23" dependencies = [ "blake3", "ed25519-dalek", -- 2.54.0 (Apple Git-157)