|
|
|
|
@ -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<i32>) -> (Vec<usize>, Vec<u32>) {
|
|
|
|
|
let n = self.n;
|
|
|
|
|
let states: Vec<usize> = (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<i32>) -> Option<f64> {
|
|
|
|
|
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<f64> = (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<f64> = 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<usize> = (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<i32>, sym: &Symmetry) -> Option<f64> {
|
|
|
|
|
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<f64> = 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::<f64>(), start.iter().map(|v| v * v).sum::<f64>());
|
|
|
|
|
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<usize> = (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<T> {
|
|
|
|
|
energies: Vec<T>,
|
|
|
|
|
derivs: Vec<T>,
|
|
|
|
|
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<T: Amplitude> Vmc<T> {
|
|
|
|
|
@ -925,7 +1285,16 @@ impl<T: Amplitude> Vmc<T> {
|
|
|
|
|
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<T: Amplitude> Vmc<T> {
|
|
|
|
|
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<T: Amplitude> Vmc<T> {
|
|
|
|
|
let top = logs.iter().copied().fold(f64::NEG_INFINITY, f64::max);
|
|
|
|
|
let raw: Vec<f64> = 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<T> {
|
|
|
|
|
@ -1011,21 +1380,31 @@ impl<T: Amplitude> Vmc<T> {
|
|
|
|
|
let thin = self.cfg.thin.max(1) * n;
|
|
|
|
|
let rbm = &self.rbm;
|
|
|
|
|
let model = &self.model;
|
|
|
|
|
let run = |chain: &mut (Vec<i8>, Rng)| -> (Vec<T>, Vec<T>, f64) {
|
|
|
|
|
let guide = (self.cfg.guide > 0.0).then(|| {
|
|
|
|
|
let logs: Vec<f64> = self.chains.iter().map(|(s, _)| rbm.log_abs_cached(s, &rbm.cache(s))).filter(|l| l.is_finite()).collect();
|
|
|
|
|
Guide { c: logs.iter().sum::<f64>() / logs.len().max(1) as f64, ln_floor: ln(self.cfg.guide) }
|
|
|
|
|
});
|
|
|
|
|
let run = |chain: &mut (Vec<i8>, Rng)| -> (Vec<T>, Vec<T>, Vec<f64>, 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<T>, Vec<T>, f64)> = if threads <= 1 || cfg!(target_arch = "wasm32") {
|
|
|
|
|
let results: Vec<(Vec<T>, Vec<T>, Vec<f64>, 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<T: Amplitude> Vmc<T> {
|
|
|
|
|
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<T: Amplitude>(b: &Batch<T>, 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::<f64>().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<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<T: Amplitude>(rbm: &Rbm<T>, 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<T: Amplitude>(rbm: &Rbm<T>, s: &mut [i8], rng: &mut Rng, moves: usize, exchange: bool, guide: Option<Guide>) -> (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<usize> = if exchange {
|
|
|
|
|
let i = rng.below(n);
|
|
|
|
|
@ -1277,13 +1673,25 @@ fn metropolis<T: Amplitude>(rbm: &Rbm<T>, 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<i8> {
|
|
|
|
|
(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}");
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|