diff --git a/wai/wai-pack/Cargo.lock b/wai/wai-pack/Cargo.lock index c2f9441f7..75709e7e7 100644 --- a/wai/wai-pack/Cargo.lock +++ b/wai/wai-pack/Cargo.lock @@ -682,7 +682,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.38" +version = "0.3.39" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index d973270ee..f98d230e3 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.38" +version = "0.3.39" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index dbbec2e3d..9bf4e3356 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.38" +version = "0.3.39" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index a35e90d40..2e7f95322 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "wai-quantum" -version = "0.3.38" +version = "0.3.39" edition = "2024" rust-version = "1.97" license = "Apache-2.0" diff --git a/wai/wai-quantum/examples/nqs_ground.rs b/wai/wai-quantum/examples/nqs_ground.rs index 8bf76217d..91d7573ca 100644 --- a/wai/wai-quantum/examples/nqs_ground.rs +++ b/wai/wai-quantum/examples/nqs_ground.rs @@ -5,12 +5,21 @@ //! ```text //! 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] -//! [--project [--kpi] [--odd]] [--dmrg-bond D] [--threads T] +//! [--complex] [--marshall] [--alpha A] [--steps N] [--samples S] [--lr R] [--shift λ] +//! [--init W] [--exact] [--guide ε] +//! [--project [--kpi] [--point-group a1|a2|b1|b2] [--odd]] [--dmrg-bond D] +//! [--threads T] [--exact-final] //! //! `--project` projects onto translations of the periodic lattice at -//! momentum 0 (or π with `--kpi`) and onto even (or `--odd`) parity under a -//! global spin flip. +//! momentum 0 (or π with `--kpi`), with `--point-group` on a square lattice +//! onto a representation of its rotations and reflections too, and onto even +//! (or `--odd`) parity under a global spin flip. Up to 20 spins, the +//! reference is then the exact ground energy in that sector. `--exact-final` +//! also evaluates the trained network with every configuration enumerated +//! (at most 24 spins), separating what training left from what sampling +//! adds. `--guide ε` samples `|ψ|²` with a floor `ε` under it, weighted +//! back (`VmcConfig::guide`); projected networks default to `0.1`, others to +//! none. //! ``` use wai_quantum::quantum_dmrg::{dmrg_mpo, DmrgConfig, Mps, C}; use wai_quantum::quantum_nqs::*; @@ -19,7 +28,7 @@ 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) { +fn train(model: SpinModel, rbm: Rbm, cfg: VmcConfig, steps: usize, reference: f64, exact_final: bool) { 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(); @@ -32,6 +41,13 @@ fn train(model: SpinModel, rbm: Rbm, cfg: VmcConfig, steps: usi vmc.cfg.samples *= 8; let f = vmc.estimate(); println!("final E {:.6} ± {:.6} variance {:.2e} relative error {:.2e}", f.energy, f.error, f.variance, (f.energy - reference) / reference.abs()); + // The trained network's own energy and variance, every configuration + // enumerated: what is left once sampling is taken out. + if exact_final { + vmc.cfg.exact = true; + let x = vmc.estimate(); + println!("exact E {:.8} variance {:.2e} relative error {:.2e} (every configuration)", x.energy, x.variance, (x.energy - reference) / reference.abs()); + } } fn main() { @@ -61,32 +77,64 @@ fn main() { }; 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 mut cfg = VmcConfig { samples: arg(&a, "--samples", 1000), chains: 16, magnetisation: sector, lr: arg(&a, "--lr", 0.02), shift: arg(&a, "--shift", 1e-3), exact: a.iter().any(|x| x == "--exact"), threads, ..VmcConfig::default() }; let steps: usize = arg(&a, "--steps", 600); let marshall = a.iter().any(|x| x == "--marshall"); + let exact_final = a.iter().any(|x| x == "--exact-final"); + let point = a.iter().position(|x| x == "--point-group").map(|i| match a.get(i + 1).map(String::as_str) { + Some("a1") => PointIrrep::A1, + Some("a2") => PointIrrep::A2, + Some("b1") => PointIrrep::B1, + Some("b2") => PointIrrep::B2, + other => panic!("--point-group: a1, a2, b1 or b2, not {other:?}"), + }); let symmetry = a.iter().any(|x| x == "--project").then(|| { assert!(periodic, "--project needs --periodic"); let kpi = a.iter().any(|x| x == "--kpi"); - let base = if ly == 1 { Symmetry::ring(n, kpi) } else { Symmetry::torus(lx, ly, kpi, kpi) }; + let base = match point { + Some(irrep) => { + assert_eq!(lx, ly, "--point-group needs a square lattice"); + Symmetry::square(lx, kpi, irrep) + } + None if ly == 1 => Symmetry::ring(n, kpi), + None => Symmetry::torus(lx, ly, kpi, kpi), + }; base.with_spin_flip(a.iter().any(|x| x == "--odd")) }); + // A projected network is held to its own sector's ground energy. + let reference = match &symmetry { + Some(sym) if n <= 20 => { + let e = model.exact_ground_energy_in(sector, sym).expect("the model is invariant under the group, and the sector is not empty"); + println!("in its sector: {e:.8}"); + e + } + _ => reference, + }; + // Parameters start uniform in ±init. A projected network that starts + // near zero is nearly the projection of a constant, the ferromagnet: an + // eigenstate, which reconfiguration barely moves. + let init: f64 = arg(&a, "--init", if symmetry.is_some() { 0.1 } else { 0.01 }); + // A projected network's amplitude is a sum that cancels near its nodes, + // where the local energy is large and |ψ|² rarely samples; a guide floor + // visits there, and training then corrects what it finds. + cfg.guide = arg(&a, "--guide", if symmetry.is_some() { 0.1 } else { 0.0 }); if a.iter().any(|x| x == "--complex") { - let mut rbm: Rbm = Rbm::complex(n, alpha, 0.01, 7); + let mut rbm: Rbm = Rbm::complex(n, alpha, init, 7); if marshall { rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); } if let Some(sym) = symmetry { rbm = rbm.with_symmetry(sym); } - train(model, rbm, cfg, steps, reference); + train(model, rbm, cfg, steps, reference, exact_final); } else { - let mut rbm = Rbm::new(n, alpha, 0.01, 7); + let mut rbm = Rbm::new(n, alpha, init, 7); if marshall { rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); } if let Some(sym) = symmetry { rbm = rbm.with_symmetry(sym); } - train(model, rbm, cfg, steps, reference); + train(model, rbm, cfg, steps, reference, exact_final); } } diff --git a/wai/wai-quantum/src/quantum_nqs.rs b/wai/wai-quantum/src/quantum_nqs.rs index e4c6ffc59..f2440f3f7 100644 --- a/wai/wai-quantum/src/quantum_nqs.rs +++ b/wai/wai-quantum/src/quantum_nqs.rs @@ -21,22 +21,35 @@ //! the network would be exactly even under a global spin flip. //! - **The sector.** [`Symmetry`] projects the network, //! `Ψ(s) = Σ_g χ(g) ψ(g·s)`, over a lattice's translations (momentum 0 or -//! π along each direction) and, optionally, the flip of every spin (even -//! or odd). Ratios and log-derivatives become weighted averages over the -//! transformed configurations, with a cache that Metropolis moves keep -//! up to date. +//! π along each direction), on a square lattice its whole space group +//! ([`Symmetry::square`]: the translations with the eight rotations and +//! reflections, in one of their four one-dimensional representations, +//! [`PointIrrep`]), and, optionally, the flip of every spin (even or odd). +//! [`Symmetry::check`] proves the elements a group and the characters +//! multiplicative, without which the sum projects onto nothing in +//! particular. Ratios and log-derivatives become weighted averages over the +//! transformed configurations, with a cache that Metropolis moves keep up +//! to date. //! - **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) -//! for expectations with no Monte Carlo noise. Each step is stochastic +//! for expectations with no Monte Carlo noise. With a guide +//! ([`VmcConfig::guide`]) the chains sample `|ψ|²` with a floor under it +//! and each sample is weighted back, so the configurations near a node, +//! where the local energy is large, are visited. Each step is stochastic //! reconfiguration, `θ ← θ − η (S + λ)⁻¹ F`, solved in parameter space or //! in sample space, whichever is smaller; the two give the same step. When //! Cholesky loses a pivot to rounding, the step comes from the //! eigendecomposition instead, with eigenvalues below `10⁻¹²` of the //! largest dropped. A projected network can take its metric there: //! amplitudes that nearly cancel give log-derivatives of order `10⁸`. -//! - **Referees.** [`SpinModel::exact_ground_energy`] is Lanczos on all `2ⁿ` -//! amplitudes (at most 22 spins). [`SpinModel::opsum`] gives the model to +//! - **Referees.** [`SpinModel::exact_ground_energy`] is Lanczos on the +//! amplitudes of a magnetisation sector, or all `2ⁿ` (at most 22 spins). +//! [`SpinModel::exact_ground_energy_in`] is the same within one symmetry +//! sector: Lanczos on `P H P`, re-projected after every product, so a +//! projected network is held to its own sector's answer. It refuses a +//! group that is not one and a model the group does not preserve +//! ([`SpinModel::commutes_with`]). [`SpinModel::opsum`] gives the model to //! DMRG on its exact MPO, for lattices past exact diagonalisation. //! //! Chains draw from their own seeded generators, and every sum runs in a @@ -60,8 +73,34 @@ //! Both spaces take the same Hermitian reconfiguration step. //! - **The projection.** It equals the explicit sum over the group, ratios //! and log-derivatives included. It carries its sector: a translation -//! multiplies `Ψ` by its character, a global flip by its parity. Through -//! accepted moves, the cache stays equal to a fresh one. +//! multiplies `Ψ` by its character, a global flip by its parity, and on the +//! 4 × 4 square a quarter turn, an axis reflection and a diagonal one by +//! their characters in each of the four representations. Through accepted +//! moves, the cache stays equal to a fresh one. +//! - **The groups.** Every square space group on 3 × 3 and 4 × 4, at both +//! momenta, in every representation, with and without the flip, closes and +//! multiplies its characters. A character that does not multiply, a set +//! that does not close, and a character other than `±1` are each refused. +//! - **Guided sampling.** On an 8-site ring, a projected network's guided +//! estimates agree with enumeration: the energy within five standard +//! errors, the variance within 20 %. Steps repeat bit for bit at any thread +//! count. +//! - **Chains start where the sector has amplitude.** A configuration fixed +//! by an element of character −1 has none, and no move from it has a +//! ratio; a chain is redrawn until its start is not one. With 32 chains at +//! momentum (π, π) in A₁ on the 4 × 4 lattice, where a quarter of the +//! configurations are such, every chain starts with amplitude. +//! - **The start.** On the 8-site ring, a projected network within `±0.01` +//! starts on the ferromagnet (`E = n/4`, variance below `10⁻³`) and is +//! within `10⁻²` of it after 40 steps; from `±0.1` it reaches the ground +//! state to `10⁻³` in 200. +//! - **The sector referee.** On a 10-site ring the sector of momentum π and +//! odd parity holds the ground state to `10⁻⁹`, and the sector of momentum +//! 0 and even parity lies above it. On the 4 × 4 lattice at `J₂ = J₁/2`, +//! the sector even under every rotation and reflection at momentum zero +//! holds the ground state to `10⁻⁹`; the one odd under a quarter turn lies +//! above. Open boundaries under ring translations, a flip out of a nonzero +//! magnetisation, and a `σᶻ` field under the flip are refused. //! - **Rings in their sector.** With every configuration enumerated, `α = 1`, //! reaches each ring's exact energy to rounding: //! @@ -125,6 +164,54 @@ //! The variance cannot catch the trap, since an excited eigenstate has none. //! Only a referee can. //! +//! A projected network has a trap of its own at the start. Parameters near +//! zero make it nearly the projection of a constant, the uniform +//! superposition, which at zero magnetisation is the ferromagnet's member: +//! `E = n/4` on a Heisenberg ring, an eigenstate, so its energy gradient +//! vanishes. On the 8-site ring (every configuration enumerated, 200 steps) +//! a `±0.01` start ends `2.8·10⁻²` above the ground state, `±0.1` ends +//! `2.6·10⁻⁴` above, and `±0.3` `5·10⁻⁷`. How long a small start stays depends +//! on the lattice and the group: on the 4 × 4 lattice the translations +//! escaped it within 60 steps and ended the same from either start, while +//! the space group from `±0.01` ended ten times further off than from `±0.1` +//! (the table under the hardest point below). Projected networks in the CLI +//! and `nqs_ground` now start at `±0.1`. +//! +//! A run's sampled variance is not its variance. The space-group network +//! started at `±0.01` reported `1.4·10⁻²` from its final 8000 samples; over +//! every configuration its variance is `0.55`. 89 % of that sits in two +//! configurations related by symmetry, each of probability `2.6·10⁻⁶` +//! (expected 0.02 times in 8000 samples), where the local energy is 310 +//! away from the mean: amplitudes near a node of the network, where the +//! local energy is a large number over a small one. +//! +//! - It is not the chains. Thinning them tenfold reported `1.9·10⁻²`; +//! 64 000 samples happened on one of the two and reported `1.55`. +//! - The network started at `±0.1` is the same: 97 % of its variance +//! (`0.075`) sits in configurations whose probabilities total +//! `2.4·10⁻³`. +//! - The energy holds better, but leans. Across both networks and the three +//! samplings, every sampled energy came out below the enumerated one, by +//! at most `4·10⁻⁴` relative: the tail that goes unsampled lies above the +//! mean. A sampled energy is not a bound. +//! +//! A guide brings the tail into view. Sampling `|ψ|²` with a floor under it +//! ([`VmcConfig::guide`]) and weighting each sample back, the `±0.1` +//! network's sampled variance came within 5 % of the enumerated one at +//! every floor from `0.01` to `1`, four seeds each; without the guide its +//! four seeds gave from `0.11` to `2.0` times it. The `±0.01` network's +//! variance sits in two configurations of 12 870, which 8000 samples cannot +//! visit reliably whatever they are drawn from: guided seeds still gave from +//! `0.1` to `2.3` times it. +//! +//! Training sees what its samples see. Unguided, the `±0.1` space-group +//! network stopped near `10⁻³`, with nearly all its remaining variance in +//! configurations whose probabilities total `2.4·10⁻³`. Guided at `0.1`, +//! the same network, start and steps reach the ground energy to +//! `1.7·10⁻¹²`. Projected networks in the CLI and `nqs_ground` now sample +//! with a guide of `0.1`. Unprojected ones gain nothing from it and sample +//! `|ψ|²`. +//! //! 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 @@ -134,23 +221,49 @@ //! //! # Honest boundaries //! -//! - **The hardest point, within a percent.** At `J₂ = 0.5` on the periodic -//! 4 × 4 lattice, exact diagonalisation gives `−8.45792335`. +//! - **The hardest point, on a small sector.** At `J₂ = 0.5` on the +//! periodic 4 × 4 lattice, exact diagonalisation gives `−8.45792335`, and +//! the sector referee places it in the sector even under every rotation +//! and reflection at momentum zero. The next state in that sector is 13 % +//! higher. //! - Unprojected at `α = 1`, 600 steps of 1000 samples leave every network -//! about 3 % above it: complex with Marshall's sign, complex without, -//! and real with it. -//! - Projected onto momentum zero and even parity, a complex `α = 1` -//! network is 1.7 % above after 600 steps. After 2000 it is -//! `7.7·10⁻³ ± 1.2·10⁻³` above (final estimate of 8000 samples). -//! - Rotations and reflections are left out of the group, and the network -//! is small; what remains is theirs to close. -//! - **Translations and the global flip only.** Rotations and reflections -//! of the lattice are not in the group, and characters are `±1`, so the -//! momenta are 0 and π. Projection multiplies the cost of every amplitude -//! by the group's size. -//! - **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. +//! about 3 % above it (sampled estimates): complex with Marshall's sign, +//! complex without, and real with it. A guide does not change that +//! (complex, `2.78·10⁻²` over every configuration). +//! - Projected, a complex `α = 1` network after 600 steps of 1000 samples, +//! its energy then taken over every configuration (`nqs_ground +//! --exact-final`): +//! +//! | projection, even parity | start | guide | relative error | +//! |----------------------------|---------|-------|----------------| +//! | translations | `±0.01` | none | `1.87·10⁻²` | +//! | translations | `±0.1` | none | `1.90·10⁻²` | +//! | space group, A₁ | `±0.01` | none | `2.09·10⁻²` | +//! | space group, A₁ | `±0.1` | none | `1.27·10⁻³` | +//! | space group, A₁ | `±0.3` | none | `1.09·10⁻³` | +//! | space group, A₁ | `±0.1` | `0.1` | `1.7·10⁻¹²` | +//! | space group, A₁ | `±0.1` | `1` | `1.2·10⁻¹²` | +//! +//! Unguided from `±0.1`, 2000 steps reach `4.60·10⁻⁴`, and `α = 2` +//! reaches `6.91·10⁻⁴` in 600. Guided, the variance over every +//! configuration is below `4·10⁻¹¹`. +//! - The sector is small: 107 configurations up to its symmetries, against +//! the network's 288 parameters. Reaching its ground state to rounding +//! says the training finds what the network can hold. It says nothing of +//! a sector larger than the network. +//! - **One-dimensional representations only.** Characters are `±1`: momenta +//! 0 and π, and the four one-dimensional representations of the square's +//! rotations and reflections. The two-dimensional one and every other +//! momentum are out of reach. Projection multiplies the cost of every +//! amplitude by the group's size, 256 on the 4 × 4 lattice with the flip; +//! 600 steps took about seven times as long as with the translations. +//! - **Error bars ignore autocorrelation and the tail.** The reported error +//! is the standard error of uncorrelated samples, from the sampled +//! variance. Chains thinned by two sweeps are correlated, and the sampled +//! variance can miss nearly all of the variance (see Measured), so the true +//! error is larger. A guide recovers a tail spread over many +//! configurations, not one held by a few. Only enumeration +//! (`--exact-final`) settles it, up to 24 spins. //! - **Dense reconfiguration.** `S` (or its sample-space twin) is formed and //! factored in full each step. That bounds the network and sample sizes //! long before iterative solvers would. @@ -300,14 +413,27 @@ impl SpinModel { s } - /// `out = H inp` on all `2ⁿ` amplitudes, site 0 the most significant bit - /// (a set bit is `|↓⟩`). - fn apply_dense(&self, inp: &[f64], out: &mut [f64]) { + /// The states with total `Σ σᶻ = magnetisation` (all `2ⁿ` when none), in + /// ascending order, and each one's place among them. + fn sector_states(&self, magnetisation: Option) -> (Vec, Vec) { + let n = self.n; + let states: Vec = (0..1usize << n).filter(|&x| magnetisation.is_none_or(|m| n as i32 - 2 * x.count_ones() as i32 == m)).collect(); + let mut place = vec![u32::MAX; 1usize << n]; + for (k, &x) in states.iter().enumerate() { + place[x] = k as u32; + } + (states, place) + } + + /// `out = H inp` on the amplitudes of `states` alone, which `H` must map + /// among themselves. Site 0 is the most significant bit of a state, and + /// a set bit is `|↓⟩`. + fn apply_sector(&self, states: &[usize], place: &[u32], inp: &[f64], out: &mut [f64]) { let n = self.n; let bit = |i: usize| 1usize << (n - 1 - i); let spin = |x: usize, i: usize| if x & bit(i) == 0 { 1.0 } else { -1.0 }; out.iter_mut().for_each(|o| *o = 0.0); - for (x, &) in inp.iter().enumerate() { + for (k, (&x, &)) in states.iter().zip(inp).enumerate() { if amp == 0.0 { continue; } @@ -318,21 +444,20 @@ impl SpinModel { for &(i, c) in &self.z { diag += c * spin(x, i); } - out[x] += diag * amp; + out[k] += diag * amp; for &(i, j, c) in &self.exchange { if spin(x, i) != spin(x, j) { - out[x ^ bit(i) ^ bit(j)] += 2.0 * c * amp; + out[place[x ^ bit(i) ^ bit(j)] as usize] += 2.0 * c * amp; } } for &(i, c) in &self.x { - out[x ^ bit(i)] += c * amp; + out[place[x ^ bit(i)] as usize] += c * amp; } } } - /// The exact ground energy by Lanczos on all `2ⁿ` amplitudes (at most 22 - /// sites), restricted to total `Σ σᶻ = magnetisation` when given (the - /// model must conserve it). + /// The exact ground energy by Lanczos (at most 22 sites), restricted to + /// total `Σ σᶻ = magnetisation` when given (the model must conserve it). pub fn exact_ground_energy(&self, magnetisation: Option) -> Option { let n = self.n; if n == 0 || n > 22 { @@ -341,13 +466,95 @@ impl SpinModel { if magnetisation.is_some() && !self.conserves_magnetisation() { return None; } - let dim = 1usize << n; - let in_sector = |x: usize| magnetisation.is_none_or(|m| n as i32 - 2 * x.count_ones() as i32 == m); - let start: Vec = (0..dim as u64).map(|x| if in_sector(x as usize) { 1.0 + (x.wrapping_mul(2_654_435_761) % 1000) as f64 / 1000.0 } else { 0.0 }).collect(); - if start.iter().all(|&v| v == 0.0) { + let (states, place) = self.sector_states(magnetisation); + if states.is_empty() { return None; } - let apply = |inp: &[f64], out: &mut [f64]| self.apply_dense(inp, out); + let start: Vec = states.iter().map(|&x| 1.0 + ((x as u64).wrapping_mul(2_654_435_761) % 1000) as f64 / 1000.0).collect(); + let apply = |inp: &[f64], out: &mut [f64]| self.apply_sector(&states, &place, inp, out); + Some(lanczos(&apply, &start, 60, 30, 1e-9, false).0) + } + + /// Whether every element of `sym` maps the model onto itself: each bond + /// and field to one of the same strength, and under a flip of every spin + /// no `σᶻ` field, which the flip reverses. + pub fn commutes_with(&self, sym: &Symmetry) -> bool { + let n = self.n; + if sym.perms.iter().any(|p| p.len() != n) { + return false; + } + let bonds = |terms: &[(usize, usize, f64)], p: &[usize]| { + let mut b: Vec<(usize, usize, u64)> = terms.iter().map(|&(i, j, c)| (p[i].min(p[j]), p[i].max(p[j]), c.to_bits())).collect(); + b.sort_unstable(); + b + }; + let fields = |terms: &[(usize, f64)], p: &[usize]| { + let mut f: Vec<(usize, u64)> = terms.iter().map(|&(i, c)| (p[i], c.to_bits())).collect(); + f.sort_unstable(); + f + }; + let identity: Vec = (0..n).collect(); + let (exchange, zz, x, z) = (bonds(&self.exchange, &identity), bonds(&self.zz, &identity), fields(&self.x, &identity), fields(&self.z, &identity)); + sym.perms.iter().zip(&sym.flips).all(|(p, &flip)| bonds(&self.exchange, p) == exchange && bonds(&self.zz, p) == zz && fields(&self.x, p) == x && fields(&self.z, p) == z && !(flip && !self.z.is_empty())) + } + + /// The exact ground energy within the sector `sym` projects onto (at most + /// 22 sites), and in the magnetisation sector when given: Lanczos on + /// `P H P` with `P = (1/|G|) Σ_g χ(g) g`, from a start in the sector and + /// projected again after every product, so rounding cannot carry it into + /// a lower sector. `None` when `sym` is not a group with characters, when + /// the model is not invariant under it, when a flip would leave the + /// magnetisation sector, or when the sector is empty. + pub fn exact_ground_energy_in(&self, magnetisation: Option, sym: &Symmetry) -> Option { + let n = self.n; + if n == 0 || n > 22 || sym.check().is_err() || !self.commutes_with(sym) { + return None; + } + if magnetisation.is_some() && !self.conserves_magnetisation() { + return None; + } + if magnetisation.is_some_and(|m| m != 0) && sym.flips.iter().any(|&f| f) { + return None; + } + let bit = |i: usize| 1usize << (n - 1 - i); + let (states, place) = self.sector_states(magnetisation); + // images[k·|G| + g]: the place of g·states[k], whose spin at i is the + // state's at perms[g][i], flipped if the element flips. + let order = sym.len(); + let mut images = Vec::with_capacity(states.len() * order); + for &x in &states { + for (p, &flip) in sym.perms.iter().zip(&sym.flips) { + let mut y = 0usize; + for (i, &pi) in p.iter().enumerate() { + if (x & bit(pi) != 0) != flip { + y |= bit(i); + } + } + images.push(place[y]); + } + } + let project = |v: &[f64], out: &mut [f64]| { + for (k, o) in out.iter_mut().enumerate() { + let row = &images[k * order..(k + 1) * order]; + let mut acc = 0.0; + for (&y, &c) in row.iter().zip(&sym.chars) { + acc += c * v[y as usize]; + } + *o = acc / order as f64; + } + }; + let raw: Vec = states.iter().map(|&x| 1.0 + ((x as u64).wrapping_mul(2_654_435_761) % 1000) as f64 / 1000.0).collect(); + let mut start = vec![0.0; states.len()]; + project(&raw, &mut start); + let (before, after) = (raw.iter().map(|v| v * v).sum::(), start.iter().map(|v| v * v).sum::()); + if after <= 1e-20 * before { + return None; + } + let apply = |inp: &[f64], out: &mut [f64]| { + let mut h = vec![0.0; states.len()]; + self.apply_sector(&states, &place, inp, &mut h); + project(&h, out); + }; Some(lanczos(&apply, &start, 60, 30, 1e-9, false).0) } } @@ -438,10 +645,43 @@ impl Amplitude for C { } } +/// One of the four one-dimensional representations of the square's eight +/// symmetries, by its character on a rotation and on a reflection. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum PointIrrep { + /// Even under every rotation and reflection. + A1, + /// Even under rotations, odd under every reflection. + A2, + /// Odd under a quarter turn and a diagonal reflection, even under a + /// reflection in an axis. + B1, + /// Odd under a quarter turn and a reflection in an axis, even under a + /// diagonal reflection. + B2, +} + +impl PointIrrep { + /// The character of `turns` quarter turns after a reflection in the + /// vertical axis if `reflect`. An odd number of turns after it makes a + /// diagonal reflection, an even number an axis one. + pub fn character(self, turns: usize, reflect: bool) -> f64 { + let rotation = if turns % 2 == 1 { -1.0 } else { 1.0 }; + let reflection = if reflect { -1.0 } else { 1.0 }; + match self { + PointIrrep::A1 => 1.0, + PointIrrep::A2 => reflection, + PointIrrep::B1 => rotation, + PointIrrep::B2 => rotation * reflection, + } + } +} + /// A symmetry group of the lattice: site permutations, each possibly /// followed by flipping every spin, each with a real character. The /// projection `Ψ(s) = Σ_g χ(g) ψ(g·s)` picks the sector: a momentum of 0 or -/// π along each direction and, with the flip, even or odd parity. +/// π along each direction, on a square lattice a representation of its +/// rotations and reflections, and with the flip even or odd parity. #[derive(Clone, Debug, PartialEq)] pub struct Symmetry { /// `perms[g][i]`: the site whose spin lands on site `i`. @@ -491,16 +731,20 @@ impl Symmetry { Symmetry::with_flips(perms, flips, chars) } - /// The translations of a ring of `n` sites, at momentum 0 or π. + /// The translations of a ring of `n` sites, at momentum 0 or π (π needs + /// `n` even: going round the ring must multiply by one). pub fn ring(n: usize, momentum_pi: bool) -> Symmetry { + assert!(!momentum_pi || n.is_multiple_of(2), "momentum π on a ring of odd length"); let perms = (0..n).map(|t| (0..n).map(|i| (i + t) % n).collect()).collect(); let chars = (0..n).map(|t| if momentum_pi && t % 2 == 1 { -1.0 } else { 1.0 }).collect(); Symmetry::new(perms, chars) } /// The translations of an `lx × ly` torus (site `y·lx + x`), at momentum - /// 0 or π along each direction. + /// 0 or π along each direction (π along one needs its length even). pub fn torus(lx: usize, ly: usize, kx_pi: bool, ky_pi: bool) -> Symmetry { + assert!(!kx_pi || lx.is_multiple_of(2), "momentum π along an odd length"); + assert!(!ky_pi || ly.is_multiple_of(2), "momentum π along an odd length"); let mut perms = Vec::new(); let mut chars = Vec::new(); for ty in 0..ly { @@ -513,6 +757,91 @@ impl Symmetry { Symmetry::new(perms, chars) } + /// The whole space group of an `l × l` torus: every translation, each + /// composed with the eight symmetries of the square about a site (quarter + /// turns, and reflections in the axes and the diagonals). Momentum 0 or + /// `(π, π)`, both of which every rotation and reflection leaves alone, so + /// a character of the eight extends to the group; `irrep` picks which. + pub fn square(l: usize, momentum_pi: bool, irrep: PointIrrep) -> Symmetry { + // At l = 2 a reflection and the identity move every site alike, and + // a character that tells them apart cannot be one. + assert!(l >= 3, "the square's symmetries are distinct from l = 3"); + assert!(!momentum_pi || l.is_multiple_of(2), "momentum π along an odd length"); + let mut perms = Vec::new(); + let mut chars = Vec::new(); + for ty in 0..l { + for tx in 0..l { + for turns in 0..4 { + for reflect in [false, true] { + perms.push( + (0..l * l) + .map(|i| { + let (mut x, mut y) = (i % l, i / l); + if reflect { + x = (l - x) % l; + } + for _ in 0..turns { + (x, y) = ((l - y) % l, x); + } + ((y + ty) % l) * l + (x + tx) % l + }) + .collect(), + ); + let momentum = if momentum_pi && (tx + ty) % 2 == 1 { -1.0 } else { 1.0 }; + chars.push(momentum * irrep.character(turns, reflect)); + } + } + } + } + Symmetry::new(perms, chars) + } + + /// Whether the elements form a group under composition and every + /// character is `±1` and multiplies as the elements do: what makes + /// `Σ_g χ(g) ψ(g·s)` a projection onto one sector, not onto a mixture or + /// onto nothing. + pub fn check(&self) -> Result<(), String> { + if self.is_empty() { + return Err("no elements".into()); + } + if let Some(c) = self.chars.iter().find(|&&c| c != 1.0 && c != -1.0) { + return Err(format!("character {c} is not ±1")); + } + let mut index: std::collections::BTreeMap<(&[usize], bool), usize> = std::collections::BTreeMap::new(); + for g in 0..self.len() { + match index.get(&(self.perms[g].as_slice(), self.flips[g])) { + Some(&h) if self.chars[h] != self.chars[g] => return Err(format!("elements {h} and {g} act alike with different characters")), + Some(_) => {} + None => { + index.insert((self.perms[g].as_slice(), self.flips[g]), g); + } + } + } + let n = self.perms[0].len(); + let identity: Vec = (0..n).collect(); + match index.get(&(identity.as_slice(), false)) { + Some(&e) if self.chars[e] == 1.0 => {} + Some(_) => return Err("the identity's character is not 1".into()), + None => return Err("no identity".into()), + } + // (g·(h·s))ᵢ = s at p_h(p_g(i)), flipped by either. + let mut q = vec![0usize; n]; + for g in 0..self.len() { + for h in 0..self.len() { + for (i, qi) in q.iter_mut().enumerate() { + *qi = self.perms[h][self.perms[g][i]]; + } + let Some(&k) = index.get(&(q.as_slice(), self.flips[g] ^ self.flips[h])) else { + return Err(format!("elements {g} and {h} compose to one outside the set")); + }; + if self.chars[k] != self.chars[g] * self.chars[h] { + return Err(format!("the characters of {g} and {h} do not multiply to their product's")); + } + } + } + Ok(()) + } + pub fn len(&self) -> usize { self.perms.len() } @@ -874,11 +1203,18 @@ pub struct VmcConfig { /// Threads over chains and configurations; the result does not depend on /// it. pub threads: usize, + /// Sample `|ψ|² + guide·|ψ̄|²` instead of `|ψ|²`, `|ψ̄|` the chains' + /// typical amplitude when the batch starts, and weight each sample by + /// `|ψ|²` over that. Near a node of the network the local energy is a + /// large number over a small one, and `|ψ|²` visits there too rarely for + /// the variance to be measured; the floor visits, the weight corrects. + /// Zero samples `|ψ|²`. + pub guide: f64, } impl Default for VmcConfig { fn default() -> VmcConfig { - VmcConfig { samples: 2000, chains: 8, burn: 50, thin: 2, exact: false, magnetisation: None, lr: 0.05, shift: 1e-3, seed: 1, threads: default_threads() } + VmcConfig { samples: 2000, chains: 8, burn: 50, thin: 2, exact: false, magnetisation: None, lr: 0.05, shift: 1e-3, seed: 1, threads: default_threads(), guide: 0.0 } } } @@ -892,7 +1228,7 @@ pub struct Estimate { /// The statistical error of `energy`, from the variance and the sample /// count, ignoring autocorrelation (zero when exact). pub error: f64, - /// Fraction of Metropolis moves accepted (one when exact). + /// Fraction of proposed Metropolis moves accepted (one when exact). pub acceptance: f64, } @@ -913,6 +1249,30 @@ struct Batch { energies: Vec, derivs: Vec, acceptance: f64, + /// Samples drawn from a guide, weighted back to `|ψ|²`. + guided: bool, +} + +/// The guide `|ψ|² + floor·e^{2c}`, by its logarithm. +#[derive(Clone, Copy, Debug)] +struct Guide { + /// The chains' mean `ln |ψ|` when the batch starts. + c: f64, + ln_floor: f64, +} + +impl Guide { + /// `ln (e^{2(l−c)} + floor)`. + fn ln_q(self, l: f64) -> f64 { + let x = 2.0 * (l - self.c); + let (hi, lo) = if x > self.ln_floor { (x, self.ln_floor) } else { (self.ln_floor, x) }; + hi + ln(1.0 + exp(lo - hi)) + } + + /// A sample's weight `|ψ|² / (|ψ|² + floor·e^{2c})`. + fn weight(self, l: f64) -> f64 { + exp(2.0 * (l - self.c) - self.ln_q(l)) + } } impl Vmc { @@ -925,7 +1285,16 @@ impl Vmc { let chains = (0..cfg.chains.max(1)) .map(|c| { let mut rng = Rng(cfg.seed ^ (c as u64).wrapping_mul(0xa076_1d64_78bd_642f)); - let s = start_configuration(n, cfg.magnetisation, &mut rng); + let mut s = start_configuration(n, cfg.magnetisation, &mut rng); + // A configuration the projection forbids has no amplitude, so + // no move from it has a ratio, and a chain started there would + // never leave. Draw again; an unprojected network never does. + let mut draws = 1; + while !rbm.log_abs_cached(&s, &rbm.cache(&s)).is_finite() { + assert!(draws < 10_000, "the projection vanishes on every start drawn: an empty sector?"); + s = start_configuration(n, cfg.magnetisation, &mut rng); + draws += 1; + } (s, rng) }) .collect(); @@ -933,7 +1302,7 @@ impl Vmc { if !cfg.exact { for c in 0..v.chains.len() { let (mut s, mut rng) = v.chains[c].clone(); - metropolis(&v.rbm, &mut s, &mut rng, cfg.burn * n, cfg.magnetisation.is_some()); + metropolis(&v.rbm, &mut s, &mut rng, cfg.burn * n, cfg.magnetisation.is_some(), None); v.chains[c] = (s, rng); } } @@ -999,7 +1368,7 @@ impl Vmc { let top = logs.iter().copied().fold(f64::NEG_INFINITY, f64::max); let raw: Vec = logs.iter().map(|l| exp(2.0 * (l - top))).collect(); let z: f64 = raw.iter().sum(); - Batch { weights: raw.iter().map(|r| r / z).collect(), energies, derivs, acceptance: 1.0 } + Batch { weights: raw.iter().map(|r| r / z).collect(), energies, derivs, acceptance: 1.0, guided: false } } fn sample(&mut self) -> Batch { @@ -1011,21 +1380,31 @@ 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 guide = (self.cfg.guide > 0.0).then(|| { + let logs: Vec = self.chains.iter().map(|(s, _)| rbm.log_abs_cached(s, &rbm.cache(s))).filter(|l| l.is_finite()).collect(); + Guide { c: logs.iter().sum::() / logs.len().max(1) as f64, ln_floor: ln(self.cfg.guide) } + }); + let run = |chain: &mut (Vec, Rng)| -> (Vec, Vec, Vec, f64) { let (s, rng) = chain; let mut energies = Vec::with_capacity(per_chain); let mut derivs = vec![T::ZERO; per_chain * p]; - let mut accepted = 0usize; + let mut weights = Vec::with_capacity(if guide.is_some() { per_chain } else { 0 }); + let (mut accepted, mut proposed) = (0usize, 0usize); for k in 0..per_chain { - accepted += metropolis(rbm, s, rng, thin, exchange); + let (a, m) = metropolis(rbm, s, rng, thin, exchange, guide); + accepted += a; + proposed += m; let cache = rbm.cache(s); + if let Some(g) = guide { + weights.push(g.weight(rbm.log_abs_cached(s, &cache))); + } energies.push(local_energy(model, rbm, s, &cache)); rbm.derivatives_cached(s, &cache, &mut derivs[k * p..(k + 1) * p]); } - (energies, derivs, accepted as f64 / (per_chain * thin) as f64) + (energies, derivs, weights, accepted as f64 / proposed.max(1) 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, Vec, f64)> = if threads <= 1 || cfg!(target_arch = "wasm32") { self.chains.iter_mut().map(run).collect() } else { let per = nc.div_ceil(threads); @@ -1044,13 +1423,21 @@ impl Vmc { let total = per_chain * nc; let mut energies = Vec::with_capacity(total); let mut derivs = Vec::with_capacity(total * p); + let mut raw = Vec::with_capacity(if guide.is_some() { total } else { 0 }); let mut acc = 0.0; - for (e, d, a) in results { + for (e, d, w, a) in results { energies.extend(e); derivs.extend(d); + raw.extend(w); acc += a; } - Batch { weights: vec![1.0 / total as f64; total], energies, derivs, acceptance: acc / nc as f64 } + let weights = if guide.is_some() { + let z: f64 = raw.iter().sum(); + raw.iter().map(|w| w / z).collect() + } else { + vec![1.0 / total as f64; total] + }; + Batch { weights, energies, derivs, acceptance: acc / nc as f64, guided: guide.is_some() } } /// Estimate the energy without changing the machine. @@ -1085,7 +1472,14 @@ fn summarise(b: &Batch, exact: bool) -> Estimate { d.scale(w).mul(d.conj()).re() }) .sum(); - let error = if exact { 0.0 } else { (variance / b.energies.len() as f64).sqrt() }; + // Weighted samples: the variance of a weighted mean, Σ w² |E − Ē|². + let error = if exact { + 0.0 + } else if b.guided { + b.weights.iter().zip(&b.energies).map(|(&w, e)| w * w * e.sub(et).norm2()).sum::().sqrt() + } else { + (variance / b.energies.len() as f64).sqrt() + }; Estimate { energy, variance, error, acceptance: b.acceptance } } @@ -1260,12 +1654,14 @@ 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 { +/// spins), accepted with probability `|ψ(s′)/ψ(s)|²`, or with a guide its +/// `q(s′)/q(s)`. Returns the number accepted and the number proposed: a draw +/// of two parallel spins proposes nothing. +fn metropolis(rbm: &Rbm, s: &mut [i8], rng: &mut Rng, moves: usize, exchange: bool, guide: Option) -> (usize, usize) { let n = s.len(); let mut cache = rbm.cache(s); - let mut accepted = 0; + let (mut accepted, mut proposed) = (0, 0); + let mut l = if guide.is_some() { rbm.log_abs_cached(s, &cache) } else { 0.0 }; for _ in 0..moves { let flips: Vec = if exchange { let i = rng.below(n); @@ -1277,13 +1673,25 @@ fn metropolis(rbm: &Rbm, s: &mut [i8], rng: &mut Rng, moves: us } else { vec![rng.below(n)] }; + proposed += 1; let r = rbm.ratio_cached(s, &cache, &flips); - if rng.uniform() < r.norm2() { + let r2 = r.norm2(); + let p = match guide { + None => r2, + // A ratio below 10⁻¹⁰ in size is a move to where the projection + // vanishes, to rounding: outside the guide, as it is outside |ψ|². + Some(_) if r2 < 1e-20 => 0.0, + Some(g) => exp(g.ln_q(l + 0.5 * ln(r2)) - g.ln_q(l)), + }; + if rng.uniform() < p { rbm.accept(s, &mut cache, &flips, r); accepted += 1; + if guide.is_some() { + l += 0.5 * ln(r2); + } } } - accepted + (accepted, proposed) } #[cfg(test)] @@ -1313,12 +1721,13 @@ mod tests { m.zz.push((0, 8, 0.2)); let dense = m.opsum().mpo().to_dense().unwrap(); let dim = 1usize << n; + let (states, place) = m.sector_states(None); let mut col = vec![0.0; dim]; let mut out = vec![0.0; dim]; for c in 0..dim { col.iter_mut().for_each(|x| *x = 0.0); col[c] = 1.0; - m.apply_dense(&col, &mut out); + m.apply_sector(&states, &place, &col, &mut out); for r in 0..dim { assert!((out[r] - dense[r * dim + c]).abs() < 1e-12, "({r}, {c})"); } @@ -1632,4 +2041,173 @@ mod tests { assert!(e.energy >= exact - 1e-9, "variational: {} vs {exact}", e.energy); assert!((e.energy - exact) / exact.abs() < 1e-8, "{} vs {exact}", e.energy); } + + #[test] + fn square_groups_are_groups_with_characters() { + for l in [3, 4] { + for momentum_pi in [false, true] { + if momentum_pi && l % 2 == 1 { + continue; + } + for irrep in [PointIrrep::A1, PointIrrep::A2, PointIrrep::B1, PointIrrep::B2] { + let sym = Symmetry::square(l, momentum_pi, irrep); + assert_eq!(sym.len(), 8 * l * l); + assert_eq!(sym.check(), Ok(()), "{l} {momentum_pi} {irrep:?}"); + } + for odd in [false, true] { + assert_eq!(Symmetry::square(l, momentum_pi, PointIrrep::B2).with_spin_flip(odd).check(), Ok(())); + } + } + } + // Characters that do not multiply, and elements that do not close. + let ring = Symmetry::ring(4, false); + assert!(Symmetry::new(ring.perms.clone(), vec![1.0, -1.0, 1.0, 1.0]).check().is_err()); + assert!(Symmetry::new(ring.perms[..2].to_vec(), vec![1.0, 1.0]).check().is_err()); + assert!(Symmetry::new(ring.perms.clone(), vec![1.0, 0.5, 1.0, 0.5]).check().is_err()); + } + + #[test] + fn a_square_projection_carries_its_representation() { + let (l, n) = (4, 16); + let quarter_turn = |x: usize, y: usize| ((l - y) % l, x); + let axis = |x: usize, y: usize| ((l - x) % l, y); + let diagonal = |x: usize, y: usize| (y, x); + let step = |x: usize, y: usize| ((x + 1) % l, y); + for irrep in [PointIrrep::A1, PointIrrep::A2, PointIrrep::B1, PointIrrep::B2] { + let rbm = Rbm::complex(n, 1, 0.3, 11).with_symmetry(Symmetry::square(l, true, irrep).with_spin_flip(false)); + let mut rng = Rng(5); + let mut seen = 0; + // A configuration some element with character −1 leaves alone + // has no amplitude in the sector; on sixteen sites that is up to + // two in five of them. + for _ in 0..40 { + if seen == 6 { + break; + } + let s = start_configuration(n, Some(0), &mut rng); + let base = projected(&rbm, &s); + if base.norm2() < 1e-16 { + continue; + } + seen += 1; + let moved = |f: &dyn Fn(usize, usize) -> (usize, usize)| -> Vec { + (0..n) + .map(|i| { + let (x, y) = f(i % l, i / l); + s[y * l + x] + }) + .collect() + }; + // Momentum (π, π): one step along an axis multiplies by −1. + let want = [(irrep.character(1, false), moved(&quarter_turn)), (irrep.character(0, true), moved(&axis)), (irrep.character(1, true), moved(&diagonal)), (-1.0, moved(&step))]; + for (chi, t) in want { + let got = projected(&rbm, &t); + assert!(got.sub(base.scale(chi)).norm2().sqrt() < 1e-10 * base.norm2().sqrt(), "{irrep:?}: {got:?} vs {chi} × {base:?}"); + } + } + assert_eq!(seen, 6, "{irrep:?}: too few configurations survive the projection"); + } + } + + #[test] + fn a_projected_start_near_zero_sits_on_the_ferromagnet() { + // An 8-site ring projected onto momentum 0 and even parity, the + // ground state's sector. Parameters within ±0.01 make the network + // nearly the projection of a constant: the uniform superposition, + // which is the ferromagnet's member at zero magnetisation, energy + // n/4, an eigenstate. Reconfiguration barely moves it. + let n = 8; + let model = SpinModel::heisenberg(n, &chain_edges(n, true), 1.0); + let exact = model.exact_ground_energy(Some(0)).unwrap(); + let ferromagnet = n as f64 / 4.0; + let cfg = VmcConfig { exact: true, magnetisation: Some(0), lr: 0.02, shift: 1e-3, ..VmcConfig::default() }; + let run = |init: f64, steps: usize| { + let rbm = Rbm::complex(n, 1, init, 7).with_symmetry(Symmetry::ring(n, false).with_spin_flip(false)); + let mut v = Vmc::new(model.clone(), rbm, cfg); + let first = v.estimate(); + for _ in 0..steps { + v.step(); + } + (first, v.estimate()) + }; + let (first, after) = run(0.01, 40); + assert!((first.energy - ferromagnet).abs() < 1e-3 && first.variance < 1e-3, "{first:?}"); + assert!((after.energy - ferromagnet).abs() < 1e-2, "{after:?}"); + let (_, after) = run(0.1, 200); + assert!((after.energy - exact) / exact.abs() < 1e-3, "{} vs {exact}", after.energy); + } + + #[test] + fn guided_sampling_agrees_with_enumeration() { + // A projected network on an 8-site ring, momentum 0 and even parity. + let n = 8; + let model = SpinModel::heisenberg(n, &chain_edges(n, true), 1.0); + let rbm = Rbm::complex(n, 1, 0.3, 5).with_symmetry(Symmetry::ring(n, false).with_spin_flip(false)); + let exact = Vmc::new(model.clone(), rbm.clone(), VmcConfig { exact: true, magnetisation: Some(0), ..VmcConfig::default() }).estimate(); + let cfg = VmcConfig { samples: 4000, chains: 8, magnetisation: Some(0), guide: 0.5, threads: 1, ..VmcConfig::default() }; + let g = Vmc::new(model.clone(), rbm.clone(), cfg).estimate(); + assert!((g.energy - exact.energy).abs() < 5.0 * g.error, "{g:?} vs {exact:?}"); + assert!((g.variance - exact.variance).abs() < 0.2 * exact.variance, "{g:?} vs {exact:?}"); + // The same bits at any thread count. + let cfg = |threads| VmcConfig { samples: 200, magnetisation: Some(0), guide: 0.5, threads, ..VmcConfig::default() }; + let mut one = Vmc::new(model.clone(), rbm.clone(), cfg(1)); + let mut many = Vmc::new(model, rbm, cfg(3)); + for _ in 0..2 { + assert_eq!(one.step(), many.step()); + } + assert_eq!(one.rbm, many.rbm); + } + + #[test] + fn no_chain_starts_where_the_projection_vanishes() { + // At momentum (π, π), even under every rotation and reflection, about + // a quarter of the 4 × 4 lattice's configurations are fixed by an + // element of character −1 and have no amplitude. A chain there has + // no move with a ratio and would never leave. + let (l, n) = (4, 16); + let model = SpinModel::heisenberg(n, &square_edges(l, l, true), 1.0); + let rbm = Rbm::complex(n, 1, 0.1, 3).with_symmetry(Symmetry::square(l, true, PointIrrep::A1).with_spin_flip(false)); + let cfg = VmcConfig { chains: 32, burn: 1, magnetisation: Some(0), threads: 1, ..VmcConfig::default() }; + let v = Vmc::new(model, rbm, cfg); + for (c, (s, _)) in v.chains.iter().enumerate() { + assert!(v.rbm.log_abs_cached(s, &v.rbm.cache(s)).is_finite(), "chain {c} sits where the projection vanishes: {s:?}"); + } + } + + #[test] + fn sector_diagonalisation_splits_the_spectrum() { + // A ring of ten: half its length is odd, so the singlet ground state + // has momentum π and odd parity under the flip, and the sector of + // momentum 0 and even parity starts higher. + let n = 10; + let model = SpinModel::heisenberg(n, &chain_edges(n, true), 1.0); + let all = model.exact_ground_energy(Some(0)).unwrap(); + let ground = model.exact_ground_energy_in(Some(0), &Symmetry::ring(n, true).with_spin_flip(true)).unwrap(); + let other = model.exact_ground_energy_in(Some(0), &Symmetry::ring(n, false).with_spin_flip(false)).unwrap(); + assert!((ground - all).abs() < 1e-9, "{ground} vs {all}"); + assert!(other > all + 0.1, "{other} vs {all}"); + // Refusals: a model the group does not preserve, a flip out of the + // magnetisation, and a field the flip reverses. + let open = SpinModel::heisenberg(n, &chain_edges(n, false), 1.0); + assert_eq!(open.exact_ground_energy_in(Some(0), &Symmetry::ring(n, false)), None); + assert_eq!(model.exact_ground_energy_in(Some(2), &Symmetry::ring(n, false).with_spin_flip(false)), None); + let mut field = model.clone(); + field.z = (0..n).map(|i| (i, 0.1)).collect(); + assert!(field.commutes_with(&Symmetry::ring(n, false))); + assert!(!field.commutes_with(&Symmetry::ring(n, false).with_spin_flip(false))); + } + + #[test] + fn the_frustrated_square_ground_state_is_fully_symmetric() { + // 4 × 4 at J₂ = J₁/2: the ground state is even under every rotation + // and reflection, at momentum zero. + let (l, n) = (4, 16); + let model = SpinModel::j1j2(n, &square_edges(l, l, true), &square_diagonals(l, l, true), 1.0, 0.5); + let all = model.exact_ground_energy(Some(0)).unwrap(); + let sector = |momentum_pi, irrep| model.exact_ground_energy_in(Some(0), &Symmetry::square(l, momentum_pi, irrep)).unwrap(); + let a1 = sector(false, PointIrrep::A1); + assert!((a1 - all).abs() < 1e-9, "{a1} vs {all}"); + let b1 = sector(false, PointIrrep::B1); + assert!(b1 > all + 0.1, "{b1} vs {all}"); + } } diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index 43a906a55..5b3c2b035 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.38" +version = "0.3.39" 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 fa6d2fd1b..3239d1cf5 100644 --- a/wai/wai-rs/src/bin/wai_quantum.rs +++ b/wai/wai-rs/src/bin/wai_quantum.rs @@ -174,8 +174,8 @@ fn usage() -> ! { 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|j1j2|ising] [--lx L] [--ly L] [--j2 J] [--h H] [--periodic]\n \ - [--complex] [--project [--kpi] [--odd]] [--alpha A] [--steps N] [--samples S] [--lr R]\n \ - [--exact] [--seed S] [--json]\n\n\ + [--complex] [--project [--kpi] [--point-group a1|a2|b1|b2] [--odd]] [--alpha A] [--steps N] [--samples S] [--lr R]\n \ + [--init W] [--guide ε] [--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 \ @@ -1190,21 +1190,53 @@ fn main() { // 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 mut 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 }; - // Projection onto translations at momentum 0 or π, and a parity - // under the global flip; the lattice must be periodic. + // Projection onto translations at momentum 0 or π, on a square + // lattice optionally its rotations and reflections too, and a + // parity under the global flip; the lattice must be periodic. + let point = opt(args, "--point-group", "").map(|g| match g.as_str() { + "a1" => PointIrrep::A1, + "a2" => PointIrrep::A2, + "b1" => PointIrrep::B1, + "b2" => PointIrrep::B2, + _ => die("--point-group: a1, a2, b1 or b2"), + }); let symmetry = flag(args, "--project").then(|| { if !periodic { die("--project needs --periodic"); } let kpi = flag(args, "--kpi"); - let base = if ly == 1 { Symmetry::ring(n, kpi) } else { Symmetry::torus(lx, ly, kpi, kpi) }; + if kpi && (lx % 2 == 1 || (ly > 1 && ly % 2 == 1)) { + die("--kpi needs even lengths"); + } + let base = match point { + Some(_) if lx != ly || lx < 3 => die("--point-group needs a square lattice, --lx = --ly ≥ 3"), + Some(irrep) => Symmetry::square(lx, kpi, irrep), + None if ly == 1 => Symmetry::ring(n, kpi), + None => Symmetry::torus(lx, ly, kpi, kpi), + }; base.with_spin_flip(flag(args, "--odd")) }); + if point.is_some() && symmetry.is_none() { + die("--point-group needs --project"); + } + // Parameters start uniform in ±init. A projected network that + // starts near zero is nearly the projection of a constant, the + // ferromagnet: an eigenstate, which reconfiguration barely moves. + let init: f64 = num(args, "--init", "", if symmetry.is_some() { 0.1 } else { 0.01 }); + // A projected network's amplitude is a sum that cancels near its + // nodes, where the local energy is large and |ψ|² rarely + // samples; a guide floor visits there. + cfg.guide = num(args, "--guide", "", if symmetry.is_some() { 0.1 } else { 0.0 }); + // The exact ground energy in the projection's own sector. + let in_sector = match &symmetry { + Some(sym) if n <= 20 => model.exact_ground_energy_in(sector, sym), + _ => None, + }; 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); @@ -1215,7 +1247,7 @@ fn main() { (params, last, vmc.estimate()) } let (params, last, fin) = if complex { - let mut rbm = Rbm::complex(n, alpha, 0.01, seed); + let mut rbm = Rbm::complex(n, alpha, init, seed); if marshall { rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); } @@ -1224,7 +1256,7 @@ fn main() { } run(model, rbm, cfg, steps) } else { - let mut rbm = Rbm::new(n, alpha, 0.01, seed); + let mut rbm = Rbm::new(n, alpha, init, seed); if marshall { rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); } @@ -1235,10 +1267,11 @@ fn main() { }; if flag(args, "--json") { println!( - "{{\"complex\":{complex},\"energy\":{},\"error\":{},\"exact\":{},\"last_step\":{},\"model\":\"{kind}\",\"parameters\":{params},\"sites\":{n},\"steps\":{steps},\"variance\":{}}}", + "{{\"complex\":{complex},\"energy\":{},\"error\":{},\"exact\":{},\"exact_in_sector\":{},\"last_step\":{},\"model\":\"{kind}\",\"parameters\":{params},\"sites\":{n},\"steps\":{steps},\"variance\":{}}}", fin.energy, fin.error, reference.map_or("null".into(), |e| e.to_string()), + in_sector.map_or("null".into(), |e| e.to_string()), last.energy, fin.variance ); @@ -1248,6 +1281,9 @@ fn main() { if let Some(e) = reference { println!("exact {e:.8} relative error {:.2e}", (fin.energy - e) / e.abs()); } + if let Some(e) = in_sector { + println!("in its sector {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 f16bbd957..34a753416 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.38" +version = "0.3.39" dependencies = [ "blake3", "ed25519-dalek",