wai-quantum 0.3.37: complex neural quantum states, and referees that cannot stop on an excited state #176
13 changed files with 614 additions and 248 deletions
2
wai/wai-quantum-serve/Cargo.lock
generated
2
wai/wai-quantum-serve/Cargo.lock
generated
|
|
@ -589,7 +589,7 @@ dependencies = [
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.36"
|
||||
version = "0.3.37"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
2
wai/wai-quantum/Cargo.lock
generated
2
wai/wai-quantum/Cargo.lock
generated
|
|
@ -415,7 +415,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.36"
|
||||
version = "0.3.37"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
|
|
@ -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"
|
||||
|
|
@ -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"]
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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<String> = 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<T: std::str::FromStr>(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<T: Amplitude>(model: SpinModel, rbm: Rbm<T>, 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::<f64>() / 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<String> = 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<C> = 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);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -509,7 +509,7 @@ pub fn fci(fd: &Fcidump) -> Option<FciResult> {
|
|||
// a little of every other determinant, so no symmetry hides the ground
|
||||
// state.
|
||||
let start: Vec<f64> = (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 })
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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<usize>]) -> (Vec<bool>, 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<f64>) {
|
||||
pub(crate) fn lanczos(apply: &dyn Fn(&[f64], &mut [f64]), start: &[f64], krylov: usize, restarts: usize, tol: f64, early: bool) -> (f64, Vec<f64>) {
|
||||
let dim = start.len();
|
||||
let norm = |v: &[f64]| v.iter().map(|x| x * x).sum::<f64>().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<f64> {
|
|||
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<f64> = (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 θᵀθ.
|
||||
|
|
|
|||
|
|
@ -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,45 @@
|
|||
//! 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.
|
||||
//! 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.
|
||||
//! 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 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.
|
||||
|
|
@ -63,8 +117,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 +172,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<usize> {
|
||||
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<bool> {
|
||||
(0..lx * ly).map(|i| (i % lx + i / lx).is_multiple_of(2)).collect()
|
||||
|
|
@ -129,6 +222,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 +310,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 +318,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<f64>`), so `ψ` is positive up to the optional sign, or complex
|
||||
/// (`Rbm<C>`), 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<T = f64> {
|
||||
pub n: usize,
|
||||
pub m: usize,
|
||||
pub a: Vec<f64>,
|
||||
pub b: Vec<f64>,
|
||||
pub a: Vec<T>,
|
||||
pub b: Vec<T>,
|
||||
/// `W[j·n + i]`.
|
||||
pub w: Vec<f64>,
|
||||
pub w: Vec<T>,
|
||||
pub sublattice: Option<Vec<bool>>,
|
||||
}
|
||||
|
||||
|
|
@ -266,50 +434,75 @@ impl Rng {
|
|||
}
|
||||
}
|
||||
|
||||
impl Rbm {
|
||||
/// `α·n` hidden units, weights uniform in `±scale`, biases zero.
|
||||
impl Rbm<f64> {
|
||||
/// 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<C> {
|
||||
/// 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<C> {
|
||||
Rbm::random(n, alpha, scale, seed)
|
||||
}
|
||||
}
|
||||
|
||||
impl<T: Amplitude> Rbm<T> {
|
||||
fn random(n: usize, alpha: usize, scale: f64, seed: u64) -> Rbm<T> {
|
||||
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<bool>) -> Rbm {
|
||||
pub fn with_marshall_sign(mut self, sublattice: Vec<bool>) -> Rbm<T> {
|
||||
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<f64> {
|
||||
fn theta(&self, s: &[i8]) -> Vec<T> {
|
||||
(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 +514,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<T: Amplitude>(model: &SpinModel, rbm: &Rbm<T>, 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 +623,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 +634,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<T = f64> {
|
||||
pub model: SpinModel,
|
||||
pub rbm: Rbm,
|
||||
pub rbm: Rbm<T>,
|
||||
pub cfg: VmcConfig,
|
||||
chains: Vec<(Vec<i8>, 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<T> {
|
||||
weights: Vec<f64>,
|
||||
energies: Vec<f64>,
|
||||
derivs: Vec<f64>,
|
||||
energies: Vec<T>,
|
||||
derivs: Vec<T>,
|
||||
acceptance: f64,
|
||||
}
|
||||
|
||||
impl Vmc {
|
||||
pub fn new(model: SpinModel, rbm: Rbm, cfg: VmcConfig) -> Vmc {
|
||||
impl<T: Amplitude> Vmc<T> {
|
||||
pub fn new(model: SpinModel, rbm: Rbm<T>, cfg: VmcConfig) -> Vmc<T> {
|
||||
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 +678,22 @@ impl Vmc {
|
|||
v
|
||||
}
|
||||
|
||||
fn batch(&mut self) -> Batch {
|
||||
fn batch(&mut self) -> Batch<T> {
|
||||
if self.cfg.exact { self.enumerate() } else { self.sample() }
|
||||
}
|
||||
|
||||
fn enumerate(&self) -> Batch {
|
||||
fn enumerate(&self) -> Batch<T> {
|
||||
let n = self.model.n;
|
||||
assert!(n <= 24, "exact enumeration takes at most 24 spins");
|
||||
let configs: Vec<usize> = (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<usize>, logs: &mut [f64], energies: &mut [f64], derivs: &mut [f64]| {
|
||||
let work = |range: core::ops::Range<usize>, logs: &mut [f64], energies: &mut [T], derivs: &mut [T]| {
|
||||
let start = range.start;
|
||||
for idx in range {
|
||||
let x = configs[idx];
|
||||
|
|
@ -516,8 +710,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 +736,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<T> {
|
||||
let n = self.model.n;
|
||||
let p = self.rbm.params();
|
||||
let nc = self.chains.len();
|
||||
|
|
@ -551,10 +745,10 @@ impl Vmc {
|
|||
let thin = self.cfg.thin.max(1) * n;
|
||||
let rbm = &self.rbm;
|
||||
let model = &self.model;
|
||||
let run = |chain: &mut (Vec<i8>, Rng)| -> (Vec<f64>, Vec<f64>, f64) {
|
||||
let run = |chain: &mut (Vec<i8>, Rng)| -> (Vec<T>, Vec<T>, 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 +759,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<f64>, Vec<f64>, f64)> = if threads <= 1 || cfg!(target_arch = "wasm32") {
|
||||
let results: Vec<(Vec<T>, Vec<T>, f64)> = if threads <= 1 || cfg!(target_arch = "wasm32") {
|
||||
self.chains.iter_mut().map(run).collect()
|
||||
} else {
|
||||
let per = nc.div_ceil(threads);
|
||||
|
|
@ -601,8 +795,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 +807,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<T: Amplitude>(b: &Batch<T>, 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<f64> {
|
||||
fn sr_direction<T: Amplitude>(b: &Batch<T>, p: usize, shift: f64) -> Vec<T> {
|
||||
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<f64> {
|
||||
{
|
||||
// (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<T: Amplitude>(o: &[T], eps: &[T], k: usize, p: usize, shift: f64) -> Vec<T> {
|
||||
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<f64>, k: usize, p: usize, shift: f64) -> Vec<f64> {
|
||||
{
|
||||
// δ = 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<T: Amplitude>(o: &[T], eps: Vec<T>, k: usize, p: usize, shift: f64) -> Vec<T> {
|
||||
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<T: Amplitude>(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 +958,7 @@ fn start_configuration(n: usize, magnetisation: Option<i32>, 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<T: Amplitude>(rbm: &Rbm<T>, 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 +974,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 +1117,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<C> = (0..k * p).map(|_| C { re: rng.uniform() - 0.5, im: rng.uniform() - 0.5 }).collect();
|
||||
let eps: Vec<C> = (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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
|
||||
|
|
|
|||
2
wai/wai-rs/Cargo.lock
generated
2
wai/wai-rs/Cargo.lock
generated
|
|
@ -3900,7 +3900,7 @@ dependencies = [
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.36"
|
||||
version = "0.3.37"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
|
|
@ -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<T: Amplitude>(model: SpinModel, rbm: Rbm<T>, 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());
|
||||
|
|
|
|||
2
wai/wai-web/Cargo.lock
generated
2
wai/wai-web/Cargo.lock
generated
|
|
@ -734,7 +734,7 @@ dependencies = [
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.36"
|
||||
version = "0.3.37"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
Loading…
Reference in a new issue