From 5cd36795cde77706b0ae91e8867747b88389a93e Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Tue, 6 Oct 2026 17:15:52 -0400 Subject: [PATCH 1/3] wai-quantum 0.3.37: complex neural quantum states, and referees that cannot stop on an excited state Fix. In 0.3.36, Lanczos stopped a restart as soon as its lowest Ritz pair converged. That is safe from DMRG's warm starts, but from the referees' cold ones that pair can be an excited state. Exact diagonalisation of the 8-site Majumdar-Ghosh ring returned -2.595 instead of -3, and a variational energy came out below it. Early stopping is now opt-in, used by the DMRG sweeps only; exact diagonalisation, full CI and the spin-model referee build every Krylov vector. The DMRG and NQS pins were re-taken. New: the RBM, VMC and stochastic reconfiguration are generic over real and complex parameters, and the real path keeps its bits. - The complex path adds exp, ln cosh and tanh, the Hermitian reconfiguration step in either space, and a Hermitian Cholesky. - New lattices: J1-J2 models, plus next-nearest bonds on chains and diagonal bonds on square lattices. - A complex network learns sign structures nobody supplied. On 8-site rings it gets within 2e-3 where a real one stays >10% off, and the CLI solves the 8-site Majumdar-Ghosh ring to rounding. On 12-site rings: Majumdar-Ghosh 2e-6, J1-J2 (0.3) 3.5e-4. - Measured and documented: - 10-site rings settle on a momentum-zero excited eigenstate, because the ground state has momentum pi. - At J2 = 0.5 on the 4x4 torus, alpha = 1 stays about 3% above exact. The NQS pin now covers a complex machine. All 22 pins pass natively and under wasm32-wasip2. The nqs CLI gains --complex, --model j1j2 and --periodic, byte-identical under WASI. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum-serve/Cargo.lock | 2 +- wai/wai-quantum/Cargo.lock | 2 +- wai/wai-quantum/Cargo.toml | 2 +- wai/wai-quantum/README.md | 7 +- wai/wai-quantum/examples/nqs_ground.rs | 103 +++-- wai/wai-quantum/examples/nqs_probe.rs | 20 + wai/wai-quantum/src/quantum_chem.rs | 2 +- wai/wai-quantum/src/quantum_dmrg.rs | 22 +- wai/wai-quantum/src/quantum_nqs.rs | 615 ++++++++++++++++++------- wai/wai-quantum/src/repro.rs | 20 +- wai/wai-quantum/tests/bit_identity.rs | 13 +- wai/wai-rs/Cargo.lock | 2 +- wai/wai-rs/src/bin/wai_quantum.rs | 61 ++- wai/wai-web/Cargo.lock | 2 +- 14 files changed, 621 insertions(+), 252 deletions(-) create mode 100644 wai/wai-quantum/examples/nqs_probe.rs diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index ea3a590c1..12d6a3727 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.36" +version = "0.3.37" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index 835fc4e62..b24857177 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.36" +version = "0.3.37" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index dac363755..a44d2ae1c 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "wai-quantum" -version = "0.3.36" +version = "0.3.37" edition = "2024" rust-version = "1.97" license = "Apache-2.0" diff --git a/wai/wai-quantum/README.md b/wai/wai-quantum/README.md index 9544d8bf0..d92896854 100644 --- a/wai/wai-quantum/README.md +++ b/wai/wai-quantum/README.md @@ -69,7 +69,7 @@ never creates information. | `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_chem` | **molecules by DMRG**: FCIDUMP integrals, the second-quantised Hamiltonian on spin orbitals compiled into an exact MPO (minimum vertex covers, no numerical compression), DMRG from Hartree–Fock with complete bases, and full configuration interaction as the referee. LiH and H₂O reach full CI to `10⁻¹³` | -| `quantum_nqs` | **neural quantum states**: spin models on any graph, a restricted Boltzmann machine (with Marshall's sign on bipartite lattices), variational Monte Carlo by Metropolis chains or exact enumeration, and stochastic reconfiguration solved in whichever of parameter or sample space is smaller. Held to exact diagonalisation and to DMRG on the same operator sum | +| `quantum_nqs` | **neural quantum states**: spin models on any graph (Heisenberg, frustrated J₁–J₂, transverse-field Ising), a restricted Boltzmann machine with real parameters (and Marshall's sign on bipartite lattices) or complex ones that learn a sign structure nobody supplied, variational Monte Carlo by Metropolis chains or exact enumeration, and stochastic reconfiguration solved in whichever of parameter or sample space is smaller. Held to exact diagonalisation and to DMRG on the same operator sum | | `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 | @@ -242,6 +242,7 @@ weight that accounts for the difference. ```sh wai-quantum nqs --lx 8 --exact --steps 200 --lr 0.1 # every configuration enumerated: no sampling noise wai-quantum nqs --model ising --lx 4 --ly 4 --h 3.044 --alpha 1 --steps 400 +wai-quantum nqs --model j1j2 --lx 8 --periodic --complex --exact --steps 250 --lr 0.1 # frustrated: the sign is learned ``` Neural-network wavefunctions are offered as evidence about systems that no @@ -249,7 +250,9 @@ exact method reaches. `quantum_nqs` trains one the standard way and holds it to exact answers. On a 4 × 4 lattice at the Ising model's critical field, it lands within its own error bar of exact diagonalisation. Past exact reach, the same Hamiltonian goes to DMRG through its exact MPO, so two independent -methods meet on one number. +methods meet on one number. For frustrated models, whose signs no rule gives, +a complex network learns them. On the Majumdar–Ghosh ring it reaches the +exact energy to rounding. ## Noisy circuits diff --git a/wai/wai-quantum/examples/nqs_ground.rs b/wai/wai-quantum/examples/nqs_ground.rs index 79cf6ca9e..4bdbbcd85 100644 --- a/wai/wai-quantum/examples/nqs_ground.rs +++ b/wai/wai-quantum/examples/nqs_ground.rs @@ -1,57 +1,76 @@ //! A neural quantum state for a spin model's ground state, by variational //! Monte Carlo with stochastic reconfiguration, against exact -//! diagonalisation (up to 22 spins) or DMRG on the same operator sum. +//! diagonalisation (up to 20 spins) or DMRG on the same operator sum. //! //! ```text -//! cargo run --release --example nqs_ground --features quantum_nqs -- heisenberg|ising LX LY ALPHA ITERS SAMPLES [h] [dmrg-bond] [threads] +//! cargo run --release --example nqs_ground --features quantum_nqs -- \ +//! --model heisenberg|j1j2|ising --lx L --ly L [--j2 J] [--h H] [--periodic] +//! [--complex] [--marshall] [--alpha A] [--steps N] [--samples S] [--lr R] +//! [--dmrg-bond D] [--threads T] //! ``` -use wai_quantum::quantum_dmrg::{dmrg_mpo, DmrgConfig, Mps}; +use wai_quantum::quantum_dmrg::{dmrg_mpo, DmrgConfig, Mps, C}; use wai_quantum::quantum_nqs::*; -fn main() { - let a: Vec = std::env::args().collect(); - let (lx, ly): (usize, usize) = (a[2].parse().unwrap(), a[3].parse().unwrap()); - let alpha: usize = a[4].parse().unwrap(); - let iters: usize = a[5].parse().unwrap(); - let samples: usize = a[6].parse().unwrap(); - let h: f64 = a.get(7).map_or(3.044, |x| x.parse().unwrap()); - let bond: usize = a.get(8).map_or(128, |x| x.parse().unwrap()); - let threads: usize = a.get(9).map_or(VmcConfig::default().threads, |x| x.parse().unwrap()); - let n = lx * ly; - let edges = if ly == 1 { chain_edges(n, false) } else { square_edges(lx, ly, false) }; - let heis = a[1] == "heisenberg"; - let model = if heis { SpinModel::heisenberg(n, &edges, 1.0) } else { SpinModel::ising(n, &edges, 1.0, h) }; - let sector = if heis { Some(0) } else { None }; - let t0 = std::time::Instant::now(); - let (reference, how) = match model.exact_ground_energy(sector) { - Some(e) if n <= 20 => (e, "exact diagonalisation".to_string()), - _ => { - let mpo = model.opsum().mpo(); - let r = dmrg_mpo(&mpo, Mps::random(n, 8, 1), &DmrgConfig { max_bond: bond, sweeps: 6, threads, ..DmrgConfig::default() }); - (r.energy, format!("DMRG, bond {bond}, discarded {:.1e}", r.discarded)) - } - }; - println!("{} {lx}×{ly}: reference {reference:.8} ({how}, {:.1}s)", a[1], t0.elapsed().as_secs_f64()); - let mut rbm = Rbm::new(n, alpha, 0.01, 7); - if heis { - rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); - } - println!("RBM α={alpha}: {} parameters, {samples} samples a step", rbm.params()); - let cfg = VmcConfig { samples, chains: 16, magnetisation: sector, lr: 0.02, shift: 1e-3, threads, ..VmcConfig::default() }; +fn arg(a: &[String], key: &str, default: T) -> T { + a.iter().position(|x| x == key).and_then(|i| a.get(i + 1)).and_then(|v| v.parse().ok()).unwrap_or(default) +} + +fn train(model: SpinModel, rbm: Rbm, cfg: VmcConfig, steps: usize, reference: f64) { + println!("RBM α={}: {} parameters, {} samples a step", rbm.m / rbm.n, rbm.params(), cfg.samples); let mut vmc = Vmc::new(model, rbm, cfg); let t0 = std::time::Instant::now(); - let mut tail = Vec::new(); - for it in 1..=iters { + for it in 1..=steps { let e = vmc.step(); - if it > iters * 4 / 5 { - tail.push(e.energy); - } - if it % (iters / 10).max(1) == 0 { + if it % (steps / 10).max(1) == 0 { println!("step {it:>5} E {:.6} ± {:.6} variance {:.2e} acceptance {:.2} relative error {:.2e} {:.0}s", e.energy, e.error, e.variance, e.acceptance, (e.energy - reference) / reference.abs(), t0.elapsed().as_secs_f64()); } } - // A long final estimate at the trained parameters. - vmc.cfg.samples = samples * 8; + vmc.cfg.samples *= 8; let f = vmc.estimate(); - println!("final E {:.6} ± {:.6} variance {:.2e} relative error {:.2e} (mean of last fifth {:.6})", f.energy, f.error, f.variance, (f.energy - reference) / reference.abs(), tail.iter().sum::() / tail.len() as f64); + println!("final E {:.6} ± {:.6} variance {:.2e} relative error {:.2e}", f.energy, f.error, f.variance, (f.energy - reference) / reference.abs()); +} + +fn main() { + let a: Vec = std::env::args().collect(); + let kind: String = arg(&a, "--model", "heisenberg".to_string()); + let (lx, ly): (usize, usize) = (arg(&a, "--lx", 4), arg(&a, "--ly", 4)); + let periodic = a.iter().any(|x| x == "--periodic"); + let threads: usize = arg(&a, "--threads", VmcConfig::default().threads); + let n = lx * ly; + let nn = if ly == 1 { chain_edges(n, periodic) } else { square_edges(lx, ly, periodic) }; + let nnn = if ly == 1 { chain_next_nearest(n, periodic) } else { square_diagonals(lx, ly, periodic) }; + let model = match kind.as_str() { + "heisenberg" => SpinModel::heisenberg(n, &nn, 1.0), + "j1j2" => SpinModel::j1j2(n, &nn, &nnn, 1.0, arg(&a, "--j2", 0.5)), + "ising" => SpinModel::ising(n, &nn, 1.0, arg(&a, "--h", 3.044)), + other => panic!("{other}"), + }; + let sector = if model.conserves_magnetisation() { Some(0) } else { None }; + let t0 = std::time::Instant::now(); + let (reference, how) = match if n <= 20 { model.exact_ground_energy(sector) } else { None } { + Some(e) => (e, "exact diagonalisation".to_string()), + None => { + let bond = arg(&a, "--dmrg-bond", 128usize); + let r = dmrg_mpo(&model.opsum().mpo(), Mps::random(n, 8, 1), &DmrgConfig { max_bond: bond, sweeps: 6, threads, ..DmrgConfig::default() }); + (r.energy, format!("DMRG, bond {bond}, discarded {:.1e}", r.discarded)) + } + }; + println!("{kind} {lx}×{ly}{}: reference {reference:.8} ({how}, {:.1}s)", if periodic { " periodic" } else { "" }, t0.elapsed().as_secs_f64()); + let alpha: usize = arg(&a, "--alpha", 2); + let cfg = VmcConfig { samples: arg(&a, "--samples", 1000), chains: 16, magnetisation: sector, lr: arg(&a, "--lr", 0.02), shift: 1e-3, threads, ..VmcConfig::default() }; + let steps: usize = arg(&a, "--steps", 600); + let marshall = a.iter().any(|x| x == "--marshall"); + if a.iter().any(|x| x == "--complex") { + let mut rbm: Rbm = Rbm::complex(n, alpha, 0.01, 7); + if marshall { + rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); + } + train(model, rbm, cfg, steps, reference); + } else { + let mut rbm = Rbm::new(n, alpha, 0.01, 7); + if marshall { + rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); + } + train(model, rbm, cfg, steps, reference); + } } diff --git a/wai/wai-quantum/examples/nqs_probe.rs b/wai/wai-quantum/examples/nqs_probe.rs new file mode 100644 index 000000000..6a9766664 --- /dev/null +++ b/wai/wai-quantum/examples/nqs_probe.rs @@ -0,0 +1,20 @@ +use wai_quantum::quantum_nqs::*; +fn main() { + for n in [10usize, 12] { + let mg = SpinModel::j1j2(n, &chain_edges(n, true), &chain_next_nearest(n, true), 1.0, 0.5); + let heis = SpinModel::heisenberg(n, &chain_edges(n, true), 1.0); + let j1j2 = SpinModel::j1j2(n, &chain_edges(n, true), &chain_next_nearest(n, true), 1.0, 0.3); + for (name, model) in [("majumdar-ghosh", mg), ("heisenberg (no sign)", heis), ("j1j2 0.3", j1j2)] { + let exact = model.exact_ground_energy(Some(0)).unwrap(); + for scale in [0.05, 0.3] { + let cfg = VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, shift: 1e-3, ..VmcConfig::default() }; + let mut cplx = Vmc::new(model.clone(), Rbm::complex(n, 2, scale, 3), cfg); + let mut ec = cplx.step(); + for _ in 0..400 { + ec = cplx.step(); + } + println!("n={n} {name:22} scale {scale}: exact {exact:.6} complex {:.6} rel {:.2e} variance {:.1e}", ec.energy, (ec.energy - exact) / exact.abs(), ec.variance); + } + } + } +} diff --git a/wai/wai-quantum/src/quantum_chem.rs b/wai/wai-quantum/src/quantum_chem.rs index 9f0805d60..7e5ca1f19 100644 --- a/wai/wai-quantum/src/quantum_chem.rs +++ b/wai/wai-quantum/src/quantum_chem.rs @@ -509,7 +509,7 @@ pub fn fci(fd: &Fcidump) -> Option { // a little of every other determinant, so no symmetry hides the ground // state. let start: Vec = (0..dim as u64).map(|x| if x == 0 { 1.0 } else { 1e-3 * (1.0 + (x.wrapping_mul(2_654_435_761) % 1000) as f64 / 1000.0) }).collect(); - let (energy, _) = lanczos(&apply, &start, 60, 40, 1e-9); + let (energy, _) = lanczos(&apply, &start, 60, 40, 1e-9, false); Some(FciResult { energy, dim }) } diff --git a/wai/wai-quantum/src/quantum_dmrg.rs b/wai/wai-quantum/src/quantum_dmrg.rs index 74b0bc162..e117c823f 100644 --- a/wai/wai-quantum/src/quantum_dmrg.rs +++ b/wai/wai-quantum/src/quantum_dmrg.rs @@ -18,9 +18,11 @@ //! placed exactly once. [`dmrg_mpo`] runs DMRG on any MPO from any start. //! - **Two-site DMRG.** [`dmrg`] sweeps a matrix-product state left and right. //! - At each bond it finds the lowest eigenvector of the effective -//! Hamiltonian by Lanczos, with full re-orthogonalisation. Lanczos stops -//! as soon as the lowest Ritz pair's residual estimate meets the -//! tolerance. +//! Hamiltonian by Lanczos, with full re-orthogonalisation. Started from +//! the current state, Lanczos stops as soon as the lowest Ritz pair's +//! residual estimate meets the tolerance. The exact referees, started +//! cold, build every Krylov vector, since there the lowest Ritz pair can +//! converge on an excited state first. //! - It splits that vector by the eigenvectors of its reduced density //! matrix. It keeps at most `max_bond` of them and drops those below a //! weight cutoff. @@ -922,7 +924,7 @@ fn min_vertex_cover(nu: usize, nv: usize, adj: &[Vec]) -> (Vec, Vec /// Lanczos for the lowest eigenpair of a symmetric operator given by /// `apply`, from `start`, with full re-orthogonalisation. -pub(crate) fn lanczos(apply: &dyn Fn(&[f64], &mut [f64]), start: &[f64], krylov: usize, restarts: usize, tol: f64) -> (f64, Vec) { +pub(crate) fn lanczos(apply: &dyn Fn(&[f64], &mut [f64]), start: &[f64], krylov: usize, restarts: usize, tol: f64, early: bool) -> (f64, Vec) { let dim = start.len(); let norm = |v: &[f64]| v.iter().map(|x| x * x).sum::().sqrt(); let mut v = start.to_vec(); @@ -955,8 +957,10 @@ pub(crate) fn lanczos(apply: &dyn Fn(&[f64], &mut [f64]), start: &[f64], krylov: } // The lowest Ritz pair's residual is β times the last component // of its vector; once that meets the tolerance, more vectors - // cannot help. - if j >= 1 { + // cannot help. Only for a warm start (DMRG's current state): from + // a cold one the lowest Ritz pair can be an excited state that + // converges first, and a referee must not stop on it. + if early && j >= 1 { let m = alpha.len(); let mut t = vec![0.0; m * m]; for i in 0..m { @@ -1058,7 +1062,7 @@ pub fn exact_ground_energy(chain: &Chain) -> Option { 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) + Some(lanczos(&apply, &start, 60, 20, 1e-10, false).0) } /// The exact ground energy of the open transverse-field Ising chain @@ -1790,7 +1794,7 @@ pub fn dmrg_mpo(mpo: &Mpo, start: Mps, cfg: &DmrgConfig) -> DmrgResult { let theta = two_site(&mps.sites[i], &mps.sites[i + 1], dl, dm, dr); let (l, r) = (&lefts[i], &rights[i + 2]); let apply = |x: &[f64], y: &mut [f64]| eng.apply(i, l, r, x, y, dl, dr, cfg.threads); - let (e, v) = lanczos(&apply, &theta, cfg.krylov, cfg.restarts, cfg.tol); + let (e, v) = lanczos(&apply, &theta, cfg.krylov, cfg.restarts, cfg.tol, true); energy = e; let (u, rest, disc, weights) = split(&v, dl * 2, 2 * dr, cfg.max_bond, cfg.cutoff); let k = weights.len(); @@ -1807,7 +1811,7 @@ pub fn dmrg_mpo(mpo: &Mpo, start: Mps, cfg: &DmrgConfig) -> DmrgResult { let theta = two_site(&mps.sites[i], &mps.sites[i + 1], dl, dm, dr); let (l, r) = (&lefts[i], &rights[i + 2]); let apply = |x: &[f64], y: &mut [f64]| eng.apply(i, l, r, x, y, dl, dr, cfg.threads); - let (e, v) = lanczos(&apply, &theta, cfg.krylov, cfg.restarts, cfg.tol); + let (e, v) = lanczos(&apply, &theta, cfg.krylov, cfg.restarts, cfg.tol, true); energy = e; // Split on the right: rows of the right factor are eigenvectors // of θᵀθ. diff --git a/wai/wai-quantum/src/quantum_nqs.rs b/wai/wai-quantum/src/quantum_nqs.rs index 8221e5250..ddc0b73e4 100644 --- a/wai/wai-quantum/src/quantum_nqs.rs +++ b/wai/wai-quantum/src/quantum_nqs.rs @@ -7,11 +7,16 @@ //! in its checkable form, deterministic and held to exact answers. //! //! - **Models.** [`SpinModel`] is a spin-½ Hamiltonian on any graph, with -//! exchange, Ising, and field terms. Presets are Heisenberg and -//! transverse-field Ising, on chains and square lattices. +//! exchange, Ising, and field terms. Presets are Heisenberg, the frustrated +//! J₁–J₂ model, and transverse-field Ising, on chains and square lattices. //! - **The network.** [`Rbm`] is a restricted Boltzmann machine with `α·n` -//! hidden units. On a bipartite lattice it can carry Marshall's sign, which -//! makes the Heisenberg ground state positive in what remains. +//! hidden units, its parameters real or complex ([`Amplitude`]). +//! - A real network is positive. On a bipartite lattice it can carry +//! Marshall's sign, which makes the Heisenberg ground state positive in +//! what remains. +//! - A complex network carries a phase, so training can find a sign +//! structure nobody supplied. Frustrated models need that, because no +//! sign rule is known for them. //! - **Variational Monte Carlo.** [`Vmc`] samples `|ψ|²` by Metropolis //! chains: single flips, or swaps of antiparallel spins at fixed //! magnetisation. Or it enumerates every configuration (at most 24 spins) @@ -37,6 +42,21 @@ //! every configuration enumerated reach the exact energy to `10⁻⁷`. //! - **Sampling.** Sampled energies agree with enumerated ones within their //! error bars, and steps repeat bit for bit at any thread count. +//! - **Complex networks.** A complex network with real parameters is the real +//! one, to `10⁻¹³`. Its ratios equal ratios of its amplitudes. Its +//! holomorphic log-derivatives agree with differences along a real step. +//! Both spaces take the same Hermitian reconfiguration step. +//! - **A sign nobody supplied.** Exact diagonalisation gives the +//! Majumdar–Ghosh ring (`J₂ = J₁/2`) its exact `−3n/8`. On 8-site rings +//! with no sign rule given, Heisenberg and Majumdar–Ghosh, a real network +//! stays more than 10 % off. A complex one comes within `2·10⁻³`, with +//! every configuration enumerated. On 12-site rings, 400 steps reach: +//! +//! | ring (no sign given) | relative error | +//! |-----------------------|----------------| +//! | Majumdar–Ghosh | 2·10⁻⁶ | +//! | J₁–J₂, `J₂ = 0.3` | 3.5·10⁻⁴ | +//! | Heisenberg | 2·10⁻³ | //! - **Against exact answers** (`nqs_ground` example; relative error of the //! final estimate, with its statistical error): //! @@ -51,11 +71,32 @@ //! an upper bound (discarded weight `2·10⁻⁵`). The network was still //! descending when it stopped there. //! +//! # Measured +//! +//! On 10-site rings every model above settled about 9 % high, on a state +//! with essentially zero variance: an eigenstate, but not the ground state. +//! When half the ring's length is odd, the ground state has momentum `π`. A +//! network started near the uniform state begins at momentum zero, and +//! training stays there and finds the lowest state it can reach. The +//! variance cannot catch this, since an excited eigenstate has none. Only a +//! referee can. +//! +//! The same trap caught the referee once. Lanczos stopped a restart once its +//! lowest Ritz pair had converged, which is safe from DMRG's warm starts. +//! From a cold start that pair can be an excited state, and the exact +//! diagonalisation of the 8-site Majumdar–Ghosh ring returned `−2.595` +//! instead of `−3`. The network's energy, a true upper bound, came out below +//! it. The referees now build every Krylov vector. +//! //! # Honest boundaries //! -//! - **A real network.** The sign structure must be known, as Marshall's is -//! on bipartite lattices. Frustrated models, whose signs are not known, -//! need complex parameters or a learned sign, which are not here. +//! - **The hardest point is not reached.** At `J₂ = 0.5` on the periodic +//! 4 × 4 lattice, exact diagonalisation gives `−8.45792335`. At `α = 1`, +//! 600 steps of 1000 samples leave all three networks about 3 % above it: +//! complex with Marshall's sign, complex without, and real with it. At +//! this size and budget the network is the limit, not the sampling. +//! - **Symmetry is not imposed.** A network is not projected onto a momentum +//! or spin sector, so it can settle in the wrong one, as above. //! - **Error bars ignore autocorrelation.** The reported error is the //! standard error of uncorrelated samples. Chains thinned by two sweeps are //! correlated, so the true error is larger. @@ -63,8 +104,8 @@ //! factored in full each step. That bounds the network and sample sizes //! long before iterative solvers would. -use crate::quantum_dmrg::{Op, OpSum, default_threads, lanczos}; -use crate::repro::{exp, ln}; +use crate::quantum_dmrg::{C, Op, OpSum, Scalar, default_threads, lanczos}; +use crate::repro::{atan2, exp, ln, sin_cos}; // --------------------------------------------------------------------------- // Spin models @@ -118,6 +159,45 @@ pub fn square_edges(lx: usize, ly: usize, periodic: bool) -> Vec<(usize, usize)> e } +/// The next-nearest-neighbour bonds of an open (or periodic) chain. +pub fn chain_next_nearest(n: usize, periodic: bool) -> Vec<(usize, usize)> { + let mut e: Vec<(usize, usize)> = (0..n.saturating_sub(2)).map(|i| (i, i + 2)).collect(); + if periodic && n > 4 { + e.push((n - 2, 0)); + e.push((n - 1, 1)); + } + e +} + +/// The diagonal (next-nearest) bonds of an `lx × ly` square lattice, open or +/// periodic in both directions. +pub fn square_diagonals(lx: usize, ly: usize, periodic: bool) -> Vec<(usize, usize)> { + let mut e = Vec::new(); + let wrap = |v: isize, l: usize| -> Option { + if v >= 0 && (v as usize) < l { + Some(v as usize) + } else if periodic && l > 2 { + Some(v.rem_euclid(l as isize) as usize) + } else { + None + } + }; + for y in 0..ly { + for x in 0..lx { + let i = y * lx + x; + for dx in [1isize, -1] { + if let (Some(xx), Some(yy)) = (wrap(x as isize + dx, lx), wrap(y as isize + 1, ly)) { + let j = yy * lx + xx; + if i != j && !e.contains(&(i, j)) && !e.contains(&(j, i)) { + e.push((i, j)); + } + } + } + } + } + e +} + /// The checkerboard sublattice of an `lx × ly` lattice (a chain is `ly = 1`). pub fn checkerboard(lx: usize, ly: usize) -> Vec { (0..lx * ly).map(|i| (i % lx + i / lx).is_multiple_of(2)).collect() @@ -129,6 +209,16 @@ impl SpinModel { SpinModel { n, exchange: edges.iter().map(|&(a, b)| (a, b, j / 4.0)).collect(), zz: edges.iter().map(|&(a, b)| (a, b, j / 4.0)).collect(), x: Vec::new(), z: Vec::new() } } + /// The frustrated `J₁ Σ_nn Sᵢ·Sⱼ + J₂ Σ_nnn Sᵢ·Sⱼ`. With `J₂ > 0` no + /// sign rule is known in general, and a real network cannot follow it. + pub fn j1j2(n: usize, nn: &[(usize, usize)], nnn: &[(usize, usize)], j1: f64, j2: f64) -> SpinModel { + let mut m = SpinModel::heisenberg(n, nn, j1); + let second = SpinModel::heisenberg(n, nnn, j2); + m.exchange.extend(second.exchange); + m.zz.extend(second.zz); + m + } + /// The transverse-field Ising model `−J Σ σᶻσᶻ − h Σ σˣ` on `edges`. pub fn ising(n: usize, edges: &[(usize, usize)], j: f64, h: f64) -> SpinModel { SpinModel { n, exchange: Vec::new(), zz: edges.iter().map(|&(a, b)| (a, b, -j)).collect(), x: (0..n).map(|i| (i, -h)).collect(), z: Vec::new() } @@ -207,7 +297,7 @@ impl SpinModel { return None; } let apply = |inp: &[f64], out: &mut [f64]| self.apply_dense(inp, out); - Some(lanczos(&apply, &start, 60, 30, 1e-9).0) + Some(lanczos(&apply, &start, 60, 30, 1e-9, false).0) } } @@ -215,31 +305,96 @@ impl SpinModel { // The restricted Boltzmann machine // --------------------------------------------------------------------------- -/// `ln cosh x`, stable for large `|x|`. -fn lncosh(x: f64) -> f64 { - let a = x.abs(); - a + ln(1.0 + exp(-2.0 * a)) - core::f64::consts::LN_2 +/// The numbers a network's parameters are: real, or complex so the network +/// can carry a sign structure nobody supplied. +pub trait Amplitude: Scalar { + fn exp(self) -> Self; + /// `ln cosh`, stable for large real parts (shifted by `ln 2`). + fn lncosh(self) -> Self; + fn tanh(self) -> Self; + /// `self / d` for a real `d`. + fn div_re(self, d: f64) -> Self; + /// An initial weight, uniform in `±scale` (each part, when complex). + fn draw(uniform: &mut dyn FnMut() -> f64, scale: f64) -> Self; } -fn tanh(x: f64) -> f64 { - let e = exp(-2.0 * x.abs()); - let t = (1.0 - e) / (1.0 + e); - if x < 0.0 { -t } else { t } +impl Amplitude for f64 { + fn exp(self) -> f64 { + exp(self) + } + fn lncosh(self) -> f64 { + let a = self.abs(); + a + ln(1.0 + exp(-2.0 * a)) - core::f64::consts::LN_2 + } + fn tanh(self) -> f64 { + let e = exp(-2.0 * self.abs()); + let t = (1.0 - e) / (1.0 + e); + if self < 0.0 { -t } else { t } + } + fn div_re(self, d: f64) -> f64 { + self / d + } + fn draw(uniform: &mut dyn FnMut() -> f64, scale: f64) -> f64 { + scale * (2.0 * uniform() - 1.0) + } } -/// A real restricted Boltzmann machine: `ψ(s) = sign(s) · exp(Σ aᵢsᵢ) · -/// Π_j cosh(b_j + Σ W_jᵢ sᵢ)`, with `m = α·n` hidden units. On a bipartite +fn c_ln(z: C) -> C { + C { re: ln((z.re * z.re + z.im * z.im).sqrt()), im: atan2(z.im, z.re) } +} + +fn c_div(a: C, b: C) -> C { + let d = b.re * b.re + b.im * b.im; + C { re: (a.re * b.re + a.im * b.im) / d, im: (a.im * b.re - a.re * b.im) / d } +} + +impl Amplitude for C { + fn exp(self) -> C { + let m = exp(self.re); + let (s, c) = sin_cos(self.im); + C { re: m * c, im: m * s } + } + fn lncosh(self) -> C { + // cosh is even; with Re z ≥ 0, ln cosh z = z + ln(1 + e^{−2z}) − ln 2. + if self.re < 0.0 { + return Scalar::scale(self, -1.0).lncosh(); + } + let e = Scalar::scale(self, -2.0).exp(); + let l = c_ln(C { re: 1.0 + e.re, im: e.im }); + C { re: self.re + l.re - core::f64::consts::LN_2, im: self.im + l.im } + } + fn tanh(self) -> C { + if self.re < 0.0 { + return Scalar::scale(Scalar::scale(self, -1.0).tanh(), -1.0); + } + let e = Scalar::scale(self, -2.0).exp(); + c_div(C { re: 1.0 - e.re, im: -e.im }, C { re: 1.0 + e.re, im: e.im }) + } + fn div_re(self, d: f64) -> C { + C { re: self.re / d, im: self.im / d } + } + fn draw(uniform: &mut dyn FnMut() -> f64, scale: f64) -> C { + let re = scale * (2.0 * uniform() - 1.0); + let im = scale * (2.0 * uniform() - 1.0); + C { re, im } + } +} + +/// A restricted Boltzmann machine: `ψ(s) = sign(s) · exp(Σ aᵢsᵢ) · +/// Π_j cosh(b_j + Σ W_jᵢ sᵢ)`, with `m = α·n` hidden units. Its parameters are +/// real (`Rbm`), so `ψ` is positive up to the optional sign, or complex +/// (`Rbm`), so `ψ` carries a phase the training finds. On a bipartite /// lattice the optional sign is Marshall's, `(−1)` to the number of up spins /// on sublattice A, which makes the Heisenberg ground state positive in the /// rest. #[derive(Clone, Debug, PartialEq)] -pub struct Rbm { +pub struct Rbm { pub n: usize, pub m: usize, - pub a: Vec, - pub b: Vec, + pub a: Vec, + pub b: Vec, /// `W[j·n + i]`. - pub w: Vec, + pub w: Vec, pub sublattice: Option>, } @@ -266,50 +421,75 @@ impl Rng { } } -impl Rbm { - /// `α·n` hidden units, weights uniform in `±scale`, biases zero. +impl Rbm { + /// A real machine: `α·n` hidden units, weights uniform in `±scale`, + /// biases zero. pub fn new(n: usize, alpha: usize, scale: f64, seed: u64) -> Rbm { + Rbm::random(n, alpha, scale, seed) + } +} + +impl Rbm { + /// A complex machine: weights with real and imaginary parts uniform in + /// `±scale`, biases zero. + pub fn complex(n: usize, alpha: usize, scale: f64, seed: u64) -> Rbm { + Rbm::random(n, alpha, scale, seed) + } +} + +impl Rbm { + fn random(n: usize, alpha: usize, scale: f64, seed: u64) -> Rbm { let m = alpha.max(1) * n; let mut rng = Rng(seed); - let w = (0..m * n).map(|_| scale * (2.0 * rng.uniform() - 1.0)).collect(); - Rbm { n, m, a: vec![0.0; n], b: vec![0.0; m], w, sublattice: None } + let mut u = || rng.uniform(); + let w = (0..m * n).map(|_| T::draw(&mut u, scale)).collect(); + Rbm { n, m, a: vec![T::ZERO; n], b: vec![T::ZERO; m], w, sublattice: None } } /// Multiply by Marshall's sign for this sublattice. - pub fn with_marshall_sign(mut self, sublattice: Vec) -> Rbm { + pub fn with_marshall_sign(mut self, sublattice: Vec) -> Rbm { assert_eq!(sublattice.len(), self.n); self.sublattice = Some(sublattice); self } - /// The number of variational parameters. + /// The number of variational parameters (complex ones counted once). pub fn params(&self) -> usize { self.n + self.m + self.n * self.m } - fn theta(&self, s: &[i8]) -> Vec { + fn theta(&self, s: &[i8]) -> Vec { (0..self.m) .map(|j| { let row = &self.w[j * self.n..(j + 1) * self.n]; let mut t = self.b[j]; for (wi, &si) in row.iter().zip(s) { - t += wi * f64::from(si); + t = t.add(wi.scale(f64::from(si))); } t }) .collect() } - /// `ln |ψ(s)|` up to a constant. - pub fn log_amplitude(&self, s: &[i8]) -> f64 { - let mut l: f64 = self.a.iter().zip(s).map(|(a, &si)| a * f64::from(si)).sum(); + /// `ln ψ(s)` up to a constant, without the sign: real part `ln |ψ|`, + /// imaginary part the phase. + pub fn log_psi(&self, s: &[i8]) -> T { + let mut l = T::ZERO; + for (a, &si) in self.a.iter().zip(s) { + l = l.add(a.scale(f64::from(si))); + } for t in self.theta(s) { - l += lncosh(t); + l = l.add(t.lncosh()); } l } - /// The sign of `ψ(s)`. + /// `ln |ψ(s)|` up to a constant. + pub fn log_amplitude(&self, s: &[i8]) -> f64 { + self.log_psi(s).re() + } + + /// The sign of `ψ(s)` from Marshall's rule (one without it). pub fn sign(&self, s: &[i8]) -> f64 { match &self.sublattice { Some(sub) => { @@ -321,72 +501,72 @@ impl Rbm { } /// `ψ(s′)/ψ(s)` where `s′` is `s` with `flips` flipped, given `θ(s)`. - fn ratio(&self, s: &[i8], theta: &[f64], flips: &[usize]) -> f64 { - let mut l = 0.0; + fn ratio(&self, s: &[i8], theta: &[T], flips: &[usize]) -> T { + let mut l = T::ZERO; for &i in flips { - l -= 2.0 * self.a[i] * f64::from(s[i]); + l = l.sub(self.a[i].scale(2.0).scale(f64::from(s[i]))); } for (j, &t) in theta.iter().enumerate() { - let mut d = 0.0; + let mut d = T::ZERO; for &i in flips { - d -= 2.0 * self.w[j * self.n + i] * f64::from(s[i]); + d = d.sub(self.w[j * self.n + i].scale(2.0).scale(f64::from(s[i]))); } - l += lncosh(t + d) - lncosh(t); + l = l.add(t.add(d).lncosh().sub(t.lncosh())); } - let mut r = exp(l); + let mut r = l.exp(); if let Some(sub) = &self.sublattice && flips.iter().filter(|&&i| sub[i]).count() % 2 == 1 { - r = -r; + r = r.scale(-1.0); } r } /// `∂ ln ψ / ∂(a, b, W)` at `s`, given `θ(s)`. - fn log_derivatives(&self, s: &[i8], theta: &[f64], out: &mut [f64]) { + fn log_derivatives(&self, s: &[i8], theta: &[T], out: &mut [T]) { let (n, m) = (self.n, self.m); for (o, &si) in out[..n].iter_mut().zip(s) { - *o = f64::from(si); + *o = T::ONE.scale(f64::from(si)); } for (j, &th) in theta.iter().enumerate() { - let t = tanh(th); + let t = th.tanh(); out[n + j] = t; for (o, &si) in out[n + m + j * n..n + m + (j + 1) * n].iter_mut().zip(s) { - *o = t * f64::from(si); + *o = t.scale(f64::from(si)); } } } - fn shift(&mut self, delta: &[f64], step: f64) { + fn shift(&mut self, delta: &[T], step: f64) { let (n, m) = (self.n, self.m); for (x, d) in self.a.iter_mut().zip(&delta[..n]) { - *x -= step * d; + *x = x.sub(d.scale(step)); } for (x, d) in self.b.iter_mut().zip(&delta[n..n + m]) { - *x -= step * d; + *x = x.sub(d.scale(step)); } for (x, d) in self.w.iter_mut().zip(&delta[n + m..]) { - *x -= step * d; + *x = x.sub(d.scale(step)); } } } /// `E_loc(s) = Σ_{s′} H_{ss′} ψ(s′)/ψ(s)`. -fn local_energy(model: &SpinModel, rbm: &Rbm, s: &[i8], theta: &[f64]) -> f64 { - let mut e = 0.0; +fn local_energy(model: &SpinModel, rbm: &Rbm, s: &[i8], theta: &[T]) -> T { + let mut e = T::ZERO; for &(i, j, c) in &model.zz { - e += c * f64::from(s[i]) * f64::from(s[j]); + e = e.add(T::ONE.scale(c * f64::from(s[i]) * f64::from(s[j]))); } for &(i, c) in &model.z { - e += c * f64::from(s[i]); + e = e.add(T::ONE.scale(c * f64::from(s[i]))); } for &(i, j, c) in &model.exchange { if s[i] != s[j] { - e += 2.0 * c * rbm.ratio(s, theta, &[i, j]); + e = e.add(rbm.ratio(s, theta, &[i, j]).scale(2.0 * c)); } } for &(i, c) in &model.x { - e += c * rbm.ratio(s, theta, &[i]); + e = e.add(rbm.ratio(s, theta, &[i]).scale(c)); } e } @@ -430,9 +610,9 @@ impl Default for VmcConfig { /// One iteration's estimate. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Estimate { - /// `⟨E_loc⟩`. + /// `Re ⟨E_loc⟩`. pub energy: f64, - /// `⟨(E_loc − E)²⟩`: zero exactly at an eigenstate. + /// `⟨|E_loc − E|²⟩`: zero exactly at an eigenstate. pub variance: f64, /// The statistical error of `energy`, from the variance and the sample /// count, ignoring autocorrelation (zero when exact). @@ -441,26 +621,27 @@ pub struct Estimate { pub acceptance: f64, } -/// Variational Monte Carlo with stochastic reconfiguration on an RBM. +/// Variational Monte Carlo with stochastic reconfiguration on an RBM, real +/// or complex. #[derive(Clone, Debug, PartialEq)] -pub struct Vmc { +pub struct Vmc { pub model: SpinModel, - pub rbm: Rbm, + pub rbm: Rbm, pub cfg: VmcConfig, chains: Vec<(Vec, Rng)>, } -/// Samples: configurations, their weights (summing to one), local energies -/// and log-derivatives (`params` per row). -struct Batch { +/// Samples: their weights (summing to one), local energies and +/// log-derivatives (`params` per row). +struct Batch { weights: Vec, - energies: Vec, - derivs: Vec, + energies: Vec, + derivs: Vec, acceptance: f64, } -impl Vmc { - pub fn new(model: SpinModel, rbm: Rbm, cfg: VmcConfig) -> Vmc { +impl Vmc { + pub fn new(model: SpinModel, rbm: Rbm, cfg: VmcConfig) -> Vmc { assert_eq!(model.n, rbm.n, "the model and the machine differ in size"); if cfg.magnetisation.is_some() { assert!(model.conserves_magnetisation(), "the model does not conserve the magnetisation"); @@ -484,22 +665,22 @@ impl Vmc { v } - fn batch(&mut self) -> Batch { + fn batch(&mut self) -> Batch { if self.cfg.exact { self.enumerate() } else { self.sample() } } - fn enumerate(&self) -> Batch { + fn enumerate(&self) -> Batch { let n = self.model.n; assert!(n <= 24, "exact enumeration takes at most 24 spins"); let configs: Vec = (0..1usize << n).filter(|&x| self.cfg.magnetisation.is_none_or(|m| n as i32 - 2 * x.count_ones() as i32 == m)).collect(); let p = self.rbm.params(); let k = configs.len(); let mut logs = vec![0.0; k]; - let mut energies = vec![0.0; k]; - let mut derivs = vec![0.0; k * p]; + let mut energies = vec![T::ZERO; k]; + let mut derivs = vec![T::ZERO; k * p]; let threads = self.cfg.threads.max(1).min(k.div_ceil(64)).max(1); let per = k.div_ceil(threads); - let work = |range: core::ops::Range, logs: &mut [f64], energies: &mut [f64], derivs: &mut [f64]| { + let work = |range: core::ops::Range, logs: &mut [f64], energies: &mut [T], derivs: &mut [T]| { let start = range.start; for idx in range { let x = configs[idx]; @@ -516,8 +697,8 @@ impl Vmc { } else { std::thread::scope(|scope| { let mut lrest: &mut [f64] = &mut logs; - let mut erest: &mut [f64] = &mut energies; - let mut drest: &mut [f64] = &mut derivs; + let mut erest: &mut [T] = &mut energies; + let mut drest: &mut [T] = &mut derivs; for t in 0..threads { let lo = t * per; let hi = ((t + 1) * per).min(k); @@ -542,7 +723,7 @@ impl Vmc { Batch { weights: raw.iter().map(|r| r / z).collect(), energies, derivs, acceptance: 1.0 } } - fn sample(&mut self) -> Batch { + fn sample(&mut self) -> Batch { let n = self.model.n; let p = self.rbm.params(); let nc = self.chains.len(); @@ -551,10 +732,10 @@ impl Vmc { let thin = self.cfg.thin.max(1) * n; let rbm = &self.rbm; let model = &self.model; - let run = |chain: &mut (Vec, Rng)| -> (Vec, Vec, f64) { + let run = |chain: &mut (Vec, Rng)| -> (Vec, Vec, f64) { let (s, rng) = chain; let mut energies = Vec::with_capacity(per_chain); - let mut derivs = vec![0.0; per_chain * p]; + let mut derivs = vec![T::ZERO; per_chain * p]; let mut accepted = 0usize; for k in 0..per_chain { accepted += metropolis(rbm, s, rng, thin, exchange); @@ -565,7 +746,7 @@ impl Vmc { (energies, derivs, accepted as f64 / (per_chain * thin) as f64) }; let threads = self.cfg.threads.max(1).min(nc); - let results: Vec<(Vec, Vec, f64)> = if threads <= 1 || cfg!(target_arch = "wasm32") { + let results: Vec<(Vec, Vec, f64)> = if threads <= 1 || cfg!(target_arch = "wasm32") { self.chains.iter_mut().map(run).collect() } else { let per = nc.div_ceil(threads); @@ -601,8 +782,9 @@ impl Vmc { /// One step of stochastic reconfiguration: `θ ← θ − η (S + λ)⁻¹ F`, with /// `S` the covariance of the log-derivatives and `F` their covariance - /// with the local energy. It is solved in parameter space or in sample - /// space, whichever is smaller; the two give the same step. + /// with the local energy (Hermitian ones, for a complex machine). It is + /// solved in parameter space or in sample space, whichever is smaller; + /// the two give the same step. pub fn step(&mut self) -> Estimate { let b = self.batch(); let est = summarise(&b, self.cfg.exact); @@ -612,130 +794,134 @@ impl Vmc { } } -fn summarise(b: &Batch, exact: bool) -> Estimate { - let energy: f64 = b.weights.iter().zip(&b.energies).map(|(w, e)| w * e).sum(); - let variance: f64 = b.weights.iter().zip(&b.energies).map(|(w, e)| w * (e - energy) * (e - energy)).sum(); +fn summarise(b: &Batch, exact: bool) -> Estimate { + let energy: f64 = b.weights.iter().zip(&b.energies).map(|(w, e)| w * e.re()).sum(); + let et = T::ONE.scale(energy); + let variance: f64 = b + .weights + .iter() + .zip(&b.energies) + .map(|(&w, e)| { + let d = e.sub(et); + d.scale(w).mul(d.conj()).re() + }) + .sum(); let error = if exact { 0.0 } else { (variance / b.energies.len() as f64).sqrt() }; Estimate { energy, variance, error, acceptance: b.acceptance } } /// `(S + λ)⁻¹ F` from a batch. -fn sr_direction(b: &Batch, p: usize, shift: f64) -> Vec { +fn sr_direction(b: &Batch, p: usize, shift: f64) -> Vec { let k = b.weights.len(); - let energy: f64 = b.weights.iter().zip(&b.energies).map(|(w, e)| w * e).sum(); - let mut mean = vec![0.0; p]; + let energy: f64 = b.weights.iter().zip(&b.energies).map(|(w, e)| w * e.re()).sum(); + let et = T::ONE.scale(energy); + let mut mean = vec![T::ZERO; p]; for (r, &w) in b.weights.iter().enumerate() { for (m, &d) in mean.iter_mut().zip(&b.derivs[r * p..(r + 1) * p]) { - *m += w * d; + *m = m.add(d.scale(w)); } } // Rows √w (O − Ō) and √w (E − Ē). - let mut o = vec![0.0; k * p]; - let mut eps = vec![0.0; k]; + let mut o = vec![T::ZERO; k * p]; + let mut eps = vec![T::ZERO; k]; for r in 0..k { let sw = b.weights[r].sqrt(); for c in 0..p { - o[r * p + c] = sw * (b.derivs[r * p + c] - mean[c]); + o[r * p + c] = b.derivs[r * p + c].sub(mean[c]).scale(sw); } - eps[r] = sw * (b.energies[r] - energy); + eps[r] = b.energies[r].sub(et).scale(sw); } if p <= k { sr_parameter_space(&o, &eps, k, p, shift) } else { sr_sample_space(&o, eps, k, p, shift) } } -/// `(OᵀO + λ)⁻¹ Oᵀε`. -fn sr_parameter_space(o: &[f64], eps: &[f64], k: usize, p: usize, shift: f64) -> Vec { - { - // (OᵀO + λ) δ = Oᵀε. - let mut s = vec![0.0; p * p]; - for r in 0..k { - let row = &o[r * p..(r + 1) * p]; - for i in 0..p { - let oi = row[i]; - if oi == 0.0 { - continue; - } - for j in i..p { - s[i * p + j] += oi * row[j]; - } - } - } +/// `(O†O + λ)⁻¹ O†ε`. +fn sr_parameter_space(o: &[T], eps: &[T], k: usize, p: usize, shift: f64) -> Vec { + let mut s = vec![T::ZERO; p * p]; + for r in 0..k { + let row = &o[r * p..(r + 1) * p]; for i in 0..p { - s[i * p + i] += shift; - for j in 0..i { - s[i * p + j] = s[j * p + i]; + let oi = row[i].conj(); + if oi.is_zero() { + continue; + } + for j in i..p { + s[i * p + j] = s[i * p + j].add(oi.mul(row[j])); } } - let mut f = vec![0.0; p]; - for r in 0..k { - for c in 0..p { - f[c] += o[r * p + c] * eps[r]; - } - } - cholesky_solve(&mut s, &mut f, p); - f } + for i in 0..p { + s[i * p + i] = s[i * p + i].add(T::ONE.scale(shift)); + for j in 0..i { + s[i * p + j] = s[j * p + i].conj(); + } + } + let mut f = vec![T::ZERO; p]; + for r in 0..k { + for c in 0..p { + f[c] = f[c].add(o[r * p + c].conj().mul(eps[r])); + } + } + cholesky_solve(&mut s, &mut f, p); + f } -/// `Oᵀ (OOᵀ + λ)⁻¹ ε`. -fn sr_sample_space(o: &[f64], eps: Vec, k: usize, p: usize, shift: f64) -> Vec { - { - // δ = Oᵀ (OOᵀ + λ)⁻¹ ε. - let mut t = vec![0.0; k * k]; - for i in 0..k { - for j in i..k { - let mut acc = 0.0; - for c in 0..p { - acc += o[i * p + c] * o[j * p + c]; - } - t[i * k + j] = acc; - t[j * k + i] = acc; - } - t[i * k + i] += shift; - } - let mut x = eps; - cholesky_solve(&mut t, &mut x, k); - let mut delta = vec![0.0; p]; - for r in 0..k { +/// `O† (OO† + λ)⁻¹ ε`. +fn sr_sample_space(o: &[T], eps: Vec, k: usize, p: usize, shift: f64) -> Vec { + let mut t = vec![T::ZERO; k * k]; + for i in 0..k { + for j in i..k { + let mut acc = T::ZERO; for c in 0..p { - delta[c] += o[r * p + c] * x[r]; + acc = acc.add(o[i * p + c].mul(o[j * p + c].conj())); } + t[i * k + j] = acc; + t[j * k + i] = acc.conj(); } - delta + t[i * k + i] = t[i * k + i].add(T::ONE.scale(shift)); } + let mut x = eps; + cholesky_solve(&mut t, &mut x, k); + let mut delta = vec![T::ZERO; p]; + for r in 0..k { + for c in 0..p { + delta[c] = delta[c].add(o[r * p + c].conj().mul(x[r])); + } + } + delta } -/// Solve `A x = b` for symmetric positive-definite `A` (overwritten by its -/// Cholesky factor), `b` overwritten by `x`. +/// Solve `A x = b` for Hermitian positive-definite `A` (overwritten by its +/// Cholesky factor `L`, `A = L L†`), `b` overwritten by `x`. #[allow(clippy::needless_range_loop)] -fn cholesky_solve(a: &mut [f64], b: &mut [f64], n: usize) { +fn cholesky_solve(a: &mut [T], b: &mut [T], n: usize) { for j in 0..n { - let mut d = a[j * n + j]; + let mut d = a[j * n + j].re(); for k in 0..j { - d -= a[j * n + k] * a[j * n + k]; + d -= a[j * n + k].norm2(); } let d = d.max(1e-300).sqrt(); - a[j * n + j] = d; + a[j * n + j] = T::ONE.scale(d); for i in j + 1..n { let mut s = a[i * n + j]; for k in 0..j { - s -= a[i * n + k] * a[j * n + k]; + s = s.sub(a[i * n + k].mul(a[j * n + k].conj())); } - a[i * n + j] = s / d; + a[i * n + j] = s.div_re(d); } } for i in 0..n { let mut s = b[i]; for k in 0..i { - s -= a[i * n + k] * b[k]; + s = s.sub(a[i * n + k].mul(b[k])); } - b[i] = s / a[i * n + i]; + b[i] = s.div_re(a[i * n + i].re()); } for i in (0..n).rev() { let mut s = b[i]; for k in i + 1..n { - s -= a[k * n + i] * b[k]; + s = s.sub(a[k * n + i].conj().mul(b[k])); } - b[i] = s / a[i * n + i]; + b[i] = s.div_re(a[i * n + i].re()); } } @@ -759,7 +945,7 @@ fn start_configuration(n: usize, magnetisation: Option, rng: &mut Rng) -> V /// `moves` Metropolis moves on `s` (single flips, or swaps of antiparallel /// spins), accepted with probability `|ψ(s′)/ψ(s)|²`. Returns the number /// accepted. -fn metropolis(rbm: &Rbm, s: &mut [i8], rng: &mut Rng, moves: usize, exchange: bool) -> usize { +fn metropolis(rbm: &Rbm, s: &mut [i8], rng: &mut Rng, moves: usize, exchange: bool) -> usize { let n = s.len(); let mut theta = rbm.theta(s); let mut accepted = 0; @@ -775,11 +961,11 @@ fn metropolis(rbm: &Rbm, s: &mut [i8], rng: &mut Rng, moves: usize, exchange: bo vec![rng.below(n)] }; let r = rbm.ratio(s, &theta, &flips); - if rng.uniform() < r * r { + if rng.uniform() < r.norm2() { for &i in &flips { let ds = -2.0 * f64::from(s[i]); for (j, t) in theta.iter_mut().enumerate() { - *t += rbm.w[j * n + i] * ds; + *t = t.add(rbm.w[j * n + i].scale(ds)); } s[i] = -s[i]; } @@ -918,5 +1104,112 @@ mod tests { assert_eq!(one.rbm, many.rbm); assert_eq!(one.chains, many.chains); } -} + #[test] + fn a_complex_machine_with_real_parameters_is_the_real_one() { + let n = 6; + let real = Rbm::new(n, 2, 0.4, 5); + let cplx = Rbm { n: real.n, m: real.m, a: real.a.iter().map(|&re| C { re, im: 0.0 }).collect(), b: real.b.iter().map(|&re| C { re, im: 0.0 }).collect(), w: real.w.iter().map(|&re| C { re, im: 0.0 }).collect(), sublattice: None }; + let mut rng = Rng(8); + for _ in 0..10 { + let s = start_configuration(n, None, &mut rng); + let (tr, tc) = (real.theta(&s), cplx.theta(&s)); + assert!((real.log_amplitude(&s) - cplx.log_amplitude(&s)).abs() < 1e-14); + assert!(cplx.log_psi(&s).im.abs() < 1e-14); + for flips in [vec![1], vec![0, 4]] { + let (r, c) = (real.ratio(&s, &tr, &flips), cplx.ratio(&s, &tc, &flips)); + assert!((r - c.re).abs() < 1e-13 * r.abs().max(1.0) && c.im.abs() < 1e-13); + } + } + } + + #[test] + fn complex_ratios_and_derivatives_are_the_machines() { + let n = 6; + let mut rbm = Rbm::complex(n, 2, 0.4, 6); + let mut rng = Rng(12); + for v in rbm.a.iter_mut().chain(rbm.b.iter_mut()) { + *v = C { re: 0.3 * (rng.uniform() - 0.5), im: 0.3 * (rng.uniform() - 0.5) }; + } + let p = rbm.params(); + for _ in 0..10 { + let s = start_configuration(n, None, &mut rng); + let theta = rbm.theta(&s); + for flips in [vec![rng.below(n)], vec![1, 5]] { + let mut t = s.clone(); + for &i in &flips { + t[i] = -t[i]; + } + let want = rbm.log_psi(&t).sub(rbm.log_psi(&s)).exp(); + let got = rbm.ratio(&s, &theta, &flips); + assert!(got.sub(want).norm2().sqrt() < 1e-12 * want.norm2().sqrt().max(1.0), "{got:?} vs {want:?}"); + } + // Holomorphic: the derivative along a real step is the derivative. + let mut d = vec![C::ZERO; p]; + rbm.log_derivatives(&s, &theta, &mut d); + let h = 1e-6; + for k in [0, n + 1, n + rbm.m + 3, p - 1] { + let mut e = vec![C::ZERO; p]; + e[k] = C { re: -h, im: 0.0 }; + let mut up = rbm.clone(); + up.shift(&e, 1.0); + e[k] = C { re: h, im: 0.0 }; + let mut dn = rbm.clone(); + dn.shift(&e, 1.0); + let fd = up.log_psi(&s).sub(dn.log_psi(&s)).scale(0.5 / h); + assert!(fd.sub(d[k]).norm2().sqrt() < 1e-7, "param {k}: {fd:?} vs {:?}", d[k]); + } + } + } + + #[test] + fn both_spaces_take_the_same_complex_step() { + let mut rng = Rng(19); + for (k, p) in [(12usize, 5usize), (5, 12)] { + let o: Vec = (0..k * p).map(|_| C { re: rng.uniform() - 0.5, im: rng.uniform() - 0.5 }).collect(); + let eps: Vec = (0..k).map(|_| C { re: rng.uniform() - 0.5, im: rng.uniform() - 0.5 }).collect(); + let a = sr_parameter_space(&o, &eps, k, p, 1e-3); + let b = sr_sample_space(&o, eps.clone(), k, p, 1e-3); + for (x, y) in a.iter().zip(&b) { + assert!(x.sub(*y).norm2().sqrt() < 1e-9 * x.norm2().sqrt().max(1.0), "{x:?} vs {y:?}"); + } + } + } + + #[test] + fn exact_diagonalisation_solves_majumdar_ghosh() { + // At J₂ = J₁/2 a ring's ground states are the two dimer coverings, + // at exactly −3n/8. + for n in [8usize, 12] { + let m = SpinModel::j1j2(n, &chain_edges(n, true), &chain_next_nearest(n, true), 1.0, 0.5); + let e = m.exact_ground_energy(Some(0)).unwrap(); + assert!((e + 3.0 * n as f64 / 8.0).abs() < 1e-9, "n={n}: {e}"); + } + } + + #[test] + fn a_complex_machine_learns_a_sign_a_real_one_cannot() { + // The Heisenberg ring and the Majumdar–Ghosh ring, with no sign rule + // given. A real network is positive and stays far off; a complex one + // finds the signs. + let n = 8; + let models = [ + SpinModel::heisenberg(n, &chain_edges(n, true), 1.0), + SpinModel::j1j2(n, &chain_edges(n, true), &chain_next_nearest(n, true), 1.0, 0.5), + ]; + for model in models { + let exact = model.exact_ground_energy(Some(0)).unwrap(); + let cfg = VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, shift: 1e-3, ..VmcConfig::default() }; + let mut real = Vmc::new(model.clone(), Rbm::new(n, 2, 0.05, 1), cfg); + let mut cplx = Vmc::new(model, Rbm::complex(n, 2, 0.05, 3), cfg); + let (mut er, mut ec) = (real.step(), cplx.step()); + for _ in 0..250 { + er = real.step(); + ec = cplx.step(); + } + assert!((er.energy - exact) / exact.abs() > 0.1, "real: {} vs {exact}", er.energy); + assert!(ec.energy >= exact - 1e-9, "variational: {} vs {exact}", ec.energy); + assert!((ec.energy - exact) / exact.abs() < 2e-3, "complex: {} vs {exact}", ec.energy); + } + } +} diff --git a/wai/wai-quantum/src/repro.rs b/wai/wai-quantum/src/repro.rs index 36bf82eb9..f0099dbf0 100644 --- a/wai/wai-quantum/src/repro.rs +++ b/wai/wai-quantum/src/repro.rs @@ -16,7 +16,8 @@ feature = "quantum_kernel", feature = "quantum_tn", feature = "quantum_sv", - feature = "quantum_tdvp" + feature = "quantum_tdvp", + feature = "quantum_nqs" ))] mod trig; #[cfg(any( @@ -28,7 +29,8 @@ mod trig; feature = "quantum_kernel", feature = "quantum_tn", feature = "quantum_sv", - feature = "quantum_tdvp" + feature = "quantum_tdvp", + feature = "quantum_nqs" ))] pub use trig::sin_cos; @@ -169,21 +171,21 @@ pub fn exp(x: f64) -> f64 { // Arctangent // --------------------------------------------------------------------------- -#[cfg(feature = "quantum_qfhrr")] +#[cfg(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))] const ATANHI: [f64; 4] = [ 4.636_476_090_008_060_935_15e-1, 7.853_981_633_974_482_789_99e-1, 9.827_937_232_473_290_540_82e-1, 1.570_796_326_794_896_558_00, ]; -#[cfg(feature = "quantum_qfhrr")] +#[cfg(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))] const ATANLO: [f64; 4] = [ 2.269_877_745_296_168_709_24e-17, 3.061_616_997_868_383_017_93e-17, 1.390_331_103_123_099_845_16e-17, 6.123_233_995_736_766_035_87e-17, ]; -#[cfg(feature = "quantum_qfhrr")] +#[cfg(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))] const AT: [f64; 11] = [ 3.333_333_333_333_293_180_27e-1, -1.999_999_999_987_648_324_76e-1, @@ -199,7 +201,7 @@ const AT: [f64; 11] = [ ]; /// The arctangent, within an ulp of the true value. -#[cfg(feature = "quantum_qfhrr")] +#[cfg(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))] pub fn atan(x: f64) -> f64 { let ix = (x.to_bits() >> 32) as u32; let negative = ix >> 31 != 0; @@ -250,14 +252,14 @@ pub fn atan(x: f64) -> f64 { if negative { -z } else { z } } -#[cfg(feature = "quantum_qfhrr")] +#[cfg(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))] const PI: f64 = 3.141_592_653_589_793_116; -#[cfg(feature = "quantum_qfhrr")] +#[cfg(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))] const PI_LO: f64 = 1.224_646_799_147_353_177_2e-16; /// The angle of the point `(x, y)`, in `(−π, π]`, within an ulp of the true /// value. -#[cfg(feature = "quantum_qfhrr")] +#[cfg(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))] pub fn atan2(y: f64, x: f64) -> f64 { if x.is_nan() || y.is_nan() { return x + y; diff --git a/wai/wai-quantum/tests/bit_identity.rs b/wai/wai-quantum/tests/bit_identity.rs index deed66bd0..0df3a85c8 100644 --- a/wai/wai-quantum/tests/bit_identity.rs +++ b/wai/wai-quantum/tests/bit_identity.rs @@ -475,7 +475,7 @@ fn dmrg_ground_states() { let (vals, vecs) = symmetric_eigen(&sym, 6); v.extend(vals); v.extend(vecs.iter().map(|x| x.abs())); - assert_eq!(digest(&v), 0xc0a7_f146_6ad9_b7e8); + assert_eq!(digest(&v), 0xcdd9_d9fb_5ba6_d71c); } #[test] @@ -546,6 +546,15 @@ fn neural_quantum_states() { v.extend([e.energy, e.variance, e.error, e.acceptance]); } v.extend(&sampled.rbm.w); - assert_eq!(digest(&v), 0x6d09_94ff_8d98_0dc1); + // A complex machine on the Majumdar–Ghosh ring, which has no sign rule: + // complex exponentials, logarithms and arctangents in every step. + let mg = SpinModel::j1j2(n, &chain_edges(n, true), &chain_next_nearest(n, true), 1.0, 0.5); + let mut cplx = Vmc::new(mg, Rbm::complex(n, 2, 0.05, 3), VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, threads: 1, ..VmcConfig::default() }); + for _ in 0..4 { + let e = cplx.step(); + v.extend([e.energy, e.variance]); + } + v.extend(cplx.rbm.w.iter().flat_map(|c| [c.re, c.im])); + assert_eq!(digest(&v), 0x0ede_030d_b0c7_54dd); } diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index c8aadced7..68a8339fe 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.36" +version = "0.3.37" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-rs/src/bin/wai_quantum.rs b/wai/wai-rs/src/bin/wai_quantum.rs index 8661fbd65..9b8910880 100644 --- a/wai/wai-rs/src/bin/wai_quantum.rs +++ b/wai/wai-rs/src/bin/wai_quantum.rs @@ -173,8 +173,8 @@ fn usage() -> ! { MOLECULES (DMRG from FCIDUMP integrals, full CI as referee):\n \ wai-quantum chem --fcidump FILE [--bond D] [--sweeps S] [--penalty P] [--no-fci] [--json]\n\n\ NEURAL QUANTUM STATES (RBM, variational Monte Carlo, stochastic reconfiguration):\n \ - wai-quantum nqs [--model heisenberg|ising] [--lx L] [--ly L] [--h H] [--alpha A] [--steps N]\n \ - [--samples S] [--lr R] [--exact] [--seed S] [--json]\n\n\ + wai-quantum nqs [--model heisenberg|j1j2|ising] [--lx L] [--ly L] [--j2 J] [--h H] [--periodic]\n \ + [--complex] [--alpha A] [--steps N] [--samples S] [--lr R] [--exact] [--seed S] [--json]\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 \ @@ -1170,36 +1170,55 @@ fn main() { if !(2..=400).contains(&n) { die("--lx × --ly: 2 to 400 spins"); } + let periodic = flag(args, "--periodic"); let kind = opt(args, "--model", "").unwrap_or_else(|| "heisenberg".into()); - let edges = if ly == 1 { chain_edges(n, false) } else { square_edges(lx, ly, false) }; - let heis = match kind.as_str() { - "heisenberg" => true, - "ising" => false, - _ => die("--model: heisenberg or ising"), + let nn = if ly == 1 { chain_edges(n, periodic) } else { square_edges(lx, ly, periodic) }; + let nnn = if ly == 1 { chain_next_nearest(n, periodic) } else { square_diagonals(lx, ly, periodic) }; + let model = match kind.as_str() { + "heisenberg" => SpinModel::heisenberg(n, &nn, 1.0), + "j1j2" => SpinModel::j1j2(n, &nn, &nnn, 1.0, num(args, "--j2", "", 0.5)), + "ising" => SpinModel::ising(n, &nn, 1.0, num(args, "--h", "", 1.0)), + _ => die("--model: heisenberg, j1j2 or ising"), }; - let model = if heis { SpinModel::heisenberg(n, &edges, 1.0) } else { SpinModel::ising(n, &edges, 1.0, num(args, "--h", "", 1.0)) }; - let sector = if heis { Some(0) } else { None }; + let sector = if model.conserves_magnetisation() { Some(0) } else { None }; let exact = flag(args, "--exact"); if exact && n > 20 { die("--exact enumerates every configuration: at most 20 spins"); } - let mut rbm = Rbm::new(n, num(args, "--alpha", "", 2), 0.01, num(args, "--seed", "", 7)); - if heis { - rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); - } + let complex = flag(args, "--complex"); + // Marshall's sign for the real network on the unfrustrated model; + // a complex network finds its own. + let marshall = kind == "heisenberg" && !complex; let cfg = VmcConfig { samples: num(args, "--samples", "", 1000), exact, magnetisation: sector, lr: num(args, "--lr", "", 0.02), seed: num(args, "--seed", "", 7), ..VmcConfig::default() }; let steps: usize = num(args, "--steps", "", 300); + let alpha: usize = num(args, "--alpha", "", 2); + let seed: u64 = num(args, "--seed", "", 7); let reference = if n <= 20 { model.exact_ground_energy(sector) } else { None }; - let params = rbm.params(); - let mut vmc = Vmc::new(model, rbm, cfg); - let mut last = vmc.estimate(); - for _ in 0..steps { - last = vmc.step(); + fn run(model: SpinModel, rbm: Rbm, cfg: VmcConfig, steps: usize) -> (usize, Estimate, Estimate) { + let params = rbm.params(); + let mut vmc = Vmc::new(model, rbm, cfg); + let mut last = vmc.estimate(); + for _ in 0..steps { + last = vmc.step(); + } + (params, last, vmc.estimate()) } - let fin = vmc.estimate(); + let (params, last, fin) = if complex { + let mut rbm = Rbm::complex(n, alpha, 0.01, seed); + if marshall { + rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); + } + run(model, rbm, cfg, steps) + } else { + let mut rbm = Rbm::new(n, alpha, 0.01, seed); + if marshall { + rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); + } + run(model, rbm, cfg, steps) + }; if flag(args, "--json") { println!( - "{{\"energy\":{},\"error\":{},\"exact\":{},\"last_step\":{},\"model\":\"{kind}\",\"parameters\":{params},\"sites\":{n},\"steps\":{steps},\"variance\":{}}}", + "{{\"complex\":{complex},\"energy\":{},\"error\":{},\"exact\":{},\"last_step\":{},\"model\":\"{kind}\",\"parameters\":{params},\"sites\":{n},\"steps\":{steps},\"variance\":{}}}", fin.energy, fin.error, reference.map_or("null".into(), |e| e.to_string()), @@ -1207,7 +1226,7 @@ fn main() { fin.variance ); } else { - println!("{kind}, {lx}×{ly}, RBM with {params} parameters, {steps} steps{}", if exact { ", exact expectations" } else { "" }); + println!("{kind}, {lx}×{ly}{}, {} RBM with {params} parameters, {steps} steps{}", if periodic { " periodic" } else { "" }, if complex { "complex" } else { "real" }, if exact { ", exact expectations" } else { "" }); println!("energy {:.8} ± {:.8} variance {:.2e}", fin.energy, fin.error, fin.variance); if let Some(e) = reference { println!("exact {e:.8} relative error {:.2e}", (fin.energy - e) / e.abs()); diff --git a/wai/wai-web/Cargo.lock b/wai/wai-web/Cargo.lock index 382807f5b..6b4cc55b2 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.36" +version = "0.3.37" dependencies = [ "blake3", "ed25519-dalek", -- 2.54.0 (Apple Git-157) From defa8a8a9422aea000b2c03e9344b7ec0432ffd5 Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Tue, 6 Oct 2026 17:27:34 -0400 Subject: [PATCH 2/3] wai-quantum 0.3.37: the 10-site trap is spin-flip parity, not momentum; drop a stray probe With zero biases the network is exactly even under a global spin flip, and training keeps it even. A singlet is odd when half the ring is odd, so on 10-site rings the ground state was out of reach. Drawn biases reach it, and the docs now say so, with both starts measured. The earlier text blamed momentum, which is wrong. Also removes examples/nqs_probe.rs, a scratch probe that the squash swept in. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum/examples/nqs_probe.rs | 20 -------------------- wai/wai-quantum/src/quantum_nqs.rs | 27 ++++++++++++++++++++------- 2 files changed, 20 insertions(+), 27 deletions(-) delete mode 100644 wai/wai-quantum/examples/nqs_probe.rs diff --git a/wai/wai-quantum/examples/nqs_probe.rs b/wai/wai-quantum/examples/nqs_probe.rs deleted file mode 100644 index 6a9766664..000000000 --- a/wai/wai-quantum/examples/nqs_probe.rs +++ /dev/null @@ -1,20 +0,0 @@ -use wai_quantum::quantum_nqs::*; -fn main() { - for n in [10usize, 12] { - let mg = SpinModel::j1j2(n, &chain_edges(n, true), &chain_next_nearest(n, true), 1.0, 0.5); - let heis = SpinModel::heisenberg(n, &chain_edges(n, true), 1.0); - let j1j2 = SpinModel::j1j2(n, &chain_edges(n, true), &chain_next_nearest(n, true), 1.0, 0.3); - for (name, model) in [("majumdar-ghosh", mg), ("heisenberg (no sign)", heis), ("j1j2 0.3", j1j2)] { - let exact = model.exact_ground_energy(Some(0)).unwrap(); - for scale in [0.05, 0.3] { - let cfg = VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, shift: 1e-3, ..VmcConfig::default() }; - let mut cplx = Vmc::new(model.clone(), Rbm::complex(n, 2, scale, 3), cfg); - let mut ec = cplx.step(); - for _ in 0..400 { - ec = cplx.step(); - } - println!("n={n} {name:22} scale {scale}: exact {exact:.6} complex {:.6} rel {:.2e} variance {:.1e}", ec.energy, (ec.energy - exact) / exact.abs(), ec.variance); - } - } - } -} diff --git a/wai/wai-quantum/src/quantum_nqs.rs b/wai/wai-quantum/src/quantum_nqs.rs index ddc0b73e4..13aa93695 100644 --- a/wai/wai-quantum/src/quantum_nqs.rs +++ b/wai/wai-quantum/src/quantum_nqs.rs @@ -75,11 +75,23 @@ //! //! On 10-site rings every model above settled about 9 % high, on a state //! with essentially zero variance: an eigenstate, but not the ground state. -//! When half the ring's length is odd, the ground state has momentum `π`. A -//! network started near the uniform state begins at momentum zero, and -//! training stays there and finds the lowest state it can reach. The -//! variance cannot catch this, since an excited eigenstate has none. Only a -//! referee can. +//! The cause is a symmetry the start imposes. +//! +//! - With its biases at zero, the network is exactly even under flipping +//! every spin. Training keeps it even, because the biases' gradients vanish +//! on an even state. +//! - A singlet is even or odd under that flip as half the ring's length is +//! even or odd. On 10-site rings the singlet ground state is therefore out +//! of reach, and training finds the lowest even state instead. +//! - Starting the biases away from zero frees it. The same 10-site runs then +//! reach `3·10⁻⁴` (Heisenberg), `1·10⁻⁶` (Majumdar–Ghosh) and `1.4·10⁻⁵` +//! (J₁–J₂ 0.3) in 400 steps. +//! - On 12-site rings, where the ground state is even, the even start is the +//! better one (the table above). Drawn biases there reach only `2·10⁻²`, +//! `6·10⁻⁷` and `3·10⁻³`. +//! +//! The variance cannot catch the trap, since an excited eigenstate has none. +//! Only a referee can. //! //! The same trap caught the referee once. Lanczos stopped a restart once its //! lowest Ritz pair had converged, which is safe from DMRG's warm starts. @@ -95,8 +107,9 @@ //! 600 steps of 1000 samples leave all three networks about 3 % above it: //! complex with Marshall's sign, complex without, and real with it. At //! this size and budget the network is the limit, not the sampling. -//! - **Symmetry is not imposed.** A network is not projected onto a momentum -//! or spin sector, so it can settle in the wrong one, as above. +//! - **Symmetry is not chosen.** A network is not projected onto a momentum +//! or spin-flip sector. Started with zero biases it stays even under a +//! global spin flip, and that settles the sector, as above. //! - **Error bars ignore autocorrelation.** The reported error is the //! standard error of uncorrelated samples. Chains thinned by two sweeps are //! correlated, so the true error is larger. -- 2.54.0 (Apple Git-157) From a6636aa3dcecfd5736f227470dd11a7cd5e6087d Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Tue, 6 Oct 2026 17:28:18 -0400 Subject: [PATCH 3/3] wai-quantum 0.3.37: nqs_ground declares its feature Without it a build of every example with default features failed on it (since 0.3.36). Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum/Cargo.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index a44d2ae1c..9a09ec8ed 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -146,6 +146,10 @@ required-features = ["quantum_qio"] name = "tdvp_quench" required-features = ["quantum_tdvp"] +[[example]] +name = "nqs_ground" +required-features = ["quantum_nqs"] + [[example]] name = "chem_dmrg" required-features = ["quantum_chem"] -- 2.54.0 (Apple Git-157)