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

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

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"]

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
 

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<usize> = 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::<f64>() / 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()
+            );
+        }
+    }
+}

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;
 

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<Vec<C>>, rows: usize, cols: usize) -> (Vec<Vec<C>>, Vec<f64>, Vec<Vec<C>>) {
+    // V starts as identity (cols × cols), stored by columns.
+    let mut v: Vec<Vec<C>> = (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; // <c_p,c_p>
+                let mut beta = 0.0; // <c_q,c_q>
+                let mut gamma = C::ZERO; // <c_p,c_q> = 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 <c_p,c_q> 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::<f64>().sqrt();
+    }
+    let mut u: Vec<Vec<C>> = 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<usize> = (0..cols).collect();
+    order.sort_by(|&a, &b| s[b].partial_cmp(&s[a]).unwrap_or(std::cmp::Ordering::Equal));
+    let s2: Vec<f64> = order.iter().map(|&i| s[i]).collect();
+    let u2: Vec<Vec<C>> = order.iter().map(|&i| u[i].clone()).collect();
+    let v2: Vec<Vec<C>> = 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<Vec<C>>, rows: usize, cols: usize) -> (Vec<Vec<C>>, Vec<f64>, Vec<Vec<C>>) {
+    // V starts as identity (cols × cols), stored by columns.
+    let mut v: Vec<Vec<C>> = (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::<f64>();
+    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; // <c_p,c_p>
+                let mut beta = 0.0; // <c_q,c_q>
+                let mut gamma = C::ZERO; // <c_p,c_q> = 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 <c_p,c_q> 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::<f64>().sqrt();
+    }
+    let mut u: Vec<Vec<C>> = 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<usize> = (0..cols).collect();
+    order.sort_by(|&a, &b| s[b].partial_cmp(&s[a]).unwrap_or(std::cmp::Ordering::Equal));
+    let s2: Vec<f64> = order.iter().map(|&i| s[i]).collect();
+    let u2: Vec<Vec<C>> = order.iter().map(|&i| u[i].clone()).collect();
+    let v2: Vec<Vec<C>> = 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<C>, Vec<C>, usize) {
+    let k = m.min(n);
+    let mut r = a.to_vec();
+    let mut vs: Vec<Vec<C>> = Vec::with_capacity(k);
+    for j in 0..k {
+        let norm = (j..m).map(|i| r[i * n + j].norm2()).sum::<f64>().sqrt();
+        let mut v: Vec<C> = (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::<f64>().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<C> = (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<C> = (0..rows * rank).map(|_| C::new(next(), next())).collect();
+            let c: Vec<C> = (0..cols * rank).map(|_| C::new(next(), next())).collect();
+            let a: Vec<C> = (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<Vec<C>> = (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");
+                    }
+                }
+            }
+        }
+    }
+}

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<Vec<C>>, rows: usize, cols: usize) -> (Vec<Vec<C>>, Vec<f64>, Vec<Vec<C>>) {
-    // V starts as identity (cols × cols), stored by columns.
-    let mut v: Vec<Vec<C>> = (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; // <c_p,c_p>
-                let mut beta = 0.0; // <c_q,c_q>
-                let mut gamma = C::ZERO; // <c_p,c_q> = 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 <c_p,c_q> 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::<f64>().sqrt();
-    }
-    let mut u: Vec<Vec<C>> = 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<usize> = (0..cols).collect();
-    order.sort_by(|&a, &b| s[b].partial_cmp(&s[a]).unwrap_or(std::cmp::Ordering::Equal));
-    let s2: Vec<f64> = order.iter().map(|&i| s[i]).collect();
-    let u2: Vec<Vec<C>> = order.iter().map(|&i| u[i].clone()).collect();
-    let v2: Vec<Vec<C>> = order.iter().map(|&i| v[i].clone()).collect();
-    (u2, s2, v2)
-}
+pub use crate::linalg::C;
+use crate::linalg::jacobi_svd;
 
 // ===========================================================================
 // MPS

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<Vec<(u32, u32)>> {
+    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<Option<usize>>> = 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::<Vec<_>>())
+        .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

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<f64> = (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);
+}

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::<f64>() / 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
 }

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

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

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!("<Z_{qubit}>     {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)])],

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;
 

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