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