mod delayed

module delayed

Delayed-acceptance basin hopping: pay for the quench only when it can matter.

Basin hopping walks the transformed surface E~(x) = E(Q(x)), where Q is a local minimisation. Every proposal costs a quench, and measured on 38 points that is 87 to 93 per cent of the whole charged budget: 25 gradient evaluations per proposal against one for the raw energy. Two attempts to shorten the quench failed and failed informatively. Extrapolating the descent geometrically mispredicts the limit by 1e4 at the step where its rule fires. Simply running fewer steps solves 0 of 4 at six steps and 1 of 4 at fifteen, against 4 of 4 at twenty-five. The chain needs a converged energy because that is the surface it walks on.

What it does not need is a converged energy for proposals it is going to reject. Most are: the accept rate runs near one half, and the rejected half pays the same 25 evaluations as the accepted half for an answer that is discarded.

The scheme

Let E^ be a cheap surrogate for E~ and let q be the proposal. One step from x:

1. Draw y ~ q(x, .) and evaluate the surrogate. Accept the first stage with .. code-block:: text

   a1 = min(1, [q(y,x) exp(-E^(y)/T)] / [q(x,y) exp(-E^(x)/T)])

   On rejection the chain stays at x and no quench is paid.

2. Only if the first stage accepted, quench y and accept with .. code-block:: text

   a2 = min(1, exp(-[(E~(y) - E~(x)) - (E^(y) - E^(x))] / T))

Why the surrogate cannot bias the answer

The composite kernel is reversible with respect to pi ∝ exp(-E~/T) for any E^ whatever, which is the Christen and Fox delayed-acceptance argument. Write the first stage as a proposal in its own right: it draws from q and keeps with a1, so its effective proposal density from x to y != x is Q1(x,y) = q(x,y) a1(x,y). That density satisfies Q1(x,y) exp(-E^(x)/T) = Q1(y,x) exp(-E^(y)/T) by construction, since a1 is a Metropolis-Hastings ratio for the surrogate target. The second stage is then a Metropolis-Hastings step with proposal Q1 against the true target, whose ratio is

[Q1(y,x) exp(-E~(y)/T)] / [Q1(x,y) exp(-E~(x)/T)]
  = exp(-[(E~(y) - E~(x)) - (E^(y) - E^(x))] / T)

because the Q1 ratio contributes exactly the surrogate difference back. So detailed balance holds for the pair, and a poor surrogate costs acceptance rate rather than correctness: a surrogate that is constant reduces the scheme to ordinary basin hopping, and a perfect one makes the second stage accept with probability one.

This is the property the crate’s existing screen does not have. That screen decides on a partial quench with a hand-set margin and lets the trial into the chain on the partial energy, which is a different chain from the one it reports.

What the surrogate is

The run labels its own training data. Every quench it pays produces a pair of an unrelaxed structure and the energy it relaxed to, and a 38-point run produces about twelve thousand of them. The surrogate regresses the quench depth E~(x) - E(x) on features of the unrelaxed structure, so the prediction needs one energy evaluation and no gradient.

Depth rather than energy because the raw energy carries most of the signal already and is available for one evaluation; asking the model only for what the relaxation adds is asking it for the part that is actually hard.

Measured, and not yet paying

The scheme is exact and the saving is real; the surrogate is not good enough for the saving to be worth what it costs. On 38 points at 2e5:

variant	stage-1 reject	stage-2 reject	hops	acceptance	accepted moves
screen (control)	–	–	5561	0.560	3114
delayed, no abstention	0.485	0.571	10058	0.223	2243
delayed, abstention at 0.5 T	0.868	0.960	21138	0.014	296

Nearly half of all proposals avoid a quench and the hop count nearly doubles, which is the mechanism working. What it buys is spent again at the second stage, where the surrogate’s error rejects more moves than the extra proposals supply, and abstention on the predictive spread makes it worse rather than better: the posterior’s own uncertainty is not calibrated against its error, so raising the bar filtered out the cases it was right about along with the ones it was wrong about.

The bottleneck is now a specific and learnable quantity rather than a structural one. The features here are the raw energy, mean coordination and closest contact, and none of them describes how far a structure sits off the manifold of relaxed structures, which is what sets how far it will fall. The gradient at the unrelaxed point is computed and discarded by the evaluation the first stage already pays for, and for a locally quadratic basin the depth goes as |g|^2 / 2 lambda, so it is the feature the model is missing and the one that costs nothing to add.

Variables

const FEATURES: usize

Features the surrogate regresses on.

Functions

fn features(x: ArrayView1<f64>, n: usize, raw: f64) -> Array1<f64>

Cheap structural summary of an unrelaxed structure.

Coordination counts and the closest contact, which is what says whether a structure has room to relax: a proposal with an overlapping pair has a large depth ahead of it, and a proposal already near a packing has little.

fn features_orthogonal(x: ArrayView1<f64>, n: usize, raw: f64, gradient: ArrayView1<f64>, from: ArrayView1<f64>) -> Array1<f64>

As features_with_gradient, splitting the displacement from the incumbent along and across the descent direction.

A single gradient norm collapses an anisotropic quantity to one number. The depth a relaxation finds is 1/2 g^T H^-1 g, so it is set by how the displacement distributes over the curvature, not by its length. The split recovers the part of that which costs nothing: the component of d = y - x along the gradient is what simple descent undoes and predicts how much energy comes back, while the orthogonal component is the part descent does not remove and is what decides which basin the trial lands in.

This is the orthogonal-deviation reading used to visualise how far a band wanders off a reaction path, applied to a quench rather than a path: the descent direction plays the role of the path tangent.

fn features_with_gradient(x: ArrayView1<f64>, n: usize, raw: f64, gnorm: f64) -> Array1<f64>

As features, with the gradient norm at the unrelaxed point.

The feature the first version was missing, and a free one: the evaluation the first stage already pays for computes the gradient and throws it away. In a locally quadratic basin the depth a relaxation will find goes as |g|^2 / 2 lambda, so the squared norm is the leading term of the quantity being predicted rather than a proxy for it, and the linear term is carried alongside because the basin is only approximately quadratic.

Structs and Unions

struct Surrogate

A learned stand-in for the quenched energy, with its own uncertainty.

warmup: usize

Quenches required before the first stage is allowed to reject.

Enforced here rather than left to the model: Screen::predict answers from whatever posterior it holds, including a prior fitted to nothing, and the warmup it carries gates a different decision. Acting on a posterior fitted to a handful of quenches would reject good proposals for no reason. The scheme is exact either way, so waiting costs only savings not yet being made.

stage_one: usize

First-stage tests run.

stage_one_rejected: usize

First-stage rejections, each of which saved a quench.

abstained: usize

First stages skipped because the posterior was too uncertain to speak.

An abstention pays the quench the run would have paid anyway, so it is the cheap way of being unsure.

stage_two: usize

Second-stage tests run.

stage_two_rejected: usize

Second-stage rejections, which are the surrogate’s mistakes.

Implementations

impl Surrogate

Functions

fn new() -> Self

A surrogate with no evidence, which accepts every first stage.

fn observe(&mut self, x: ArrayView1<f64>, n: usize, raw: f64, quenched: f64)

Records a quench: the structure before it and the energy after.

fn observe_features(&mut self, f: ArrayView1<f64>, raw: f64, quenched: f64)

Records a quench from a feature vector the caller built.

fn observe_full(&mut self, x: ArrayView1<f64>, n: usize, raw: f64, gnorm: f64, quenched: f64)

As Surrogate::observe, with the gradient norm the prediction used.

The training features must match the prediction features or the model is fitted on one thing and consulted about another.

fn predict(&self, x: ArrayView1<f64>, n: usize, raw: f64) -> Option<f64>

Predicted quenched energy of an unrelaxed structure, or None while the posterior has too little evidence to be worth consulting.

raw is the structure’s own energy, which the caller has already paid for and which the model corrects rather than replaces.

fn predict_at(&self, x: ArrayView1<f64>, n: usize, raw: f64, tolerance: f64) -> Option<f64>

As Surrogate::predict, abstaining when the predictive spread exceeds tolerance.

Abstention is what an uncertainty is for here. The two ways of being wrong do not cost the same: abstaining pays one quench, which is what the run would have paid anyway, while guessing wrongly passes a proposal the second stage then rejects, and that costs an accepted move. So the first stage should only speak where the posterior is sharp relative to the temperature that scales the acceptance ratio.

Measured without it, the second stage rejected 57 per cent of what the first passed, and composite acceptance fell from 0.56 to 0.223: the surrogate was confidently wrong often enough to lose more moves than the extra proposals bought.

fn predict_features(&self, f: ArrayView1<f64>, raw: f64, tolerance: f64) -> Option<f64>

Prediction from a feature vector the caller built, which is how the orthogonal split is supplied.

fn predict_full(&self, x: ArrayView1<f64>, n: usize, raw: f64, gnorm: f64, tolerance: f64) -> Option<f64>

As Surrogate::predict_at, given the gradient norm at x.

fn saved_share(&self) -> f64

Quenches avoided as a share of proposals, which is the whole point.

fn seen(&self) -> usize

Quenches the surrogate has been trained on.

fn stage_one_probability(&self, surrogate_x: f64, surrogate_y: f64, t: f64) -> f64

First-stage acceptance probability for a symmetric proposal.

Symmetric because every move in this crate’s library is: a displacement, a relocation, a twin and a symmetrisation all propose y from x with the same density as x from y. Under asymmetry the ratio q(y,x) / q(x,y) multiplies in, and this returns the wrong number rather than a slightly wrong one, which is why it says so here.

fn stage_two_probability(&self, surrogate_x: f64, surrogate_y: f64, true_x: f64, true_y: f64, t: f64) -> f64

Second-stage acceptance probability, which corrects the first.

The surrogate difference is subtracted back out, so what remains is the error the surrogate made on this pair. A surrogate that was right about the difference gives one.

Traits implemented

impl Default for Surrogate