diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index 275c89018..e719b744d 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.34" +version = "0.3.35" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index 613315311..9368b760c 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.34" +version = "0.3.35" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index 0f863de09..e4419e5aa 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "wai-quantum" -version = "0.3.34" +version = "0.3.35" edition = "2024" rust-version = "1.97" license = "Apache-2.0" @@ -82,6 +82,9 @@ quantum_noisy = ["quantum_sv"] # Ground states of spin chains by two-site DMRG, held to exact diagonalisation # and to the free-fermion solution of the transverse-field Ising chain. quantum_dmrg = [] +# Time evolution of spin chains by the time-dependent variational principle, +# held to exact evolution and to the free-fermion solution of a quench. +quantum_tdvp = ["quantum_dmrg"] # 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"] @@ -121,7 +124,7 @@ quantum_phasor_meter = ["_signing"] # Everything. full = [ - "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "quantum_sv", "quantum_arch", "quantum_ldpc", "quantum_qio", "quantum_uf", "quantum_noisy", "quantum_dmrg", "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_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "quantum_sv", "quantum_arch", "quantum_ldpc", "quantum_qio", "quantum_uf", "quantum_noisy", "quantum_dmrg", "quantum_tdvp", "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", diff --git a/wai/wai-quantum/README.md b/wai/wai-quantum/README.md index bb1adfdd8..5fa68ca94 100644 --- a/wai/wai-quantum/README.md +++ b/wai/wai-quantum/README.md @@ -67,6 +67,7 @@ never creates information. | `quantum_qio` | **quantum-inspired optimisation**: Ising and Max-Cut problems (G-set reader, spin-glass generator), simulated bifurcation in its ballistic, discrete and heated forms (batched eight trials per pass, bit-identical to one at a time), simulated annealing, single-flip descent, an exact Gray-code referee up to 34 spins, and step-to-solution. Reaches the best-known cut on seven G-set instances | | `quantum_uf` | **union-find decoding in near-linear time** over the same matching graph: event-driven cluster growth on integer weights (a radix heap orders the fill events), merges, and peeling. Measured against exact matching on the same shots: within 0.92–1.18× the failures, up to 77× faster at `d = 11` | | `quantum_dmrg` | **ground states of spin chains by DMRG**: matrix-product operators for nearest-neighbour chains (XXZ Heisenberg in a field, transverse-field Ising), two-site sweeps with Lanczos and density-matrix truncation, entanglement entropies, a symmetric eigensolver, and two exact referees: full diagonalisation up to 22 sites, and the free-fermion solution of the Ising chain at any length. The critical Ising chain is reproduced to about `10⁻¹¹` per site at 200 sites | +| `quantum_tdvp` | **time evolution of spin chains by TDVP**: two-site sweeps forward and back, Lanczos exponentials that halve a step they cannot converge, bonds split by singular value decomposition, and every basis completed with null vectors so a product state's small bonds do not hold the evolution back. Measures site and bond observables, entanglement at every bond, energy and the Loschmidt echo. Held to two exact referees: Krylov evolution of all amplitudes up to 20 sites, and the free-fermion solution of an Ising quench at any length. On a 100-site critical quench at bond 64 it stays within `10⁻⁹` of the exact answer until the entanglement fills the bond, and past that its error follows the discarded weight it reports | | `quantum_noisy` | **noisy circuits**: quantum trajectories on the state vector (depolarizing, dephasing, amplitude damping, readout), an exact density-matrix referee up to 12 qubits that reuses the same kernels, linear cross-entropy fidelity (sampled, exact and the circuit's own), and the digital error model. Measures where that model is off: below the true fidelity by a margin that grows with the error rate, not with depth | | `quantum_ldpc` | **qLDPC memory at circuit level**: bivariate-bicycle codes from their two polynomials (`k` by rank, logical bases), the depth-8 syndrome cycle with circuit-level noise, its detector error model, and BP-OSD decoding (min-sum with adaptive scaling, ordered statistics with a combination sweep). The error model has the same columns as the source's own decoding matrices, and on identical syndromes the decoder makes the source's decisions shot for shot | | `quantum_arch` | **architecture-level estimation**: yoked surface-code storage (row and grid parity checks over patches, hot and cold), magic-state cultivation read off the published simulation curves, CCZ factories, and a clock set by whichever is slowest of lattice surgery, the magic-state supply and the control loop. Assembles whole-machine estimates; its storage search matches the source's footprint tool layout for layout | @@ -197,6 +198,28 @@ At 200 sites the critical Ising chain comes out within about `10⁻¹¹` per sit of its exact energy. Runs repeat bit for bit at any thread count, natively and as WebAssembly. +## Time evolution + +```sh +wai-quantum tdvp -n 100 --bond 64 --dt 0.1 --t 4 # critical Ising quench, against free fermions +wai-quantum tdvp --chain heisenberg -n 100 --bond 64 --t 3 # the Néel state melting, entanglement included +``` + +A quantum simulation of dynamics is checkable when a classical method can +follow it. `quantum_tdvp` evolves a matrix-product state by the +time-dependent variational principle and reports what it gave up: the +discarded weight, the energy drift and the entanglement at every bond. It is +held to exact evolution of every amplitude on short chains, and to the +free-fermion solution of the Ising quench at any length, including the +Loschmidt echo. + +On the 100-site critical Ising quench at bond 64, every site's `σˣ` stays +within `10⁻⁹` of the exact answer through `t = 3`, with the Loschmidt rate +agreeing to seven digits. When the entanglement outgrows the bond, the error +rises with the discarded weight, and the run says so. A chain small enough to +fit the bond evolves exactly, to rounding, at any time step. Runs repeat bit +for bit at any thread count, natively and as WebAssembly. + ## Noisy circuits ```sh diff --git a/wai/wai-quantum/examples/tdvp_quench.rs b/wai/wai-quantum/examples/tdvp_quench.rs new file mode 100644 index 000000000..89aa30542 --- /dev/null +++ b/wai/wai-quantum/examples/tdvp_quench.rs @@ -0,0 +1,54 @@ +//! A quench of the transverse-field Ising chain by TDVP, held at every +//! report to the exact free-fermion solution: the chain starts with every +//! σˣ = +1 and evolves under `−J Σ σᶻσᶻ − h Σ σˣ`. +//! +//! ```text +//! cargo run --release --example tdvp_quench --features quantum_tdvp -- N BOND DT T [h] [cutoff] [threads] +//! ``` +use wai_quantum::quantum_dmrg::{Chain, Mps, Op, Scalar, C}; +use wai_quantum::quantum_tdvp::{ising_quench, Evolution, TdvpConfig}; + +fn main() { + let a: Vec = std::env::args().collect(); + let n: usize = a[1].parse().unwrap(); + let bond: usize = a[2].parse().unwrap(); + let dt: f64 = a[3].parse().unwrap(); + let t_end: f64 = a[4].parse().unwrap(); + let h: f64 = a.get(5).map_or(1.0, |x| x.parse().unwrap()); + let cutoff: f64 = a.get(6).map_or(1e-20, |x| x.parse().unwrap()); + let threads: usize = a.get(7).map_or(TdvpConfig::default().threads, |x| x.parse().unwrap()); + let chain = Chain::ising(n, 1.0, h); + let plus = Mps::product(&vec![[C { re: 1.0, im: 0.0 }, C { re: 1.0, im: 0.0 }]; n]); + let cfg = TdvpConfig { dt, max_bond: bond, cutoff, threads, ..TdvpConfig::default() }; + let mut ev = Evolution::new(&chain, &plus, cfg); + let e0 = ev.energy(); + let steps = (t_end / dt).round() as usize; + let every = ((0.5 / dt).round() as usize).max(1); + let t0 = std::time::Instant::now(); + let mut discarded = 0.0; + println!("ising n={n} h={h} bond={bond} dt={dt} cutoff={cutoff:e}"); + println!(" t max|Δσˣ| max|Δσᶻσᶻ| σˣ(mid) exact rate exact rate S(mid) bond |ΔE| discarded seconds"); + for s in 1..=steps { + discarded += ev.step().discarded; + if s % every == 0 || s == steps { + let t = ev.time(); + let ff = ising_quench(n, 1.0, h, t); + let x = ev.expect(Op::x()); + let zz = ev.bond_expect(Op::z(), Op::z()); + let dx = x.iter().zip(&ff.x).map(|(p, q)| (p - q).abs()).fold(0.0, f64::max); + let dzz = zz.iter().zip(&ff.zz).map(|(p, q)| (p - q).abs()).fold(0.0, f64::max); + let echo = ev.overlap(&plus).norm2(); + let rate = -echo.ln() / n as f64; + let exact_rate = -ff.echo.ln() / n as f64; + let s_mid = ev.entropies()[n / 2 - 1]; + println!( + "{t:5.2} {dx:.2e} {dzz:.2e} {:+.10} {:+.10} {rate:.8} {exact_rate:.8} {s_mid:.4} {:4} {:.1e} {discarded:.1e} {:.1}", + x[n / 2], + ff.x[n / 2], + ev.max_bond(), + (ev.energy() - e0).abs(), + t0.elapsed().as_secs_f64() + ); + } + } +} diff --git a/wai/wai-quantum/src/lib.rs b/wai/wai-quantum/src/lib.rs index fcfeb6d2d..af3bb2bf9 100644 --- a/wai/wai-quantum/src/lib.rs +++ b/wai/wai-quantum/src/lib.rs @@ -17,7 +17,7 @@ #[cfg(feature = "_signing")] mod merkle; -#[cfg(any(feature = "quantum_mps", feature = "quantum_bptn", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_sv"))] +#[cfg(any(feature = "quantum_mps", feature = "quantum_bptn", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_sv", feature = "quantum_dmrg"))] mod linalg; // Gates as matrices, shared by the state vector and the tensor networks. @@ -38,7 +38,8 @@ pub mod gates; feature = "quantum_sv", feature = "quantum_arch", feature = "quantum_qio", - feature = "quantum_dmrg" + feature = "quantum_dmrg", + feature = "quantum_tdvp" ))] mod repro; @@ -123,6 +124,9 @@ pub mod quantum_noisy; #[cfg(feature = "quantum_dmrg")] pub mod quantum_dmrg; +#[cfg(feature = "quantum_tdvp")] +pub mod quantum_tdvp; + #[cfg(feature = "quantum_frame")] pub mod quantum_match; diff --git a/wai/wai-quantum/src/linalg.rs b/wai/wai-quantum/src/linalg.rs index 235c3cc86..b9ed6925b 100644 --- a/wai/wai-quantum/src/linalg.rs +++ b/wai/wai-quantum/src/linalg.rs @@ -8,7 +8,7 @@ // Complex f64 (deterministic: only + - * / sqrt, no FMA/transcendentals) // =========================================================================== -#[derive(Clone, Copy, Debug)] +#[derive(Clone, Copy, Debug, PartialEq)] pub struct C { pub re: f64, pub im: f64, @@ -148,7 +148,7 @@ pub(crate) fn jacobi_svd(mut m: Vec>, rows: usize, cols: usize) -> (Vec>, rows: usize, cols: usize) -> (Vec>, Vec, Vec>) { // V starts as identity (cols × cols), stored by columns. let mut v: Vec> = (0..cols) diff --git a/wai/wai-quantum/src/quantum_dmrg.rs b/wai/wai-quantum/src/quantum_dmrg.rs index 723b0a750..451c2c2c6 100644 --- a/wai/wai-quantum/src/quantum_dmrg.rs +++ b/wai/wai-quantum/src/quantum_dmrg.rs @@ -26,16 +26,19 @@ //! free fermions: the ground energy of an open chain of any length, from //! the eigenvalues of an `n × n` matrix. //! -//! The eigensolver is Householder tridiagonalisation followed by the implicit -//! QL iteration, with hypotenuses taken in-crate. The effective Hamiltonian +//! The engine is generic over [`Scalar`]: real here, complex for time +//! evolution. The eigensolver is Householder tridiagonalisation followed by +//! the implicit QL iteration, with hypotenuses taken in-crate; +//! [`hermitian_eigen`] reduces a Hermitian matrix to a real tridiagonal one +//! first and shares the iteration. The effective Hamiltonian //! is applied in blocks of rows across threads, each row in the same order. //! Every sum runs in a fixed order, so a run is the same bits on every //! machine and at any thread count. //! //! # Checked //! -//! - **The eigensolver.** It reconstructs random symmetric matrices up to -//! 40 × 40 to `10⁻¹²` with orthonormal vectors. +//! - **The eigensolvers.** They reconstruct random symmetric and Hermitian +//! matrices up to 40 × 40 to `10⁻¹²` with orthonormal vectors. //! - **The free-fermion solution.** It equals exact diagonalisation of the //! Ising chain at three couplings, for 4, 7 and 10 sites. //! - **DMRG against exact diagonalisation.** On Heisenberg chains (isotropic, @@ -60,17 +63,124 @@ //! # Honest boundaries //! //! - **One dimension, nearest neighbours, real Hamiltonians.** Longer-range -//! couplings, complex terms and two-dimensional lattices are not handled, -//! and nor is time evolution. +//! couplings, complex terms and two-dimensional lattices are not handled. +//! Time evolution is the TDVP module's, on this engine. //! - **No symmetry sectors.** Conserved quantities are not used to block the //! matrices, so a run costs more than a symmetric code would. //! - **Variational.** The energy is an upper bound, converging from above as //! `max_bond` grows. Critical chains converge slowest. +pub use crate::linalg::C; use crate::repro::ln; // --------------------------------------------------------------------------- -// The symmetric eigensolver +// Scalars +// --------------------------------------------------------------------------- + +/// The numbers a matrix-product state is made of: real for the ground states +/// of real Hamiltonians, complex for time evolution. Real arithmetic runs in +/// exactly the order it did before the engine was generic, so its bits are +/// unchanged. +pub trait Scalar: Copy + PartialEq + core::fmt::Debug + Send + Sync + 'static { + const ZERO: Self; + const ONE: Self; + fn add(self, o: Self) -> Self; + fn sub(self, o: Self) -> Self; + fn mul(self, o: Self) -> Self; + fn conj(self) -> Self; + /// `self · c` for a real `c`. + fn scale(self, c: f64) -> Self; + /// `|self|²`. + fn norm2(self) -> f64; + fn re(self) -> f64; + fn is_zero(self) -> bool; + /// Eigenvalues (ascending) and eigenvectors (column `i` of the returned + /// row-major matrix) of a Hermitian `n × n` matrix. + fn eigh(a: &[Self], n: usize) -> (Vec, Vec); +} + +impl Scalar for f64 { + const ZERO: f64 = 0.0; + const ONE: f64 = 1.0; + #[inline] + fn add(self, o: f64) -> f64 { + self + o + } + #[inline] + fn sub(self, o: f64) -> f64 { + self - o + } + #[inline] + fn mul(self, o: f64) -> f64 { + self * o + } + #[inline] + fn conj(self) -> f64 { + self + } + #[inline] + fn scale(self, c: f64) -> f64 { + self * c + } + #[inline] + fn norm2(self) -> f64 { + self * self + } + #[inline] + fn re(self) -> f64 { + self + } + #[inline] + fn is_zero(self) -> bool { + self == 0.0 + } + fn eigh(a: &[f64], n: usize) -> (Vec, Vec) { + symmetric_eigen(a, n) + } +} + +impl Scalar for C { + const ZERO: C = C { re: 0.0, im: 0.0 }; + const ONE: C = C { re: 1.0, im: 0.0 }; + #[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, c: f64) -> C { + C { re: self.re * c, im: self.im * c } + } + #[inline] + fn norm2(self) -> f64 { + self.re * self.re + self.im * self.im + } + #[inline] + fn re(self) -> f64 { + self.re + } + #[inline] + fn is_zero(self) -> bool { + self.re == 0.0 && self.im == 0.0 + } + fn eigh(a: &[C], n: usize) -> (Vec, Vec) { + hermitian_eigen(a, n) + } +} + +// --------------------------------------------------------------------------- +// The symmetric and Hermitian eigensolvers // --------------------------------------------------------------------------- fn hypot(a: f64, b: f64) -> f64 { @@ -188,7 +298,17 @@ pub fn symmetric_eigen(a: &[f64], n: usize) -> (Vec, Vec) { } v[n - 1][n - 1] = 1.0; e[0] = 0.0; - // Implicit QL on the tridiagonal matrix. + tql2(&mut d, &mut e, &mut v); + (d, v.into_iter().flatten().collect()) +} + +/// The implicit QL iteration on the tridiagonal matrix with diagonal `d` and +/// subdiagonal `e[i] = T[i][i − 1]` (`e[0]` unused), applied to the rows of +/// `v`. On return `d` holds the eigenvalues ascending and the columns of `v` +/// the eigenvectors, carried through `v`'s starting transformation. +#[allow(clippy::needless_range_loop)] +fn tql2(d: &mut [f64], e: &mut [f64], v: &mut [Vec]) { + let n = d.len(); for i in 1..n { e[i - 1] = e[i]; } @@ -271,7 +391,98 @@ pub fn symmetric_eigen(a: &[f64], n: usize) -> (Vec, Vec) { } } } - (d, v.into_iter().flatten().collect()) +} + +/// Eigenvalues (ascending) and eigenvectors of a Hermitian `n × n` matrix +/// given row-major. Eigenvector `i` is column `i` of the returned row-major +/// matrix. +/// +/// Householder reflections chosen to leave a real subdiagonal reduce the +/// matrix to a real symmetric tridiagonal one, `Q† A Q = T`. The QL iteration +/// diagonalises `T`, and the reflections carry its eigenvectors back. +#[allow(clippy::needless_range_loop)] +pub fn hermitian_eigen(a: &[C], n: usize) -> (Vec, Vec) { + if n == 0 { + return (Vec::new(), Vec::new()); + } + let mut a = a.to_vec(); + let mut d = vec![0.0; n]; + let mut e = vec![0.0; n]; + // (tau, v) for each reflection H = I − tau v v†, acting on rows i + 1..n. + let mut reflections: Vec<(C, Vec)> = Vec::with_capacity(n.saturating_sub(1)); + for i in 0..n - 1 { + let m = n - i - 1; + let alpha = a[(i + 1) * n + i]; + let xnorm2: f64 = (i + 2..n).map(|r| a[r * n + i].norm2()).sum(); + if xnorm2 == 0.0 && alpha.im == 0.0 { + // Already real and tridiagonal in this column. + e[i + 1] = alpha.re; + reflections.push((C::ZERO, Vec::new())); + } else { + let norm = (alpha.re * alpha.re + alpha.im * alpha.im + xnorm2).sqrt(); + let beta = if alpha.re >= 0.0 { -norm } else { norm }; + let tau = C { re: (beta - alpha.re) / beta, im: -alpha.im / beta }; + // v = (1, x / (alpha − beta)). + let z = C { re: alpha.re - beta, im: alpha.im }; + let zz = z.re * z.re + z.im * z.im; + let inv = C { re: z.re / zz, im: -z.im / zz }; + let mut v = Vec::with_capacity(m); + v.push(C::ONE); + for r in i + 2..n { + v.push(Scalar::mul(a[r * n + i], inv)); + } + e[i + 1] = beta; + // p = tau · B v over the trailing block B = A[i+1.., i+1..]. + let mut p = vec![C::ZERO; m]; + for r in 0..m { + let row = &a[(i + 1 + r) * n + i + 1..(i + 1 + r) * n + n]; + let mut s = C::ZERO; + for (x, y) in row.iter().zip(&v) { + s = Scalar::add(s, Scalar::mul(*x, *y)); + } + p[r] = Scalar::mul(tau, s); + } + // w = p − ½ tau (p† v) v + let mut pv = C::ZERO; + for (x, y) in p.iter().zip(&v) { + pv = Scalar::add(pv, Scalar::mul(Scalar::conj(*x), *y)); + } + let k = Scalar::scale(Scalar::mul(tau, pv), -0.5); + let w: Vec = p.iter().zip(&v).map(|(x, y)| Scalar::add(*x, Scalar::mul(k, *y))).collect(); + // B −= v w† + w v† + for r in 0..m { + for c in 0..m { + let upd = Scalar::add(Scalar::mul(v[r], Scalar::conj(w[c])), Scalar::mul(w[r], Scalar::conj(v[c]))); + let at = (i + 1 + r) * n + i + 1 + c; + a[at] = Scalar::sub(a[at], upd); + } + } + reflections.push((tau, v)); + } + d[i] = a[i * n + i].re; + } + d[n - 1] = a[(n - 1) * n + n - 1].re; + let mut z: Vec> = (0..n).map(|r| (0..n).map(|c| if r == c { 1.0 } else { 0.0 }).collect()).collect(); + tql2(&mut d, &mut e, &mut z); + // X = H_0 H_1 … H_{n−2} Z. + let mut x: Vec = z.into_iter().flatten().map(|re| C { re, im: 0.0 }).collect(); + for (i, (tau, v)) in reflections.iter().enumerate().rev() { + if v.is_empty() { + continue; + } + for c in 0..n { + let mut s = C::ZERO; + for (t, vt) in v.iter().enumerate() { + s = Scalar::add(s, Scalar::mul(Scalar::conj(*vt), x[(i + 1 + t) * n + c])); + } + let ts = Scalar::mul(*tau, s); + for (t, vt) in v.iter().enumerate() { + let at = (i + 1 + t) * n + c; + x[at] = Scalar::sub(x[at], Scalar::mul(*vt, ts)); + } + } + } + (d, x) } // --------------------------------------------------------------------------- @@ -434,6 +645,51 @@ fn lanczos(apply: &dyn Fn(&[f64], &mut [f64]), start: &[f64], krylov: usize, res (energy, v) } +/// `out = H inp` for `chain` on all `2ⁿ` amplitudes, site 0 the most +/// significant bit. +pub(crate) fn apply_dense(chain: &Chain, inp: &[f64], out: &mut [f64]) { + let n = chain.n; + let bit = |x: usize, i: usize| (x >> (n - 1 - i)) & 1; + out.iter_mut().for_each(|o| *o = 0.0); + for (x, &) in inp.iter().enumerate() { + if amp == 0.0 { + continue; + } + for i in 0..n { + let s = bit(x, i); + for &(c, o) in &chain.sites { + for sp in 0..2 { + let m = o.0[sp * 2 + s]; + if m != 0.0 { + let y = (x & !(1 << (n - 1 - i))) | (sp << (n - 1 - i)); + out[y] += c * m * amp; + } + } + } + if i + 1 < n { + let s2 = bit(x, i + 1); + for &(c, a, b) in &chain.bonds { + for ap in 0..2 { + let ma = a.0[ap * 2 + s]; + if ma == 0.0 { + continue; + } + for bp in 0..2 { + let mb = b.0[bp * 2 + s2]; + if mb == 0.0 { + continue; + } + let mask = (1 << (n - 1 - i)) | (1 << (n - 2 - i)); + let y = (x & !mask) | (ap << (n - 1 - i)) | (bp << (n - 2 - i)); + out[y] += c * ma * mb * amp; + } + } + } + } + } + } +} + /// The exact ground energy of `chain` by Lanczos on all `2ⁿ` states (at most /// 22 sites). pub fn exact_ground_energy(chain: &Chain) -> Option { @@ -442,47 +698,7 @@ pub fn exact_ground_energy(chain: &Chain) -> Option { return None; } let dim = 1usize << n; - let bit = |x: usize, i: usize| (x >> (n - 1 - i)) & 1; - let apply = |inp: &[f64], out: &mut [f64]| { - out.iter_mut().for_each(|o| *o = 0.0); - for (x, &) in inp.iter().enumerate() { - if amp == 0.0 { - continue; - } - for i in 0..n { - let s = bit(x, i); - for &(c, o) in &chain.sites { - for sp in 0..2 { - let m = o.0[sp * 2 + s]; - if m != 0.0 { - let y = (x & !(1 << (n - 1 - i))) | (sp << (n - 1 - i)); - out[y] += c * m * amp; - } - } - } - if i + 1 < n { - let s2 = bit(x, i + 1); - for &(c, a, b) in &chain.bonds { - for ap in 0..2 { - let ma = a.0[ap * 2 + s]; - if ma == 0.0 { - continue; - } - for bp in 0..2 { - let mb = b.0[bp * 2 + s2]; - if mb == 0.0 { - continue; - } - let mask = (1 << (n - 1 - i)) | (1 << (n - 2 - i)); - let y = (x & !mask) | (ap << (n - 1 - i)) | (bp << (n - 2 - i)); - out[y] += c * ma * mb * amp; - } - } - } - } - } - } - }; + let apply = |inp: &[f64], out: &mut [f64]| apply_dense(chain, inp, out); // u64, not usize: the start vector must not depend on the pointer width. let start: Vec = (0..dim as u64).map(|x| 1.0 + (x.wrapping_mul(2_654_435_761) % 1000) as f64 / 1000.0).collect(); Some(lanczos(&apply, &start, 60, 20, 1e-10).0) @@ -533,15 +749,16 @@ pub fn ising_exact_energy(n: usize, j: f64, h: f64) -> f64 { // Matrix-product states and DMRG // --------------------------------------------------------------------------- -/// A real matrix-product state: site `i` is `left × 2 × right`, stored -/// `[(a·2 + s)·right + b]`. +/// A matrix-product state: site `i` is `left × 2 × right`, stored +/// `[(a·2 + s)·right + b]`. Real for ground states, complex ([`C`]) for time +/// evolution. #[derive(Clone, Debug, PartialEq)] -pub struct Mps { +pub struct Mps { pub dims: Vec, - pub sites: Vec>, + pub sites: Vec>, } -impl Mps { +impl Mps { /// A deterministic pseudo-random state with bond dimension `bond` (capped /// by what the chain allows). pub fn random(n: usize, bond: usize, seed: u64) -> Mps { @@ -562,9 +779,71 @@ impl Mps { Mps { dims, sites } } + /// The same state with complex entries. + pub fn to_complex(&self) -> Mps { + Mps { dims: self.dims.clone(), sites: self.sites.iter().map(|s| s.iter().map(|&re| C { re, im: 0.0 }).collect()).collect() } + } +} + +impl Mps { + /// The product state `⊗ (a_i |↑⟩ + b_i |↓⟩)`, each factor normalised. + pub fn product(factors: &[[C; 2]]) -> Mps { + let sites = factors + .iter() + .map(|f| { + let norm = (f[0].norm2() + f[1].norm2()).sqrt(); + assert!(norm > 0.0, "a product factor must be non-zero"); + vec![Scalar::scale(f[0], 1.0 / norm), Scalar::scale(f[1], 1.0 / norm)] + }) + .collect(); + Mps { dims: vec![1; factors.len() + 1], sites } + } + + /// The basis state with site `i` down where `down[i]`. + pub fn basis(down: &[bool]) -> Mps { + let f: Vec<[C; 2]> = down.iter().map(|&d| if d { [C::ZERO, C::ONE] } else { [C::ONE, C::ZERO] }).collect(); + Mps::product(&f) + } +} + +impl Mps { pub fn n(&self) -> usize { self.sites.len() } + + /// All `2ⁿ` amplitudes (at most 24 sites), site 0 the most significant + /// bit and `|↑⟩ = 0`. + pub fn to_dense(&self) -> Option> { + let n = self.n(); + if n == 0 || n > 24 { + return None; + } + // Rows are the basis states of the sites so far; columns the bond. + let mut cur = vec![T::ONE]; + let mut rows = 1usize; + for i in 0..n { + let (dl, dr) = (self.dims[i], self.dims[i + 1]); + let site = &self.sites[i]; + let mut next = vec![T::ZERO; rows * 2 * dr]; + for r in 0..rows { + for a in 0..dl { + let c = cur[r * dl + a]; + if c.is_zero() { + continue; + } + for s in 0..2 { + for b in 0..dr { + let at = (r * 2 + s) * dr + b; + next[at] = next[at].add(c.mul(site[(a * 2 + s) * dr + b])); + } + } + } + } + cur = next; + rows *= 2; + } + Some(cur) + } } /// Settings for [`dmrg`]. @@ -587,10 +866,14 @@ pub struct DmrgConfig { pub threads: usize, } +/// The machine's thread count; one in WebAssembly. +pub(crate) fn default_threads() -> usize { + if cfg!(target_arch = "wasm32") { 1 } else { std::thread::available_parallelism().map_or(1, |n| n.get()) } +} + impl Default for DmrgConfig { fn default() -> DmrgConfig { - let threads = if cfg!(target_arch = "wasm32") { 1 } else { std::thread::available_parallelism().map_or(1, |n| n.get()) }; - DmrgConfig { max_bond: 64, cutoff: 1e-12, sweeps: 6, krylov: 24, restarts: 4, tol: 1e-8, seed: 1, threads } + DmrgConfig { max_bond: 64, cutoff: 1e-12, sweeps: 6, krylov: 24, restarts: 4, tol: 1e-8, seed: 1, threads: default_threads() } } } @@ -608,13 +891,18 @@ pub struct DmrgResult { pub mps: Mps, } -/// Contract environments and the two-site effective Hamiltonian. -struct Engine<'a> { +/// Contract environments and effective Hamiltonians. Environments are +/// `L[a', w, a]` and `R[b', w, b]`, bra index first. +pub(crate) struct Engine<'a> { mpo: &'a Mpo, w: usize, } impl Engine<'_> { + pub(crate) fn new(mpo: &Mpo) -> Engine<'_> { + Engine { mpo, w: mpo.w } + } + fn op(&self, a: usize, b: usize) -> &Op { &self.mpo.site[a * self.w + b] } @@ -623,30 +911,44 @@ impl Engine<'_> { self.op(a, b).0.iter().any(|&v| v != 0.0) } - /// `L'[b', w2, b] = Σ A[a',s',b'] L[a',w1,a] W[w1,w2,s',s] A[a,s,b]`. - fn grow_left(&self, l: &[f64], a: &[f64], dl: usize, dr: usize) -> Vec { + /// The environment at the left edge. + pub(crate) fn left_edge(&self) -> Vec { + let mut l = vec![T::ZERO; self.w]; + l[0] = T::ONE; + l + } + + /// The environment at the right edge. + pub(crate) fn right_edge(&self) -> Vec { + let mut r = vec![T::ZERO; self.w]; + r[self.w - 1] = T::ONE; + r + } + + /// `L'[b', w2, b] = Σ A*[a',s',b'] L[a',w1,a] W[w1,w2,s',s] A[a,s,b]`. + pub(crate) fn grow_left(&self, l: &[T], a: &[T], dl: usize, dr: usize) -> Vec { let w = self.w; // T[a', w1, s, b] = Σ_a L[a', w1, a] A[a, s, b] - let mut t = vec![0.0; dl * w * 2 * dr]; + let mut t = vec![T::ZERO; dl * w * 2 * dr]; for ap in 0..dl { for w1 in 0..w { let lrow = &l[(ap * w + w1) * dl..(ap * w + w1 + 1) * dl]; for (aa, &lv) in lrow.iter().enumerate() { - if lv == 0.0 { + if lv.is_zero() { continue; } for s in 0..2 { let src = &a[(aa * 2 + s) * dr..(aa * 2 + s + 1) * dr]; let dst = &mut t[((ap * w + w1) * 2 + s) * dr..((ap * w + w1) * 2 + s + 1) * dr]; for (d, &x) in dst.iter_mut().zip(src) { - *d += lv * x; + *d = d.add(lv.mul(x)); } } } } } // U[a', w2, s', b] = Σ_{w1, s} W[w1, w2, s', s] T[a', w1, s, b] - let mut u = vec![0.0; dl * w * 2 * dr]; + let mut u = vec![T::ZERO; dl * w * 2 * dr]; for ap in 0..dl { for w1 in 0..w { for w2 in 0..w { @@ -663,27 +965,28 @@ impl Engine<'_> { let src = ((ap * w + w1) * 2 + s) * dr; let dst = ((ap * w + w2) * 2 + sp) * dr; for b in 0..dr { - u[dst + b] += m * t[src + b]; + u[dst + b] = u[dst + b].add(t[src + b].scale(m)); } } } } } } - // L'[b', w2, b] = Σ_{a', s'} A[a', s', b'] U[a', w2, s', b] - let mut out = vec![0.0; dr * w * dr]; + // L'[b', w2, b] = Σ_{a', s'} A*[a', s', b'] U[a', w2, s', b] + let mut out = vec![T::ZERO; dr * w * dr]; for ap in 0..dl { for sp in 0..2 { let arow = &a[(ap * 2 + sp) * dr..(ap * 2 + sp + 1) * dr]; for (bp, &av) in arow.iter().enumerate() { - if av == 0.0 { + if av.is_zero() { continue; } + let av = av.conj(); for w2 in 0..w { let src = &u[((ap * w + w2) * 2 + sp) * dr..((ap * w + w2) * 2 + sp + 1) * dr]; let dst = &mut out[(bp * w + w2) * dr..(bp * w + w2 + 1) * dr]; for (d, &x) in dst.iter_mut().zip(src) { - *d += av * x; + *d = d.add(av.mul(x)); } } } @@ -692,20 +995,20 @@ impl Engine<'_> { out } - /// `R'[a', w1, a] = Σ B[a',s',b'] W[w1,w2,s',s] R[b',w2,b] B[a,s,b]`. - fn grow_right(&self, r: &[f64], bt: &[f64], dl: usize, dr: usize) -> Vec { + /// `R'[a', w1, a] = Σ B*[a',s',b'] W[w1,w2,s',s] R[b',w2,b] B[a,s,b]`. + pub(crate) fn grow_right(&self, r: &[T], bt: &[T], dl: usize, dr: usize) -> Vec { let w = self.w; // T[a, s, b', w2] = Σ_b B[a, s, b] R[b', w2, b] - let mut t = vec![0.0; dl * 2 * dr * w]; + let mut t = vec![T::ZERO; dl * 2 * dr * w]; for a in 0..dl { for s in 0..2 { let brow = &bt[(a * 2 + s) * dr..(a * 2 + s + 1) * dr]; for bp in 0..dr { for w2 in 0..w { let rrow = &r[(bp * w + w2) * dr..(bp * w + w2 + 1) * dr]; - let mut acc = 0.0; + let mut acc = T::ZERO; for (x, y) in brow.iter().zip(rrow) { - acc += x * y; + acc = acc.add(x.mul(*y)); } t[((a * 2 + s) * dr + bp) * w + w2] = acc; } @@ -713,7 +1016,7 @@ impl Engine<'_> { } } // U[a, s', b', w1] = Σ_{w2, s} W[w1, w2, s', s] T[a, s, b', w2] - let mut u = vec![0.0; dl * 2 * dr * w]; + let mut u = vec![T::ZERO; dl * 2 * dr * w]; for w1 in 0..w { for w2 in 0..w { if !self.nonzero(w1, w2) { @@ -728,25 +1031,28 @@ impl Engine<'_> { } for a in 0..dl { for bp in 0..dr { - u[((a * 2 + sp) * dr + bp) * w + w1] += m * t[((a * 2 + s) * dr + bp) * w + w2]; + let at = ((a * 2 + sp) * dr + bp) * w + w1; + u[at] = u[at].add(t[((a * 2 + s) * dr + bp) * w + w2].scale(m)); } } } } } } - // R'[a', w1, a] = Σ_{s', b'} B[a', s', b'] U[a, s', b', w1] - let mut out = vec![0.0; dl * w * dl]; + // R'[a', w1, a] = Σ_{s', b'} B*[a', s', b'] U[a, s', b', w1] + let mut out = vec![T::ZERO; dl * w * dl]; for ap in 0..dl { for sp in 0..2 { for bp in 0..dr { let bv = bt[(ap * 2 + sp) * dr + bp]; - if bv == 0.0 { + if bv.is_zero() { continue; } + let bv = bv.conj(); for a in 0..dl { for w1 in 0..w { - out[(ap * w + w1) * dl + a] += bv * u[((a * 2 + sp) * dr + bp) * w + w1]; + let at = (ap * w + w1) * dl + a; + out[at] = out[at].add(bv.mul(u[((a * 2 + sp) * dr + bp) * w + w1])); } } } @@ -755,33 +1061,33 @@ impl Engine<'_> { out } - /// Rows `rows` (bra indices `a'`) of `H_eff θ`, written to `out` (which - /// holds just those rows). + /// Rows `rows` (bra indices `a'`) of `H_eff θ` for the two-site tensor + /// `θ[a, s1, s2, b]`, written to `out` (which holds just those rows). #[allow(clippy::too_many_arguments)] - fn apply_rows(&self, l: &[f64], r: &[f64], theta: &[f64], out: &mut [f64], dl: usize, dr: usize, rows: core::ops::Range) { + fn apply_rows(&self, l: &[T], r: &[T], theta: &[T], out: &mut [T], dl: usize, dr: usize, rows: core::ops::Range) { let w = self.w; let first = rows.start; let nrows = rows.len(); // T1[a', w0, s1, s2, b] = Σ_a L[a', w0, a] θ[a, s1, s2, b] let blk = 4 * dr; - let mut t1 = vec![0.0; nrows * w * blk]; + let mut t1 = vec![T::ZERO; nrows * w * blk]; for ap in 0..nrows { for w0 in 0..w { let lrow = &l[((first + ap) * w + w0) * dl..((first + ap) * w + w0 + 1) * dl]; let dst = &mut t1[(ap * w + w0) * blk..(ap * w + w0 + 1) * blk]; for (a, &lv) in lrow.iter().enumerate() { - if lv == 0.0 { + if lv.is_zero() { continue; } for (d, &x) in dst.iter_mut().zip(&theta[a * blk..(a + 1) * blk]) { - *d += lv * x; + *d = d.add(lv.mul(x)); } } } } // T2[a', w1, s1', s2, b] = Σ_{w0, s1} W[w0, w1, s1', s1] T1[a', w0, s1, s2, b] let half = 2 * dr; - let mut t2 = vec![0.0; nrows * w * blk]; + let mut t2 = vec![T::ZERO; nrows * w * blk]; for ap in 0..nrows { for w0 in 0..w { for w1 in 0..w { @@ -798,7 +1104,7 @@ impl Engine<'_> { let src = (ap * w + w0) * blk + s1 * half; let dst = (ap * w + w1) * blk + s1p * half; for k in 0..half { - t2[dst + k] += m * t1[src + k]; + t2[dst + k] = t2[dst + k].add(t1[src + k].scale(m)); } } } @@ -806,7 +1112,7 @@ impl Engine<'_> { } } // T3[a', s1', w2, s2', b] = Σ_{w1, s2} W[w1, w2, s2', s2] T2[a', w1, s1', s2, b] - let mut t3 = vec![0.0; nrows * 2 * w * 2 * dr]; + let mut t3 = vec![T::ZERO; nrows * 2 * w * 2 * dr]; for ap in 0..nrows { for w1 in 0..w { for w2 in 0..w { @@ -824,7 +1130,7 @@ impl Engine<'_> { let src = (ap * w + w1) * blk + s1p * half + s2 * dr; let dst = (((ap * 2 + s1p) * w + w2) * 2 + s2p) * dr; for b in 0..dr { - t3[dst + b] += m * t2[src + b]; + t3[dst + b] = t3[dst + b].add(t2[src + b].scale(m)); } } } @@ -837,12 +1143,12 @@ impl Engine<'_> { for s1p in 0..2 { for s2p in 0..2 { for bp in 0..dr { - let mut acc = 0.0; + let mut acc = T::ZERO; for w2 in 0..w { let trow = &t3[(((ap * 2 + s1p) * w + w2) * 2 + s2p) * dr..(((ap * 2 + s1p) * w + w2) * 2 + s2p + 1) * dr]; let rrow = &r[(bp * w + w2) * dr..(bp * w + w2 + 1) * dr]; for (x, y) in trow.iter().zip(rrow) { - acc += x * y; + acc = acc.add(x.mul(*y)); } } out[((ap * 2 + s1p) * 2 + s2p) * dr + bp] = acc; @@ -852,19 +1158,85 @@ impl Engine<'_> { } } - /// `out = H_eff θ` for `θ[a, s1, s2, b]` between `l` (dl) and `r` (dr). - /// Blocks of bra rows run on separate threads; each row is computed in - /// the same order whatever the split, so the result is the same bits. + /// Rows `rows` of `H_eff C` for the one-site tensor `C[a, s, b]`. + #[cfg(feature = "quantum_tdvp")] #[allow(clippy::too_many_arguments)] - fn apply(&self, l: &[f64], r: &[f64], theta: &[f64], out: &mut [f64], dl: usize, dr: usize, threads: usize) { - let row = 4 * dr; - // At least eight rows a thread; small bonds stay on one thread, where - // spawning would cost more than it saves. - let threads = if cfg!(target_arch = "wasm32") { 1 } else { threads.max(1).min(dl / 8) }; - if threads <= 1 || dl * dr < 2048 { + fn apply1_rows(&self, l: &[T], r: &[T], c: &[T], out: &mut [T], dl: usize, dr: usize, rows: core::ops::Range) { + let w = self.w; + let first = rows.start; + let nrows = rows.len(); + // T1[a', w0, s, b] = Σ_a L[a', w0, a] C[a, s, b] + let blk = 2 * dr; + let mut t1 = vec![T::ZERO; nrows * w * blk]; + for ap in 0..nrows { + for w0 in 0..w { + let lrow = &l[((first + ap) * w + w0) * dl..((first + ap) * w + w0 + 1) * dl]; + let dst = &mut t1[(ap * w + w0) * blk..(ap * w + w0 + 1) * blk]; + for (a, &lv) in lrow.iter().enumerate() { + if lv.is_zero() { + continue; + } + for (d, &x) in dst.iter_mut().zip(&c[a * blk..(a + 1) * blk]) { + *d = d.add(lv.mul(x)); + } + } + } + } + // T2[a', s', w1, b] = Σ_{w0, s} W[w0, w1, s', s] T1[a', w0, s, b] + let mut t2 = vec![T::ZERO; nrows * 2 * w * dr]; + for ap in 0..nrows { + for w0 in 0..w { + for w1 in 0..w { + if !self.nonzero(w0, w1) { + continue; + } + let o = self.op(w0, w1).0; + for sp in 0..2 { + for s in 0..2 { + let m = o[sp * 2 + s]; + if m == 0.0 { + continue; + } + let src = (ap * w + w0) * blk + s * dr; + let dst = ((ap * 2 + sp) * w + w1) * dr; + for b in 0..dr { + t2[dst + b] = t2[dst + b].add(t1[src + b].scale(m)); + } + } + } + } + } + } + // out[a', s', b'] = Σ_{w1, b} T2[a', s', w1, b] R[b', w1, b] + for ap in 0..nrows { + for sp in 0..2 { + for bp in 0..dr { + let mut acc = T::ZERO; + for w1 in 0..w { + let trow = &t2[((ap * 2 + sp) * w + w1) * dr..((ap * 2 + sp) * w + w1 + 1) * dr]; + let rrow = &r[(bp * w + w1) * dr..(bp * w + w1 + 1) * dr]; + for (x, y) in trow.iter().zip(rrow) { + acc = acc.add(x.mul(*y)); + } + } + out[(ap * 2 + sp) * dr + bp] = acc; + } + } + } + } + + /// `out = H_eff θ` for the two-site `θ[a, s1, s2, b]` between `l` (dl) + /// and `r` (dr). Blocks of bra rows run on separate threads; each row is + /// computed in the same order whatever the split, so the result is the + /// same bits. + #[allow(clippy::too_many_arguments)] + pub(crate) fn apply(&self, l: &[T], r: &[T], theta: &[T], out: &mut [T], dl: usize, dr: usize, threads: usize) { + let threads = Self::split_threads(threads, dl, dr); + if threads <= 1 { self.apply_rows(l, r, theta, out, dl, dr, 0..dl); return; } + let row = 4 * dr; let per = dl.div_ceil(threads); std::thread::scope(|scope| { for (c, chunk) in out.chunks_mut(per * row).enumerate() { @@ -874,25 +1246,53 @@ impl Engine<'_> { } }); } + + /// `out = H_eff C` for the one-site `C[a, s, b]`, threaded like + /// [`Engine::apply`]. + #[cfg(feature = "quantum_tdvp")] + #[allow(clippy::too_many_arguments)] + pub(crate) fn apply1(&self, l: &[T], r: &[T], c: &[T], out: &mut [T], dl: usize, dr: usize, threads: usize) { + let threads = Self::split_threads(threads, dl, dr); + if threads <= 1 { + self.apply1_rows(l, r, c, out, dl, dr, 0..dl); + return; + } + let row = 2 * dr; + let per = dl.div_ceil(threads); + std::thread::scope(|scope| { + for (k, chunk) in out.chunks_mut(per * row).enumerate() { + let start = k * per; + let end = (start + per).min(dl); + scope.spawn(move || self.apply1_rows(l, r, c, chunk, dl, dr, start..end)); + } + }); + } + + /// At least eight rows a thread; small bonds stay on one thread, where + /// spawning would cost more than it saves. + fn split_threads(threads: usize, dl: usize, dr: usize) -> usize { + let threads = if cfg!(target_arch = "wasm32") { 1 } else { threads.max(1).min(dl / 8) }; + if dl * dr < 2048 { 1 } else { threads } + } } -/// Split `m` (rows × cols) by the eigenvectors of `m mᵀ`, keeping at most -/// `max_bond` and dropping weight below `cutoff`. Returns the kept vectors -/// (rows × k), `uᵀ m` (k × cols), the discarded weight and the kept -/// normalised weights. -fn split(m: &[f64], rows: usize, cols: usize, max_bond: usize, cutoff: f64) -> (Vec, Vec, f64, Vec) { - let mut rho = vec![0.0; rows * rows]; +/// Split `m` (rows × cols) by the eigenvectors of `m m†`, keeping at most +/// `max_bond` (and never more than `m` has rank for) and dropping weight +/// below `cutoff`. Returns the kept vectors (rows × k), `u† m` (k × cols), +/// the discarded weight and the kept normalised weights. +pub(crate) fn split(m: &[T], rows: usize, cols: usize, max_bond: usize, cutoff: f64) -> (Vec, Vec, f64, Vec) { + let mut rho = vec![T::ZERO; rows * rows]; for i in 0..rows { for j in i..rows { - let mut s = 0.0; + let mut s = T::ZERO; for c in 0..cols { - s += m[i * cols + c] * m[j * cols + c]; + s = s.add(m[i * cols + c].mul(m[j * cols + c].conj())); } rho[i * rows + j] = s; - rho[j * rows + i] = s; + rho[j * rows + i] = s.conj(); } } - let (vals, vecs) = symmetric_eigen(&rho, rows); + let (vals, vecs) = T::eigh(&rho, rows); let total: f64 = vals.iter().map(|v| v.max(0.0)).sum(); // Largest first. let order: Vec = (0..rows).rev().collect(); @@ -900,14 +1300,14 @@ fn split(m: &[f64], rows: usize, cols: usize, max_bond: usize, cutoff: f64) -> ( let mut kept_weight = 0.0; for &idx in &order { let w = vals[idx].max(0.0) / total; - if k >= max_bond || (k > 0 && w < cutoff) { + if k >= max_bond.min(cols) || (k > 0 && w < cutoff) { break; } kept_weight += w; k += 1; } let k = k.max(1); - let mut u = vec![0.0; rows * k]; + let mut u = vec![T::ZERO; rows * k]; let mut weights = Vec::with_capacity(k); for (c, &idx) in order.iter().take(k).enumerate() { weights.push(vals[idx].max(0.0) / total); @@ -915,23 +1315,24 @@ fn split(m: &[f64], rows: usize, cols: usize, max_bond: usize, cutoff: f64) -> ( u[r * k + c] = vecs[r * rows + idx]; } } - let mut rest = vec![0.0; k * cols]; + let mut rest = vec![T::ZERO; k * cols]; for c in 0..k { for r in 0..rows { let uv = u[r * k + c]; - if uv == 0.0 { + if uv.is_zero() { continue; } + let uv = uv.conj(); for j in 0..cols { - rest[c * cols + j] += uv * m[r * cols + j]; + rest[c * cols + j] = rest[c * cols + j].add(uv.mul(m[r * cols + j])); } } } (u, rest, (1.0 - kept_weight).max(0.0), weights) } -fn transpose(m: &[f64], rows: usize, cols: usize) -> Vec { - let mut t = vec![0.0; rows * cols]; +pub(crate) fn transpose(m: &[T], rows: usize, cols: usize) -> Vec { + let mut t = vec![T::ZERO; rows * cols]; for r in 0..rows { for c in 0..cols { t[c * rows + r] = m[r * cols + c]; @@ -940,48 +1341,32 @@ fn transpose(m: &[f64], rows: usize, cols: usize) -> Vec { t } -fn entropy(weights: &[f64]) -> f64 { +pub(crate) fn entropy(weights: &[f64]) -> f64 { -weights.iter().filter(|&&p| p > 0.0).map(|&p| p * ln(p)).sum::() } -/// Two-site DMRG for the ground state of `chain`. -pub fn dmrg(chain: &Chain, cfg: &DmrgConfig) -> DmrgResult { - let n = chain.n; - assert!(n >= 2, "DMRG needs at least two sites"); - let mpo = chain.mpo(); - let eng = Engine { mpo: &mpo, w: mpo.w }; - let w = mpo.w; - let mut mps = Mps::random(n, cfg.max_bond.min(8), cfg.seed); - // Bring the state into right-canonical form from the right end, and - // build the right environments. - let left_edge = { - let mut l = vec![0.0; w]; - l[0] = 1.0; - l - }; - let right_edge = { - let mut r = vec![0.0; w]; - r[w - 1] = 1.0; - r - }; +/// Bring `mps` into right-canonical form: every site but the first has +/// orthonormal rows, and the first carries the whole state. +pub(crate) fn right_canonicalise(mps: &mut Mps) { + let n = mps.n(); for i in (1..n).rev() { let (dl, dr) = (mps.dims[i], mps.dims[i + 1]); - // Site as dl × (2·dr); orthonormalise its rows. + // Site as dl × (2·dr); orthonormalise its rows through the + // eigenvectors of its transpose. let m = mps.sites[i].clone(); let (u, rest, _, _) = split(&transpose(&m, dl, 2 * dr), 2 * dr, dl, dl, 0.0); - // rest is k × dl = (uᵀ mᵀ); new site = uᵀ (k × 2dr), absorb mᵀ... keep it simple: let k = rest.len() / dl; let new_site = transpose(&u, 2 * dr, k); let carry = transpose(&rest, k, dl); // dl × k - // Previous site absorbs `carry` on its right index. + // The previous site absorbs `carry` on its right index. let (pl, pr) = (mps.dims[i - 1], mps.dims[i]); let prev = &mps.sites[i - 1]; - let mut np = vec![0.0; pl * 2 * k]; + let mut np = vec![T::ZERO; pl * 2 * k]; for a in 0..pl * 2 { for c in 0..k { - let mut s = 0.0; + let mut s = T::ZERO; for b in 0..pr { - s += prev[a * pr + b] * carry[b * k + c]; + s = s.add(prev[a * pr + b].mul(carry[b * k + c])); } np[a * k + c] = s; } @@ -990,13 +1375,25 @@ pub fn dmrg(chain: &Chain, cfg: &DmrgConfig) -> DmrgResult { mps.sites[i] = new_site; mps.dims[i] = k; } +} + +/// Two-site DMRG for the ground state of `chain`. +pub fn dmrg(chain: &Chain, cfg: &DmrgConfig) -> DmrgResult { + let n = chain.n; + assert!(n >= 2, "DMRG needs at least two sites"); + let mpo = chain.mpo(); + let eng = Engine::new(&mpo); + let mut mps = Mps::random(n, cfg.max_bond.min(8), cfg.seed); + // Bring the state into right-canonical form from the right end, and + // build the right environments. + right_canonicalise(&mut mps); let mut rights: Vec> = vec![Vec::new(); n + 1]; - rights[n] = right_edge.clone(); + rights[n] = eng.right_edge(); for i in (1..n).rev() { rights[i] = eng.grow_right(&rights[i + 1], &mps.sites[i], mps.dims[i], mps.dims[i + 1]); } let mut lefts: Vec> = vec![Vec::new(); n + 1]; - lefts[0] = left_edge; + lefts[0] = eng.left_edge(); let mut energy = 0.0; let mut sweep_energies = Vec::new(); let mut entropies = vec![0.0; n - 1]; @@ -1046,19 +1443,19 @@ pub fn dmrg(chain: &Chain, cfg: &DmrgConfig) -> DmrgResult { } /// `θ[a, s1, s2, b] = Σ_c A[a, s1, c] B[c, s2, b]`. -fn two_site(a: &[f64], b: &[f64], dl: usize, dm: usize, dr: usize) -> Vec { - let mut t = vec![0.0; dl * 4 * dr]; +pub(crate) fn two_site(a: &[T], b: &[T], dl: usize, dm: usize, dr: usize) -> Vec { + let mut t = vec![T::ZERO; dl * 4 * dr]; for x in 0..dl * 2 { for c in 0..dm { let av = a[x * dm + c]; - if av == 0.0 { + if av.is_zero() { continue; } for s2 in 0..2 { let src = &b[(c * 2 + s2) * dr..(c * 2 + s2 + 1) * dr]; let dst = &mut t[(x * 2 + s2) * dr..(x * 2 + s2 + 1) * dr]; for (d, &v) in dst.iter_mut().zip(src) { - *d += av * v; + *d = d.add(av.mul(v)); } } } @@ -1098,6 +1495,68 @@ mod tests { } } + #[test] + fn the_hermitian_eigensolver_reconstructs_its_matrix() { + let mut s = 11u64; + let mut next = move || { + s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + ((s >> 11) as f64 / 9_007_199_254_740_992.0) - 0.5 + }; + for n in [1usize, 2, 3, 6, 17, 40] { + let mut a = vec![C::ZERO; n * n]; + for i in 0..n { + a[i * n + i] = C { re: next(), im: 0.0 }; + for j in i + 1..n { + let v = C { re: next(), im: next() }; + a[i * n + j] = v; + a[j * n + i] = Scalar::conj(v); + } + } + // A column already tridiagonal and real must pass through too. + if n == 6 { + for r in 2..n { + a[r * n] = C::ZERO; + a[r] = C::ZERO; + } + a[n] = C { re: 0.25, im: 0.0 }; + a[1] = C { re: 0.25, im: 0.0 }; + } + let (vals, vecs) = hermitian_eigen(&a, n); + for w in vals.windows(2) { + assert!(w[0] <= w[1]); + } + for r in 0..n { + for c in 0..n { + let mut rec = C::ZERO; + let mut orth = C::ZERO; + for k in 0..n { + rec = Scalar::add(rec, Scalar::scale(Scalar::mul(vecs[r * n + k], Scalar::conj(vecs[c * n + k])), vals[k])); + orth = Scalar::add(orth, Scalar::mul(Scalar::conj(vecs[k * n + r]), vecs[k * n + c])); + } + let d = Scalar::sub(rec, a[r * n + c]); + assert!(d.norm2().sqrt() < 1e-12, "n={n} reconstruct ({r},{c})"); + let id = if r == c { 1.0 } else { 0.0 }; + assert!((orth.re - id).abs() < 1e-12 && orth.im.abs() < 1e-12, "n={n} orthonormal ({r},{c})"); + } + } + } + // On a real symmetric matrix it agrees with the real solver. + let n = 9; + let mut a = vec![0.0; n * n]; + for i in 0..n { + for j in i..n { + let v = next(); + a[i * n + j] = v; + a[j * n + i] = v; + } + } + let (real, _) = symmetric_eigen(&a, n); + let (herm, _) = hermitian_eigen(&a.iter().map(|&re| C { re, im: 0.0 }).collect::>(), n); + for (x, y) in real.iter().zip(&herm) { + assert!((x - y).abs() < 1e-13); + } + } + #[test] fn free_fermions_solve_the_ising_chain() { for n in [4usize, 7, 10] { diff --git a/wai/wai-quantum/src/quantum_tdvp.rs b/wai/wai-quantum/src/quantum_tdvp.rs new file mode 100644 index 000000000..74fffe8c6 --- /dev/null +++ b/wai/wai-quantum/src/quantum_tdvp.rs @@ -0,0 +1,1076 @@ +//! Time evolution of spin chains by TDVP — `wai.quantum.tdvp`. +//! +//! A claim that a quantum device simulates the dynamics of matter beyond +//! classical reach is a claim about time evolution. In one dimension the +//! time-dependent variational principle carries a matrix-product state +//! forward with its error held by the bond it keeps; until the entanglement +//! outgrows that bond, the dynamics are followed to rounding. This module is +//! that method, deterministic and held to exact answers. +//! +//! - **Two-site TDVP.** [`Evolution`] sweeps a matrix-product state left and +//! right once a step, a symmetric second-order scheme. +//! - Each bond's two-site tensor goes forward half a step under its +//! effective Hamiltonian; the new centre then goes back half a step under +//! its one-site one. +//! - The exponentials are Lanczos with full re-orthogonalisation +//! ([`expm_krylov`]). A step that does not converge in `krylov` vectors +//! is halved. +//! - Bonds are split by a singular value decomposition (one-sided Jacobi), +//! so amplitudes are resolved to rounding, not to its square root as a +//! density-matrix split would resolve them. +//! - **Complete bases.** TDVP's projection error comes from bases too small +//! to express `H|ψ⟩`, and a product state's bases are as small as they +//! come. [`Evolution::new`] completes every bond's basis with orthonormal +//! null vectors up to `max_bond`: the state is unchanged and the tangent +//! space grows. With `cutoff = 0` the bases stay complete, and a chain that +//! fits within `max_bond` then evolves with every projector the identity: +//! exact to rounding at any step size. +//! - **Measurements.** Site and bond expectations, the entanglement entropy +//! at every bond, the energy, and the overlap with any other state, which +//! gives the Loschmidt echo. +//! - **Exact referees.** +//! - [`exact_evolve`] evolves all `2ⁿ` amplitudes by Krylov steps (at most +//! 20 sites). +//! - [`ising_quench`] solves the transverse-field Ising chain started with +//! every `σˣ = +1` as free fermions, at any length: `σˣ` on every site, +//! `σᶻσᶻ` on every bond, and the Loschmidt echo as the overlap of two +//! Gaussian states. +//! +//! The engine is the DMRG module's, on complex numbers. Every sum runs in a +//! fixed order and the effective Hamiltonians split across threads by rows, +//! so a run is the same bits on every machine and at any thread count. +//! +//! # Checked +//! +//! - **The referees agree.** The free-fermion solution matches exact +//! evolution to `10⁻¹⁰` on 8 sites for three quenches, echo included. The +//! Krylov exponential matches the dense one to `10⁻¹¹`, including a long +//! step it has to halve. +//! - **Exact when the basis is complete.** On 8 sites with complete bases, +//! four steps of 0.5 match exact evolution to `10⁻¹¹` in the state. Two +//! cases are checked: Heisenberg from the Néel state, and Ising from a +//! tilted product state. +//! - **Truncated, against free fermions.** 16 sites at bond 8 follow the +//! critical Ising quench to `10⁻⁶` in `σˣ`, `σᶻσᶻ` and the echo. +//! - **At scale** (`tdvp_quench` example; the largest error in `σˣ` over +//! every site, against the free-fermion answer): +//! +//! | sites | bond | t | error in σˣ | discarded weight | entropy, mid | +//! |-------|------|-----|-------------|------------------|--------------| +//! | 100 | 64 | 1 | 5.6·10⁻¹¹ | 1·10⁻¹³ | 0.74 | +//! | 100 | 64 | 2 | 2.5·10⁻¹⁰ | 3·10⁻¹³ | 1.39 | +//! | 100 | 64 | 3 | 9.0·10⁻¹⁰ | 4·10⁻¹⁰ | 2.02 | +//! | 100 | 64 | 3.5 | 6.7·10⁻⁷ | 8·10⁻⁷ | 2.34 | +//! | 100 | 64 | 4 | 1.6·10⁻⁴ | 1.5·10⁻⁴ | 2.66 | +//! | 100 | 32 | 2.5 | 1.8·10⁻⁸ | 9·10⁻⁹ | 1.70 | +//! | 100 | 32 | 3 | 1.5·10⁻⁵ | 8·10⁻⁶ | 2.02 | +//! | 100 | 32 | 3.5 | 1.9·10⁻³ | 1.7·10⁻³ | 2.34 | +//! +//! The critical quench (`J = h = 1`) at step 0.1. Until the entanglement +//! fills the bond, the error stays below `10⁻⁹` and the Loschmidt rate +//! agrees with the exact one to seven digits. Once it does, the error +//! follows the discarded weight up, within a factor of three, which is why +//! that weight is reported. The bond buys time: at `t = 3` the error is +//! `1.5·10⁻⁵` at bond 32 and `9·10⁻¹⁰` at bond 64. +//! +//! - **The tests bite.** Lengthening the backward half-step by 0.1 % fails +//! both accuracy tests. +//! - **Conserved and repeatable.** Energy holds to `10⁻⁹` over ten truncated +//! steps. Runs repeat bit for bit at any thread count, and the same bits +//! come out natively and under `wasm32-wasip2`. +//! +//! # Measured, and not kept +//! +//! A fourth-order composition of the symmetric step (Yoshida's three-stage +//! one) was built and measured on the 40-site critical quench at bond 32. At +//! time steps 0.05, 0.1 and 0.2 it was less accurate than the second-order +//! step: `1.1·10⁻⁹` against `2.3·10⁻¹⁰` at `t = 2` for step 0.05. It was also +//! three to ten times slower. Its backward middle stage entangles the state +//! the next split must truncate, and the time-step error it removes is not +//! where TDVP's error lies. +//! +//! # Honest boundaries +//! +//! - **The same Hamiltonians as DMRG**: one dimension, nearest neighbours, +//! real terms. +//! - **Entanglement sets the reach.** After a quench it grows linearly in +//! time, so the bond needed grows exponentially. Past that point a result +//! carries truncation error. The discarded weight is reported; it is not +//! converted into a bound on observables. +//! - **The cutoff trades accuracy for speed.** A cutoff drops basis vectors +//! the state does not yet use, and with them some of the tangent space. +//! On the 40-site critical quench at bond 32 and step 0.05, the error in +//! `σˣ` at `t = 2` is `4·10⁻⁶` at cutoff `10⁻¹²`, `3·10⁻⁸` at `10⁻¹⁶` and +//! `2·10⁻¹⁰` at `10⁻²⁰`, the default. At zero it is `4·10⁻¹⁰`: by then the +//! bond is full and truncation sets the error whatever the cutoff, at three +//! times the cost. + +use crate::linalg::jacobi_svd_strict; +use crate::quantum_dmrg::{C, Chain, Engine, Mpo, Mps, Op, Scalar, apply_dense, default_threads, entropy, hermitian_eigen, right_canonicalise, symmetric_eigen, transpose, two_site}; +use crate::repro::sin_cos; + +// --------------------------------------------------------------------------- +// Dense helpers +// --------------------------------------------------------------------------- + +fn dotc(a: &[C], b: &[C]) -> C { + let mut s = C::ZERO; + for (x, y) in a.iter().zip(b) { + s = s.add(x.conj().mul(*y)); + } + s +} + +fn norm(v: &[C]) -> f64 { + v.iter().map(|x| x.norm2()).sum::().sqrt() +} + +// --------------------------------------------------------------------------- +// The Krylov exponential +// --------------------------------------------------------------------------- + +/// `exp(−i τ H) v` for a Hermitian `H` given by `apply`: Lanczos with full +/// re-orthogonalisation, then the exponential of the small tridiagonal +/// matrix. It stops once the next Krylov vector would carry less than `tol` +/// of the result; if `krylov` vectors are not enough, the step is halved and +/// taken twice. +pub fn expm_krylov(apply: &dyn Fn(&[C], &mut [C]), v: &[C], tau: f64, krylov: usize, tol: f64) -> Vec { + expm_halving(apply, v, tau, krylov, tol, 0) +} + +fn expm_halving(apply: &dyn Fn(&[C], &mut [C]), v: &[C], tau: f64, krylov: usize, tol: f64, depth: u32) -> Vec { + match expm_once(apply, v, tau, krylov, tol) { + Some(out) => out, + None => { + // Each halving shrinks the error estimate by at least half once + // two vectors are allowed, so this depth means the input is not + // finite or the operator not Hermitian. + assert!(depth < 60, "the Krylov exponential did not converge: is the state finite and the operator Hermitian?"); + let half = expm_halving(apply, v, tau / 2.0, krylov, tol, depth + 1); + expm_halving(apply, &half, tau / 2.0, krylov, tol, depth + 1) + } + } +} + +fn expm_once(apply: &dyn Fn(&[C], &mut [C]), v: &[C], tau: f64, krylov: usize, tol: f64) -> Option> { + let dim = v.len(); + let beta0 = norm(v); + if beta0 == 0.0 || tau == 0.0 { + return Some(v.to_vec()); + } + // Two vectors at least: with one, the error estimate does not fall as + // the step does. + let k = krylov.max(2).min(dim); + let mut basis: Vec> = vec![v.iter().map(|x| x.scale(1.0 / beta0)).collect()]; + let (mut alpha, mut beta): (Vec, Vec) = (Vec::new(), Vec::new()); + let mut w = vec![C::ZERO; dim]; + loop { + let j = basis.len() - 1; + apply(&basis[j], &mut w); + alpha.push(dotc(&basis[j], &w).re); + // Full re-orthogonalisation, twice. + for _ in 0..2 { + for b in &basis { + let c = dotc(b, &w); + for (x, y) in w.iter_mut().zip(b) { + *x = x.sub(c.mul(*y)); + } + } + } + let bn = norm(&w); + let m = alpha.len(); + let coeffs = exp_tridiagonal(&alpha, &beta, tau); + let err = bn * coeffs[m - 1].norm2().sqrt(); + if err <= tol || m == dim { + let mut out = vec![C::ZERO; dim]; + for (c, b) in coeffs.iter().zip(&basis) { + let c = c.scale(beta0); + for (o, x) in out.iter_mut().zip(b) { + *o = o.add(c.mul(*x)); + } + } + return Some(out); + } + if m == k { + return None; + } + beta.push(bn); + basis.push(w.iter().map(|x| x.scale(1.0 / bn)).collect()); + } +} + +/// `exp(−i τ T) e₁` for the real symmetric tridiagonal `T` with diagonal +/// `alpha` and off-diagonal `beta`. +fn exp_tridiagonal(alpha: &[f64], beta: &[f64], tau: f64) -> Vec { + let m = alpha.len(); + let mut t = vec![0.0; m * m]; + for i in 0..m { + t[i * m + i] = alpha[i]; + if i + 1 < m { + t[i * m + i + 1] = beta[i]; + t[(i + 1) * m + i] = beta[i]; + } + } + let (vals, vecs) = symmetric_eigen(&t, m); + let phases: Vec = vals + .iter() + .map(|&l| { + let (s, c) = sin_cos(tau * l); + C { re: c, im: -s } + }) + .collect(); + (0..m) + .map(|r| { + let mut s = C::ZERO; + for (l, p) in phases.iter().enumerate() { + s = s.add(p.scale(vecs[r * m + l] * vecs[l])); + } + s + }) + .collect() +} + +// --------------------------------------------------------------------------- +// TDVP +// --------------------------------------------------------------------------- + +/// Settings for [`Evolution`]. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct TdvpConfig { + /// The time step. + pub dt: f64, + pub max_bond: usize, + /// Discard Schmidt weight below this at each split. Zero keeps every + /// basis vector up to `max_bond`, null ones included, so the tangent + /// space is as large as the bond allows. + pub cutoff: f64, + /// Lanczos vectors per exponential before its step is halved. + pub krylov: usize, + /// The error allowed in each exponential, relative to the state. + pub tol: f64, + /// Threads for the effective Hamiltonians; the result does not depend + /// on it. + pub threads: usize, +} + +impl Default for TdvpConfig { + fn default() -> TdvpConfig { + TdvpConfig { dt: 0.05, max_bond: 64, cutoff: 1e-20, krylov: 30, tol: 1e-12, threads: default_threads() } + } +} + +/// What one step did. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct StepReport { + /// Density-matrix weight discarded, summed over every split of the step. + pub discarded: f64, + /// The largest bond after the step. + pub max_bond: usize, +} + +/// A matrix-product state evolving under a chain's Hamiltonian by the +/// two-site time-dependent variational principle. Between steps the state is +/// right-canonical with its centre on site 0. +#[derive(Clone, Debug, PartialEq)] +pub struct Evolution { + mpo: Mpo, + mps: Mps, + /// `rights[i]` is the environment of sites `i..n`. + rights: Vec>, + cfg: TdvpConfig, + steps: u64, +} + +impl Evolution { + /// Start from `state` (normalised here) under `chain`. + pub fn new(chain: &Chain, state: &Mps, cfg: TdvpConfig) -> Evolution { + let n = chain.n; + assert!(n >= 2, "TDVP needs at least two sites"); + assert_eq!(state.n(), n, "the state and the chain differ in length"); + let mpo = chain.mpo(); + let mut mps = state.clone(); + right_canonicalise(&mut mps); + let nrm = norm(&mps.sites[0]); + assert!(nrm > 0.0, "the state is zero"); + mps.sites[0].iter_mut().for_each(|x| *x = x.scale(1.0 / nrm)); + complete_bases(&mut mps, cfg.max_bond); + let eng = Engine::new(&mpo); + let mut rights: Vec> = vec![Vec::new(); n + 1]; + rights[n] = eng.right_edge(); + for i in (1..n).rev() { + rights[i] = eng.grow_right(&rights[i + 1], &mps.sites[i], mps.dims[i], mps.dims[i + 1]); + } + Evolution { mpo, mps, rights, cfg, steps: 0 } + } + + /// Advance by one time step. + pub fn step(&mut self) -> StepReport { + let discarded = self.sweep_pair(self.cfg.dt); + self.steps += 1; + StepReport { discarded, max_bond: self.max_bond() } + } + + /// One symmetric step of length `tau`: left to right and back, each + /// bond forward by `tau / 2` and each centre in between back by as much. + /// Returns the weight discarded. + fn sweep_pair(&mut self, tau: f64) -> f64 { + let n = self.mps.n(); + let cfg = self.cfg; + let half = tau / 2.0; + let eng = Engine::new(&self.mpo); + let mps = &mut self.mps; + let rights = &mut self.rights; + let mut lefts: Vec> = vec![Vec::new(); n + 1]; + lefts[0] = eng.left_edge(); + let mut discarded = 0.0; + for i in 0..n - 1 { + let (dl, dm, dr) = (mps.dims[i], mps.dims[i + 1], mps.dims[i + 2]); + let theta = two_site(&mps.sites[i], &mps.sites[i + 1], dl, dm, dr); + let (l, r) = (&lefts[i], &rights[i + 2]); + let theta = expm_krylov(&|x: &[C], y: &mut [C]| eng.apply(l, r, x, y, dl, dr, cfg.threads), &theta, half, cfg.krylov, cfg.tol); + let (u, mut rest, disc, weights) = svd_split(&theta, dl * 2, 2 * dr, cfg.max_bond, cfg.cutoff); + if disc > 0.0 { + renormalise(&mut rest, norm(&theta)); + } + discarded += disc; + let k = weights.len(); + mps.sites[i] = u; + mps.dims[i + 1] = k; + lefts[i + 1] = eng.grow_left(&lefts[i], &mps.sites[i], dl, k); + mps.sites[i + 1] = if i + 2 < n { + let (l, r) = (&lefts[i + 1], &rights[i + 2]); + expm_krylov(&|x: &[C], y: &mut [C]| eng.apply1(l, r, x, y, k, dr, cfg.threads), &rest, -half, cfg.krylov, cfg.tol) + } else { + rest + }; + } + for i in (0..n - 1).rev() { + let (dl, dm, dr) = (mps.dims[i], mps.dims[i + 1], mps.dims[i + 2]); + let theta = two_site(&mps.sites[i], &mps.sites[i + 1], dl, dm, dr); + let (l, r) = (&lefts[i], &rights[i + 2]); + let theta = expm_krylov(&|x: &[C], y: &mut [C]| eng.apply(l, r, x, y, dl, dr, cfg.threads), &theta, half, cfg.krylov, cfg.tol); + let (u, rest, disc, weights) = svd_split(&transpose(&theta, dl * 2, 2 * dr), 2 * dr, dl * 2, cfg.max_bond, cfg.cutoff); + discarded += disc; + let k = weights.len(); + mps.sites[i + 1] = transpose(&u, 2 * dr, k); + mps.dims[i + 1] = k; + rights[i + 1] = eng.grow_right(&rights[i + 2], &mps.sites[i + 1], k, dr); + let mut centre = transpose(&rest, k, dl * 2); + if disc > 0.0 { + renormalise(&mut centre, norm(&theta)); + } + mps.sites[i] = if i > 0 { + let (l, r) = (&lefts[i], &rights[i + 1]); + expm_krylov(&|x: &[C], y: &mut [C]| eng.apply1(l, r, x, y, dl, k, cfg.threads), ¢re, -half, cfg.krylov, cfg.tol) + } else { + centre + }; + } + discarded + } + + /// The time evolved so far. + pub fn time(&self) -> f64 { + self.steps as f64 * self.cfg.dt + } + + pub fn state(&self) -> &Mps { + &self.mps + } + + pub fn max_bond(&self) -> usize { + self.mps.dims.iter().copied().max().unwrap_or(1) + } + + /// `‖ψ‖`: one, less whatever truncation has not put back. + pub fn norm(&self) -> f64 { + norm(&self.mps.sites[0]) + } + + /// `⟨ψ|H|ψ⟩`, from the centre and its environments. + pub fn energy(&self) -> f64 { + let eng = Engine::new(&self.mpo); + let c = &self.mps.sites[0]; + let (dl, dr) = (self.mps.dims[0], self.mps.dims[1]); + let mut hc = vec![C::ZERO; c.len()]; + eng.apply1(&eng.left_edge(), &self.rights[1], c, &mut hc, dl, dr, self.cfg.threads); + dotc(c, &hc).re / dotc(c, c).re + } + + /// `⟨O_i⟩` on every site. For a real symmetric `op` this is its + /// expectation; otherwise its real part. + pub fn expect(&self, op: Op) -> Vec { + let mut out = Vec::with_capacity(self.mps.n()); + let mut env = vec![C::ONE]; + for i in 0..self.mps.n() { + let (dl, dr) = (self.mps.dims[i], self.mps.dims[i + 1]); + let m = &self.mps.sites[i]; + let t = env_times(&env, m, dl, 2 * dr); + let mut v = C::ZERO; + for ap in 0..dl { + for sp in 0..2 { + for s in 0..2 { + let o = op.0[sp * 2 + s]; + if o == 0.0 { + continue; + } + for b in 0..dr { + v = v.add(m[(ap * 2 + sp) * dr + b].conj().mul(t[(ap * 2 + s) * dr + b]).scale(o)); + } + } + } + } + out.push(v.re); + env = env_close(m, &t, dl, dr); + } + out + } + + /// `⟨A_i B_{i+1}⟩` on every bond, real parts as for [`Evolution::expect`]. + pub fn bond_expect(&self, a: Op, b: Op) -> Vec { + let n = self.mps.n(); + let mut out = Vec::with_capacity(n - 1); + let mut env = vec![C::ONE]; + for i in 0..n - 1 { + let (dl, dm, dr) = (self.mps.dims[i], self.mps.dims[i + 1], self.mps.dims[i + 2]); + let theta = two_site(&self.mps.sites[i], &self.mps.sites[i + 1], dl, dm, dr); + let t = env_times(&env, &theta, dl, 4 * dr); + let mut v = C::ZERO; + for ap in 0..dl { + for s1p in 0..2 { + for s2p in 0..2 { + for s1 in 0..2 { + let oa = a.0[s1p * 2 + s1]; + if oa == 0.0 { + continue; + } + for s2 in 0..2 { + let ob = b.0[s2p * 2 + s2]; + if ob == 0.0 { + continue; + } + for c in 0..dr { + let bra = theta[((ap * 2 + s1p) * 2 + s2p) * dr + c].conj(); + v = v.add(bra.mul(t[((ap * 2 + s1) * 2 + s2) * dr + c]).scale(oa * ob)); + } + } + } + } + } + } + out.push(v.re); + let m = &self.mps.sites[i]; + let t1 = env_times(&env, m, dl, 2 * dm); + env = env_close(m, &t1, dl, dm); + } + out + } + + /// Von Neumann entanglement entropy at each of the `n − 1` bonds: the + /// state's Schmidt weights there are the eigenvalues of the Gram matrix + /// of its left part, since its right part is orthonormal. + pub fn entropies(&self) -> Vec { + let n = self.mps.n(); + let mut out = Vec::with_capacity(n - 1); + let mut env = vec![C::ONE]; + for i in 0..n - 1 { + let (dl, dr) = (self.mps.dims[i], self.mps.dims[i + 1]); + let m = &self.mps.sites[i]; + let t = env_times(&env, m, dl, 2 * dr); + env = env_close(m, &t, dl, dr); + let (vals, _) = hermitian_eigen(&env, dr); + let total: f64 = vals.iter().map(|v| v.max(0.0)).sum(); + let weights: Vec = vals.iter().map(|v| v.max(0.0) / total).collect(); + out.push(entropy(&weights)); + } + out + } + + /// `⟨φ|ψ⟩` with another state of the same length. + pub fn overlap(&self, phi: &Mps) -> C { + overlap(phi, &self.mps) + } +} + +/// Complete the orthonormal rows of every right-canonical site, from the +/// right, up to the bond the chain allows or `max_bond`, and give the site to +/// the left zero columns for them. The state does not change; the tangent +/// space it spans grows, so a state with small bonds (a product state, say) +/// is not held back by bases that cannot yet express `H|ψ⟩`. +fn complete_bases(mps: &mut Mps, max_bond: usize) { + let n = mps.n(); + for i in (1..n).rev() { + let (k, dr) = (mps.dims[i], mps.dims[i + 1]); + let cols = 2 * dr; + let left_space = 1usize.checked_shl(i.min(40) as u32).unwrap_or(usize::MAX); + let target = cols.min(left_space).min(max_bond.max(k)); + if target <= k { + continue; + } + let mut rows: Vec> = mps.sites[i].chunks(cols).map(|r| r.to_vec()).collect(); + for unit in 0..cols { + if rows.len() == target { + break; + } + let mut v = vec![C::ZERO; cols]; + v[unit] = C::ONE; + // Gram–Schmidt against the rows so far, twice. + for _ in 0..2 { + for r in &rows { + let c = dotc(r, &v); + for (x, y) in v.iter_mut().zip(r) { + *x = x.sub(c.mul(*y)); + } + } + } + let nv = norm(&v); + if nv > 0.5 { + rows.push(v.iter().map(|x| x.scale(1.0 / nv)).collect()); + } + } + let added = rows.len(); + mps.sites[i] = rows.into_iter().flatten().collect(); + // The site to the left gains zero columns. + let pl = mps.dims[i - 1]; + let prev = &mps.sites[i - 1]; + let mut wider = vec![C::ZERO; pl * 2 * added]; + for r in 0..pl * 2 { + wider[r * added..r * added + k].copy_from_slice(&prev[r * k..(r + 1) * k]); + } + mps.sites[i - 1] = wider; + mps.dims[i] = added; + } +} + +/// Split `m` (rows × cols) by its singular value decomposition, keeping at +/// most `max_bond` (and never more than its rank allows) and dropping weight +/// below `cutoff`. Returns the kept left vectors (rows × k), `u† m` +/// (k × cols), the discarded weight and the kept normalised weights. +/// +/// One-sided Jacobi runs on `m†`, so the left vectors are its accumulated +/// rotations: orthonormal to rounding however small the singular value, and +/// a complete basis to draw null vectors from when the cutoff is zero. +/// Amplitudes are resolved to rounding, not to its square root as a +/// density-matrix split would. +fn svd_split(m: &[C], rows: usize, cols: usize, max_bond: usize, cutoff: f64) -> (Vec, Vec, f64, Vec) { + // Columns of m† are the conjugated rows of m. + let adj: Vec> = (0..rows).map(|r| m[r * cols..(r + 1) * cols].iter().map(|x| x.conj()).collect()).collect(); + let (_, s, v) = jacobi_svd_strict(adj, cols, rows); + let total: f64 = s.iter().map(|x| x * x).sum(); + let cap = max_bond.min(rows).min(cols).max(1); + let mut k = 0; + let mut kept = 0.0; + for &sv in &s { + let w = if total > 0.0 { sv * sv / total } else { 0.0 }; + if k >= cap || (k > 0 && w < cutoff) { + break; + } + kept += w; + k += 1; + } + let k = k.max(1); + let mut u = vec![C::ZERO; rows * k]; + for (c, col) in v.iter().take(k).enumerate() { + for r in 0..rows { + u[r * k + c] = col[r]; + } + } + let mut rest = vec![C::ZERO; k * cols]; + for c in 0..k { + for r in 0..rows { + let uv = u[r * k + c]; + if uv.is_zero() { + continue; + } + let uv = uv.conj(); + for j in 0..cols { + rest[c * cols + j] = rest[c * cols + j].add(uv.mul(m[r * cols + j])); + } + } + } + let weights = s.iter().take(k).map(|x| if total > 0.0 { x * x / total } else { 0.0 }).collect(); + (u, rest, (1.0 - kept).max(0.0), weights) +} + +fn renormalise(v: &mut [C], target: f64) { + let nv = norm(v); + if nv > 0.0 { + v.iter_mut().for_each(|x| *x = x.scale(target / nv)); + } +} + +/// `t[a', x] = Σ_a env[a', a] m[a, x]` for `m` with `dl` rows of `cols`. +fn env_times(env: &[C], m: &[C], dl: usize, cols: usize) -> Vec { + let mut t = vec![C::ZERO; dl * cols]; + for ap in 0..dl { + let dst = &mut t[ap * cols..(ap + 1) * cols]; + for a in 0..dl { + let e = env[ap * dl + a]; + if e.is_zero() { + continue; + } + for (d, &x) in dst.iter_mut().zip(&m[a * cols..(a + 1) * cols]) { + *d = d.add(e.mul(x)); + } + } + } + t +} + +/// `env'[b', b] = Σ_{a', s} m*[a', s, b'] t[a', s, b]`. +fn env_close(m: &[C], t: &[C], dl: usize, dr: usize) -> Vec { + let mut out = vec![C::ZERO; dr * dr]; + for ap in 0..dl { + for s in 0..2 { + for bp in 0..dr { + let bra = m[(ap * 2 + s) * dr + bp]; + if bra.is_zero() { + continue; + } + let bra = bra.conj(); + let dst = &mut out[bp * dr..(bp + 1) * dr]; + for (d, &x) in dst.iter_mut().zip(&t[(ap * 2 + s) * dr..(ap * 2 + s + 1) * dr]) { + *d = d.add(bra.mul(x)); + } + } + } + } + out +} + +/// `⟨φ|ψ⟩` for two matrix-product states of the same length. +pub fn overlap(phi: &Mps, psi: &Mps) -> C { + assert_eq!(phi.n(), psi.n(), "the states differ in length"); + // env[a', a]: a' indexes φ's bond, a indexes ψ's. + let mut env = vec![C::ONE]; + for i in 0..psi.n() { + let (pl, pr) = (phi.dims[i], phi.dims[i + 1]); + let (dl, dr) = (psi.dims[i], psi.dims[i + 1]); + // t[a', s, b] = Σ_a env[a', a] ψ[a, s, b] + let mut t = vec![C::ZERO; pl * 2 * dr]; + for ap in 0..pl { + for a in 0..dl { + let e = env[ap * dl + a]; + if e.is_zero() { + continue; + } + for (d, &x) in t[ap * 2 * dr..(ap + 1) * 2 * dr].iter_mut().zip(&psi.sites[i][a * 2 * dr..(a + 1) * 2 * dr]) { + *d = d.add(e.mul(x)); + } + } + } + let mut next = vec![C::ZERO; pr * dr]; + for ap in 0..pl { + for s in 0..2 { + for bp in 0..pr { + let bra = phi.sites[i][(ap * 2 + s) * pr + bp]; + if bra.is_zero() { + continue; + } + let bra = bra.conj(); + for (d, &x) in next[bp * dr..(bp + 1) * dr].iter_mut().zip(&t[(ap * 2 + s) * dr..(ap * 2 + s + 1) * dr]) { + *d = d.add(bra.mul(x)); + } + } + } + } + env = next; + } + env[0] +} + +// --------------------------------------------------------------------------- +// Exact referees +// --------------------------------------------------------------------------- + +/// `exp(−iHt) ψ` for `chain` on all `2ⁿ` amplitudes (at most 20 sites), in +/// `steps` Krylov steps. +pub fn exact_evolve(chain: &Chain, psi: &[C], t: f64, steps: usize) -> Vec { + let n = chain.n; + assert!((1..=20).contains(&n), "exact evolution takes 1 to 20 sites"); + assert_eq!(psi.len(), 1 << n); + let apply = |x: &[C], y: &mut [C]| { + let re: Vec = x.iter().map(|c| c.re).collect(); + let im: Vec = x.iter().map(|c| c.im).collect(); + let (mut hre, mut him) = (vec![0.0; x.len()], vec![0.0; x.len()]); + apply_dense(chain, &re, &mut hre); + apply_dense(chain, &im, &mut him); + for ((o, r), i) in y.iter_mut().zip(hre).zip(him) { + *o = C { re: r, im: i }; + } + }; + let steps = steps.max(1); + let mut v = psi.to_vec(); + for _ in 0..steps { + v = expm_krylov(&apply, &v, t / steps as f64, 40, 1e-14); + } + v +} + +/// `⟨ψ|O_site|ψ⟩` on a dense state of `n` sites (real part). +pub fn dense_expect(psi: &[C], n: usize, site: usize, op: Op) -> f64 { + let shift = n - 1 - site; + let mut v = C::ZERO; + for (x, &) in psi.iter().enumerate() { + let s = (x >> shift) & 1; + for sp in 0..2 { + let o = op.0[sp * 2 + s]; + if o != 0.0 { + let y = (x & !(1 << shift)) | (sp << shift); + v = v.add(psi[y].conj().mul(amp).scale(o)); + } + } + } + v.re +} + +/// `⟨ψ|A_site B_{site+1}|ψ⟩` on a dense state of `n` sites (real part). +pub fn dense_bond_expect(psi: &[C], n: usize, site: usize, a: Op, b: Op) -> f64 { + let (sa, sb) = (n - 1 - site, n - 2 - site); + let mut v = C::ZERO; + for (x, &) in psi.iter().enumerate() { + let (s1, s2) = ((x >> sa) & 1, (x >> sb) & 1); + for s1p in 0..2 { + let oa = a.0[s1p * 2 + s1]; + if oa == 0.0 { + continue; + } + for s2p in 0..2 { + let ob = b.0[s2p * 2 + s2]; + if ob == 0.0 { + continue; + } + let y = (x & !((1 << sa) | (1 << sb))) | (s1p << sa) | (s2p << sb); + v = v.add(psi[y].conj().mul(amp).scale(oa * ob)); + } + } + } + v.re +} + +/// The transverse-field Ising chain after a quench, solved exactly. +#[derive(Clone, Debug, PartialEq)] +pub struct IsingQuench { + /// `⟨σˣ_i(t)⟩` on every site. + pub x: Vec, + /// `⟨σᶻ_i σᶻ_{i+1}(t)⟩` on every bond. + pub zz: Vec, + /// The Loschmidt echo `|⟨ψ(0)|ψ(t)⟩|²`. + pub echo: f64, +} + +/// The open chain `H = −J Σ σᶻσᶻ − h Σ σˣ`, started with every `σˣ = +1` +/// and evolved for time `t`, as free fermions. +/// +/// With Majorana operators `γ_{2i} = (Π_{k IsingQuench { + assert!(n >= 1, "the chain needs a site"); + let m = 2 * n; + let mut a = vec![0.0; m * m]; + for i in 0..n { + a[(2 * i) * m + 2 * i + 1] = -2.0 * h * t; + a[(2 * i + 1) * m + 2 * i] = 2.0 * h * t; + if i + 1 < n { + a[(2 * i + 1) * m + 2 * i + 2] = -2.0 * j * t; + a[(2 * i + 2) * m + 2 * i + 1] = 2.0 * j * t; + } + } + let r = expm_real(&a, m); + // Γ(t)_ab = Σ_k (R_{a,2k} R_{b,2k+1} − R_{a,2k+1} R_{b,2k}). + let gamma = |p: usize, q: usize| -> f64 { + let mut s = 0.0; + for k in 0..n { + s += r[p * m + 2 * k] * r[q * m + 2 * k + 1] - r[p * m + 2 * k + 1] * r[q * m + 2 * k]; + } + s + }; + let x = (0..n).map(|i| gamma(2 * i, 2 * i + 1)).collect(); + let zz = (0..n.saturating_sub(1)).map(|i| gamma(2 * i + 1, 2 * i + 2)).collect(); + // (Γ(0) + Γ(t)) / 2. + let mut s = vec![0.0; m * m]; + for p in 0..m { + for q in 0..m { + if p != q { + s[p * m + q] = 0.5 * gamma(p, q); + } + } + } + for i in 0..n { + s[(2 * i) * m + 2 * i + 1] += 0.5; + s[(2 * i + 1) * m + 2 * i] -= 0.5; + } + let echo = determinant(s, m).abs().sqrt(); + IsingQuench { x, zz, echo } +} + +/// `e^A` for a real `n × n` matrix: Taylor's series to degree 18 after +/// halving `A` until its 1-norm is at most ½, then squaring back. +fn expm_real(a: &[f64], n: usize) -> Vec { + let norm1 = (0..n).map(|c| (0..n).map(|r| a[r * n + c].abs()).sum::()).fold(0.0, f64::max); + let mut halvings = 0; + let mut scale = 1.0; + while norm1 * scale > 0.5 { + scale *= 0.5; + halvings += 1; + } + let b: Vec = a.iter().map(|x| x * scale).collect(); + let identity = |m: &mut [f64]| { + for i in 0..n { + m[i * n + i] += 1.0; + } + }; + // Horner: I + B (I + B/2 (I + … (I + B/18))). + let mut e = vec![0.0; n * n]; + identity(&mut e); + for k in (1..=18).rev() { + let mut next = matmul(&b, &e, n); + next.iter_mut().for_each(|x| *x /= k as f64); + identity(&mut next); + e = next; + } + for _ in 0..halvings { + e = matmul(&e, &e, n); + } + e +} + +fn matmul(a: &[f64], b: &[f64], n: usize) -> Vec { + let mut c = vec![0.0; n * n]; + for i in 0..n { + for k in 0..n { + let av = a[i * n + k]; + if av == 0.0 { + continue; + } + for (d, &x) in c[i * n..(i + 1) * n].iter_mut().zip(&b[k * n..(k + 1) * n]) { + *d += av * x; + } + } + } + c +} + +/// The determinant by Gaussian elimination with partial pivoting. +fn determinant(mut a: Vec, n: usize) -> f64 { + let mut det = 1.0; + for c in 0..n { + let mut p = c; + for r in c + 1..n { + if a[r * n + c].abs() > a[p * n + c].abs() { + p = r; + } + } + if a[p * n + c] == 0.0 { + return 0.0; + } + if p != c { + for k in 0..n { + a.swap(c * n + k, p * n + k); + } + det = -det; + } + let piv = a[c * n + c]; + det *= piv; + for r in c + 1..n { + let f = a[r * n + c] / piv; + if f == 0.0 { + continue; + } + for k in c..n { + a[r * n + k] -= f * a[c * n + k]; + } + } + } + det +} + +#[cfg(test)] +mod tests { + use super::*; + + fn rng(seed: u64) -> impl FnMut() -> f64 { + let mut s = seed; + move || { + s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + ((s >> 11) as f64 / 9_007_199_254_740_992.0) - 0.5 + } + } + + #[test] + fn the_krylov_exponential_matches_the_dense_one() { + let mut next = rng(5); + let n = 24; + let mut h = vec![C::ZERO; n * n]; + for i in 0..n { + h[i * n + i] = C { re: 3.0 * next(), im: 0.0 }; + for j in i + 1..n { + let v = C { re: next(), im: next() }; + h[i * n + j] = v; + h[j * n + i] = v.conj(); + } + } + let v: Vec = (0..n).map(|_| C { re: next(), im: next() }).collect(); + let (vals, vecs) = hermitian_eigen(&h, n); + let apply = |x: &[C], y: &mut [C]| { + for r in 0..n { + let mut s = C::ZERO; + for c in 0..n { + s = s.add(h[r * n + c].mul(x[c])); + } + y[r] = s; + } + }; + // Short steps converge in a few vectors; a long one with few vectors + // has to halve itself. + for (tau, krylov) in [(0.1, 30), (-0.7, 30), (6.0, 8)] { + let got = expm_krylov(&apply, &v, tau, krylov, 1e-13); + // Dense: X e^{−iτΛ} X† v. + for r in 0..n { + let mut want = C::ZERO; + for k in 0..n { + let mut proj = C::ZERO; + for c in 0..n { + proj = proj.add(vecs[c * n + k].conj().mul(v[c])); + } + let (s, co) = sin_cos(tau * vals[k]); + want = want.add(vecs[r * n + k].mul(C { re: co, im: -s }).mul(proj)); + } + assert!(got[r].sub(want).norm2().sqrt() < 1e-11, "tau={tau} row {r}"); + } + } + } + + #[test] + fn free_fermions_solve_the_quench() { + let n = 8; + let plus = Mps::product(&vec![[C::ONE, C::ONE]; n]).to_dense().unwrap(); + for (j, h, t) in [(1.0, 1.0, 0.7), (1.0, 0.4, 2.0), (0.6, 1.3, 1.1)] { + let chain = Chain::ising(n, j, h); + let psi = exact_evolve(&chain, &plus, t, 8); + let ff = ising_quench(n, j, h, t); + for i in 0..n { + let x = dense_expect(&psi, n, i, Op::x()); + assert!((x - ff.x[i]).abs() < 1e-10, "J={j} h={h} t={t} x[{i}]: {x} vs {}", ff.x[i]); + } + for i in 0..n - 1 { + let zz = dense_bond_expect(&psi, n, i, Op::z(), Op::z()); + assert!((zz - ff.zz[i]).abs() < 1e-10, "zz[{i}]: {zz} vs {}", ff.zz[i]); + } + let echo = dotc(&plus, &psi).norm2() / dotc(&plus, &plus).re.powi(2); + assert!((echo - ff.echo).abs() < 1e-10, "echo: {echo} vs {}", ff.echo); + } + } + + #[test] + fn tdvp_is_exact_when_the_basis_is_complete() { + // Completed bases make every projector the identity: no projection + // error and no splitting error, at any step, from product states. + let n = 8; + let neel: Vec = (0..n).map(|i| i % 2 == 1).collect(); + let mut next = rng(9); + let tilted: Vec<[C; 2]> = (0..n).map(|_| [C { re: 1.0, im: 0.0 }, C { re: next(), im: next() }]).collect(); + let cases = [(Chain::heisenberg(n, 1.0, 0.6, 0.2), Mps::basis(&neel)), (Chain::ising(n, 1.0, 0.7), Mps::product(&tilted))]; + for (chain, start) in cases { + let cfg = TdvpConfig { dt: 0.5, max_bond: 16, cutoff: 0.0, ..TdvpConfig::default() }; + let mut ev = Evolution::new(&chain, &start, cfg); + for _ in 0..4 { + ev.step(); + } + let exact = exact_evolve(&chain, &start.to_dense().unwrap(), ev.time(), 8); + let got = ev.state().to_dense().unwrap(); + let err = got.iter().zip(&exact).map(|(a, b)| a.sub(*b).norm2()).sum::().sqrt(); + assert!(err < 1e-11, "{chain:?}: state error {err:e}"); + let z = ev.expect(Op::z()); + for (i, zi) in z.iter().enumerate() { + assert!((zi - dense_expect(&exact, n, i, Op::z())).abs() < 1e-11); + } + } + } + + #[test] + fn truncated_tdvp_follows_the_quench() { + // Sixteen sites at bond 8: the state no longer fits, so this is + // truncation and projection error, held to the free-fermion answer. + let n = 16; + let chain = Chain::ising(n, 1.0, 1.0); + let plus = Mps::product(&vec![[C::ONE, C::ONE]; n]); + let cfg = TdvpConfig { dt: 0.1, max_bond: 8, ..TdvpConfig::default() }; + let mut ev = Evolution::new(&chain, &plus, cfg); + for _ in 0..10 { + ev.step(); + } + let ff = ising_quench(n, 1.0, 1.0, ev.time()); + let x = ev.expect(Op::x()); + let zz = ev.bond_expect(Op::z(), Op::z()); + for (i, (got, want)) in x.iter().zip(&ff.x).enumerate() { + assert!((got - want).abs() < 1e-6, "x[{i}]: {got} vs {want}"); + } + for (i, (got, want)) in zz.iter().zip(&ff.zz).enumerate() { + assert!((got - want).abs() < 1e-6, "zz[{i}]: {got} vs {want}"); + } + let echo = ev.overlap(&plus).norm2(); + assert!((echo - ff.echo).abs() < 1e-6, "echo: {echo} vs {}", ff.echo); + // Entanglement peaks mid-chain. + let s = ev.entropies(); + assert!(s[n / 2 - 1] > s[0] && s[n / 2 - 1] > s[n - 2]); + } + + #[test] + fn energy_is_conserved_and_runs_repeat() { + let n = 12; + let chain = Chain::heisenberg(n, 1.0, 0.8, 0.1); + let start = Mps::basis(&(0..n).map(|i| i % 2 == 0).collect::>()); + let cfg = TdvpConfig { dt: 0.1, max_bond: 12, ..TdvpConfig::default() }; + let mut ev = Evolution::new(&chain, &start, cfg); + let e0 = ev.energy(); + let want = 0.8 * 0.25 * -(n as f64 - 1.0); + assert!((e0 - want).abs() < 1e-12, "{e0} vs {want}"); + for _ in 0..10 { + ev.step(); + } + assert!((ev.energy() - e0).abs() < 1e-9, "{} vs {e0}", ev.energy()); + assert!((ev.norm() - 1.0).abs() < 1e-12); + // Same inputs, same bits, at any thread count. + let mut one = Evolution::new(&chain, &start, TdvpConfig { threads: 1, ..cfg }); + let mut many = Evolution::new(&chain, &start, TdvpConfig { threads: 5, ..cfg }); + for _ in 0..10 { + one.step(); + many.step(); + } + assert_eq!(one.state(), many.state()); + assert_eq!(one.state(), ev.state()); + } + + #[test] + fn threaded_complex_effective_hamiltonians_are_bit_identical() { + let mpo = Chain::heisenberg(8, 1.0, 0.7, 0.2).mpo(); + let eng = Engine::new(&mpo); + let (dl, dr, w) = (56usize, 48usize, mpo.w); + let mut next = rng(3); + let mut c = || C { re: next(), im: next() }; + let l: Vec = (0..dl * w * dl).map(|_| c()).collect(); + let r: Vec = (0..dr * w * dr).map(|_| c()).collect(); + let theta: Vec = (0..dl * 4 * dr).map(|_| c()).collect(); + let site: Vec = (0..dl * 2 * dr).map(|_| c()).collect(); + let r1: Vec = r.clone(); + let mut one = vec![C::ZERO; theta.len()]; + let mut one1 = vec![C::ZERO; site.len()]; + eng.apply(&l, &r, &theta, &mut one, dl, dr, 1); + eng.apply1(&l, &r1, &site, &mut one1, dl, dr, 1); + for t in [2, 3, 7] { + let mut many = vec![C::ZERO; theta.len()]; + let mut many1 = vec![C::ZERO; site.len()]; + eng.apply(&l, &r, &theta, &mut many, dl, dr, t); + eng.apply1(&l, &r1, &site, &mut many1, dl, dr, t); + assert_eq!(one, many, "threads={t}"); + assert_eq!(one1, many1, "threads={t}"); + } + } +} diff --git a/wai/wai-quantum/src/repro.rs b/wai/wai-quantum/src/repro.rs index 7233a6977..ee1f70afa 100644 --- a/wai/wai-quantum/src/repro.rs +++ b/wai/wai-quantum/src/repro.rs @@ -15,7 +15,8 @@ feature = "quantum_phasor", feature = "quantum_kernel", feature = "quantum_tn", - feature = "quantum_sv" + feature = "quantum_sv", + feature = "quantum_tdvp" ))] mod trig; #[cfg(any( @@ -26,7 +27,8 @@ mod trig; feature = "quantum_phasor", feature = "quantum_kernel", feature = "quantum_tn", - feature = "quantum_sv" + feature = "quantum_sv", + feature = "quantum_tdvp" ))] pub use trig::sin_cos; diff --git a/wai/wai-quantum/tests/bit_identity.rs b/wai/wai-quantum/tests/bit_identity.rs index 4b9aab75b..752ad9ce4 100644 --- a/wai/wai-quantum/tests/bit_identity.rs +++ b/wai/wai-quantum/tests/bit_identity.rs @@ -25,7 +25,8 @@ feature = "quantum_qio", feature = "quantum_uf", feature = "quantum_noisy", - feature = "quantum_dmrg" + feature = "quantum_dmrg", + feature = "quantum_tdvp" ))] use wai_quantum::*; @@ -474,3 +475,36 @@ fn dmrg_ground_states() { v.extend(vecs.iter().map(|x| x.abs())); assert_eq!(digest(&v), 0x3b28_8bdd_be3f_48a2); } + +#[test] +fn tdvp_time_evolution() { + use quantum_dmrg::{Chain, Mps, Op, Scalar, C}; + use quantum_tdvp::*; + let mut v = Vec::new(); + // A truncated Heisenberg quench from the Néel state, and the Ising quench + // with its free-fermion referee. + let n = 10; + let neel = Mps::basis(&(0..n).map(|i| i % 2 == 1).collect::>()); + let mut ev = Evolution::new(&Chain::heisenberg(n, 1.0, 0.7, 0.1), &neel, TdvpConfig { dt: 0.1, max_bond: 6, threads: 1, ..TdvpConfig::default() }); + for _ in 0..6 { + v.push(ev.step().discarded); + } + v.extend(ev.expect(Op::z())); + v.extend(ev.bond_expect(Op::sz(), Op::sz())); + v.extend(ev.entropies()); + v.push(ev.energy()); + let o = ev.overlap(&neel); + v.extend([o.re, o.im]); + let plus = Mps::product(&vec![[C::ONE, C::ONE]; 8]); + let mut ev = Evolution::new(&Chain::ising(8, 1.0, 0.9), &plus, TdvpConfig { dt: 0.2, max_bond: 16, cutoff: 0.0, threads: 1, ..TdvpConfig::default() }); + for _ in 0..3 { + ev.step(); + } + v.extend(ev.expect(Op::x())); + let ff = ising_quench(8, 1.0, 0.9, ev.time()); + v.extend(ff.x); + v.extend(ff.zz); + v.push(ff.echo); + assert_eq!(digest(&v), 0xe1b4_fd28_af23_50cd); +} + diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index a3d176385..c89a431b9 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.34" +version = "0.3.35" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-rs/Cargo.toml b/wai/wai-rs/Cargo.toml index 593158e50..45a21bd67 100644 --- a/wai/wai-rs/Cargo.toml +++ b/wai/wai-rs/Cargo.toml @@ -368,7 +368,7 @@ quantum_arch = ["dep:wai-quantum", "wai-quantum/quantum_arch"] # qLDPC memory at circuit level: bivariate-bicycle codes, BP-OSD. quantum_ldpc = ["quantum_frame", "wai-quantum/quantum_ldpc"] # Quantum-inspired optimisation: simulated bifurcation, annealing, exact referee. -quantum_qio = ["dep:wai-quantum", "wai-quantum/quantum_qio", "wai-quantum/quantum_dmrg"] +quantum_qio = ["dep:wai-quantum", "wai-quantum/quantum_qio", "wai-quantum/quantum_dmrg", "wai-quantum/quantum_tdvp"] 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 diff --git a/wai/wai-rs/src/bin/wai_quantum.rs b/wai/wai-rs/src/bin/wai_quantum.rs index 905a4a2d2..a4f2c7105 100644 --- a/wai/wai-rs/src/bin/wai_quantum.rs +++ b/wai/wai-rs/src/bin/wai_quantum.rs @@ -167,7 +167,9 @@ fn usage() -> ! { wai-quantum gbs … --check SAMPLES.txt (observed patterns vs exact marginals)\n\n\ SPIN CHAINS (DMRG ground states, exact referees):\n \ wai-quantum dmrg [--chain ising|heisenberg] [-n N] [--j J] [--jz JZ] [--h H] [--bond D] [--sweeps S]\n \ - [--cutoff C] [--exact] [--json]\n\n\ + [--cutoff C] [--exact] [--json]\n \ + wai-quantum tdvp [--chain ising|heisenberg] [-n N] [--j J] [--jz JZ] [--h H] [--start plus|neel|up]\n \ + [--bond D] [--dt DT] [--t T] [--every K] [--cutoff C] [--json] (time evolution by TDVP)\n\n\ OPTIMISATION (simulated bifurcation, annealing, exact):\n \ wai-quantum qio (--gset FILE | --sk N [--instance S]) [--method dsb|bsb|hdsb|hbsb|sa|exact]\n \ [--steps N] [--trials T] [--dt D] [--c1 C] [--gamma G] [--beta0 B --beta1 B] [--descend]\n \ @@ -1030,6 +1032,90 @@ fn main() { } } + // Time evolution of a spin chain by TDVP from a product state. The + // Ising chain started with every σˣ = +1 is held at each report to + // its free-fermion solution. + "tdvp" => { + use wai::quantum_dmrg::{Chain, Mps, Op, Scalar, C}; + use wai::quantum_tdvp::*; + let n: usize = num(args, "--sites", "-n", 40); + if n < 2 { + die("-n: at least 2 sites"); + } + let kind = opt(args, "--chain", "").unwrap_or_else(|| "ising".into()); + let j: f64 = num(args, "--j", "", 1.0); + let h: f64 = num(args, "--h", "", if kind == "heisenberg" { 0.0 } else { 1.0 }); + let chain = match kind.as_str() { + "ising" => Chain::ising(n, j, h), + "heisenberg" => Chain::heisenberg(n, j, num(args, "--jz", "", j), h), + _ => die("--chain: ising or heisenberg"), + }; + let start_kind = opt(args, "--start", "").unwrap_or_else(|| if kind == "ising" { "plus".into() } else { "neel".into() }); + let start = match start_kind.as_str() { + "plus" => Mps::product(&vec![[C::ONE, C::ONE]; n]), + "neel" => Mps::basis(&(0..n).map(|i| i % 2 == 1).collect::>()), + "up" => Mps::basis(&vec![false; n]), + _ => die("--start: plus, neel or up"), + }; + let dt: f64 = num(args, "--dt", "", 0.05); + let t_end: f64 = num(args, "--t", "", 2.0); + if dt.is_nan() || dt <= 0.0 || t_end.is_nan() || t_end < 0.0 { + die("--dt must be positive and --t not negative"); + } + let cfg = TdvpConfig { dt, max_bond: num(args, "--bond", "", 32), cutoff: num(args, "--cutoff", "", 1e-20), ..TdvpConfig::default() }; + let steps = (t_end / dt).round() as usize; + let every: usize = num(args, "--every", "", ((0.5 / dt).round() as usize).max(1)); + // The observable at the middle site: σˣ for the Ising chain, σᶻ + // for the Heisenberg chain. + let (obs, obs_name) = if kind == "ising" { (Op::x(), "x") } else { (Op::z(), "z") }; + let referee = kind == "ising" && start_kind == "plus"; + let mut ev = Evolution::new(&chain, &start, cfg); + let e0 = ev.energy(); + let mut discarded = 0.0; + let mut rows: Vec = Vec::new(); + let json = flag(args, "--json"); + if !json { + println!("{kind} chain, {n} sites, start {start_kind}, bond ≤ {}, dt {dt}", cfg.max_bond); + print!(" t {obs_name}(mid) "); + if referee { + print!("exact max error rate exact rate "); + } + println!("S(mid) bond |ΔE| discarded"); + } + for s in 1..=steps { + discarded += ev.step().discarded; + if s % every != 0 && s != steps { + continue; + } + let t = ev.time(); + let o = ev.expect(obs); + let s_mid = ev.entropies()[n / 2 - 1]; + let de = (ev.energy() - e0).abs(); + let rate = -ev.overlap(&start).norm2().ln() / n as f64; + let ff = referee.then(|| ising_quench(n, j, h, t)); + let max_err = ff.as_ref().map(|f| o.iter().zip(&f.x).map(|(p, q)| (p - q).abs()).fold(0.0, f64::max)); + if json { + rows.push(format!( + "{{\"bond\":{},\"discarded\":{discarded},\"energy_drift\":{de},\"entropy\":{s_mid},\"exact\":{},\"exact_rate\":{},\"max_error\":{},\"mid\":{},\"rate\":{rate},\"t\":{t}}}", + ev.max_bond(), + ff.as_ref().map_or("null".into(), |f| f.x[n / 2].to_string()), + ff.as_ref().map_or("null".into(), |f| (-f.echo.ln() / n as f64).to_string()), + max_err.map_or("null".into(), |e| e.to_string()), + o[n / 2] + )); + } else { + print!("{t:6.2} {:+.10} ", o[n / 2]); + if let (Some(f), Some(e)) = (&ff, max_err) { + print!("{:+.10} {e:.2e} {rate:.8} {:.8} ", f.x[n / 2], -f.echo.ln() / n as f64); + } + println!("{s_mid:.4} {:4} {de:.1e} {discarded:.1e}", ev.max_bond()); + } + } + if json { + println!("{{\"chain\":\"{kind}\",\"dt\":{dt},\"rows\":[{}],\"sites\":{n},\"start\":\"{start_kind}\"}}", rows.join(",")); + } + } + // Quantum-inspired optimisation: Max-Cut (G-set text) or a spin // glass, by simulated bifurcation, annealing, or exact enumeration. // Max-Cut reports cuts; a spin glass reports energies. diff --git a/wai/wai-rs/src/lib.rs b/wai/wai-rs/src/lib.rs index 03fc06905..399e9fe68 100644 --- a/wai/wai-rs/src/lib.rs +++ b/wai/wai-rs/src/lib.rs @@ -663,6 +663,9 @@ pub use wai_quantum::quantum_qio; #[cfg(feature = "quantum_qio")] pub use wai_quantum::quantum_dmrg; +#[cfg(feature = "quantum_qio")] +pub use wai_quantum::quantum_tdvp; + #[cfg(any(feature = "quantum_tn", feature = "quantum_sv"))] pub use wai_quantum::gates; #[cfg(feature = "quantum_qir")] diff --git a/wai/wai-web/Cargo.lock b/wai/wai-web/Cargo.lock index 4b7475b6b..fbbb16335 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.34" +version = "0.3.35" dependencies = [ "blake3", "ed25519-dalek",