Estimate noise and stop on the floor

Computational noise sets the smallest finite-difference step that is worth taking, and the gradient size at which a tighter tolerance is no longer a stationary point. The estimator is Moré and Wild, SIAM J. Sci. Comput. 33, 1292 (2011), https://doi.org/10.1137/100786125. The steps are Moré and Wild, ACM Trans. Math. Softw. 38, 19 (2012), https://doi.org/10.1145/2168773.2168777. The difference table itself is ECnoise.

Estimate

use ndarray::array;
use rgmin::{ecnoise, sample_gradient_noise, sample_noise, NoiseFlag, DEFAULT_M};

let values = [1.0, 1.01, 1.0, 1.01, 1.0, 1.01, 1.0];
let table = ecnoise(&values).unwrap();
assert_eq!(table.ok, NoiseFlag::Detected);

let energy = sample_noise(&objective, x.view(), direction.view(), 1e-4, DEFAULT_M).unwrap();
let grad = sample_gradient_noise(
    &objective,
    x.view(),
    direction.view(),
    1e-4,
    DEFAULT_M,
    seed,
)
.unwrap();

h is the spacing between samples. m = 0 selects DEFAULT_M (8). inform 2 retries at 100 h (next_spacing); inform 3 retries at h / 100. NoiseEstimate::declared(eps) records a value the backend already knows and skips the samples.

sample_gradient_noise projects the force onto one unit vector w. Component i of w comes from a counter-based stream keyed by (seed, i), not from a thread-local generator. The normal draw uses the pure-Rust libm logarithm and cosine, so the bits do not depend on the platform libm. The dot product is taken in index order. With the par feature the m + 1 evaluations run in parallel and are stored at the sample index, so the noise bits match a one-thread run.

sample_gradient_noise_k repeats that projection on k directions keyed by (seed, j) and returns the root-mean-square of the detected eps_g². k = 1 is the single-direction estimate. For any noise covariance Σ, one uniform direction has E[eps_g²] = tr(Σ) / n, because the sphere average of w wᵀ is I / n. The identity does not need Σ to be isotropic.

Choose the step

An ECnoise table is accepted because its high differences are noise. Those differences are O(eps / h^k), so they are not a derivative. |f''| and |T| come from Moré and Wild, ACM Trans. Math. Softw. 38, 19 (2012), Algorithm 5.1, https://doi.org/10.1145/2168773.2168777.

estimate_second(oracle, x, p, eps_f) builds a central second difference at its own spacing. The first try is ha = eps_f^{1/4}. It is accepted when |Δ²f| ≥ 100 eps_f and |f(x±h) − f(x)| ≤ 0.1 max(|f(x)|, |f(x±h)|). One retry uses hb = (eps_f / µa)^{1/4}. That retry is kept when it meets the same two tests, and also when |µa − µb| ≤ µb / 2. Otherwise the estimate is None.

estimate_third applies the same algorithm to the directional derivative φ(t) = ∇f(x + t û) · û, with eps_g in place of eps_f. The second difference of φ is the third-derivative scale |T| of f along û. The two tests, the spacing, and the one retry are unchanged apart from that substitution.

forward_step(eps_f, second, point_scale) returns 8^{1/4} sqrt(eps_f / |f''|). That step minimizes (h/2)|f''| + sqrt(2) eps_f / h.

gradient_dr(eps_g, third, point_scale) returns 2 sqrt(eps_g / |T|). That step minimizes (dr/2)|T| + 2 eps_g / dr.

point_scale is max_i |x_i|. The clamp is [max(1e-16, 4 ulp(point_scale)), 1]. A step below one ulp of a coordinate near 10 Å (1.8e-15) leaves x + h == x, so the lower end follows the coordinate. A zero noise returns the lower end. A zero scale, and a missing estimate (None), with positive noise return the upper end.

Stop on the floor

The floor is off until Control::noise_floor is set. gtol is tested first, so a noise-free quadratic that meets gtol reports Status::Tolerance.

let mut control = Control::default()
    .with_noise_floor(NoiseFloor::new(grad.eps).with_c(3.0));
control.gtol = 1e-12;

NoiseFloor::new sets c = 3. The iteration stops when the stationarity measure is below c * sqrt(n_free) * eps_g and still at or above gtol. n_free is the number of coordinates that measure runs over. Frozen or fixed coordinates are left out. On the session path it is the length of the horizontal or projected gradient, or the count of box coordinates with positive width.

One projection has variance tr(Σ) / n, so the noise in ||g||_2 is sqrt(n) times eps_g. Comparing ||g||_2 with eps_g alone sits that factor too low. The rms of a noisy gradient whose true norm is r is sqrt(r² + tr Σ). That equals the threshold c sqrt(tr Σ) when r = sqrt(c² − 1) sqrt(tr Σ). At c = 3 the 50% firing radius is sqrt(c² − 1) sqrt(tr Σ) ≈ 2.83 sqrt(tr Σ). With an exact eps_g, a radius of 2 sqrt(tr Σ) already fires on 93% of draws at n = 3 and on all of them at n = 21 and n = 114. A measured eps_g sits lower. That spread, not the chi-square tail, is what a tighter radius is covering. On a 100-fold spread of component σ, one direction’s factor sqrt(n wᵀΣw / tr Σ) runs from about 0.56 to 1.45. At k = 1 the measured floor at 1.7 sqrt(tr Σ) fired on 112 of 200 seeds. At k = 8 it cleared 90%. At a true gradient of 5 sqrt(tr Σ) the measured floor fired on at most 5 of 200 seeds, for k = 1 and k = 8, isotropic and anisotropic.

eps_f is the energy noise. sample_value_gradient_noise measures it on the same stencil as eps_g, through one value_and_gradient per sample. NoiseFloor::with_eps_f records it. The line search then stops with Status::NoiseFloor, and keeps the pre-step point, when no trial decreases the energy by at least c * eps_f. A gradient floor does not see that stall: the gradient can sit far above 3 sqrt(n) eps_g while every trial is inside the energy roundoff.

The stop is Status::NoiseFloor. Spending the iteration budget above both tests is Status::Budget. SCG’s own step and objective tests, taken while ||g|| is still at or above gtol, are Status::StepTolerance.

The bench noises the forces only on the Gaussian series, so a budget run there descends on an exact energy. A second series adds an independent energy draw of standard deviation σ · 1 Å at the same point. It is not the gradient of a potential. f32 rounds both the energy and the forces. calls includes the stencil. sample_calls is that stencil alone, m + 1 = 9 gradients for each spacing tried.

The all-ones direction is a translation. LJ38’s energy and forces do not change under it, so an f32 table is all zeros and ECnoise reports StepTooSmall. The measurement then uses the first coordinate axis, which changes a pair distance. On that axis at h = 1e-4, LJ38 f32 gives eps_g = 2.77e-6 and eps_f = 5.32e-7, both detected. Without eps_f the line search first sticks at ||g|| = 1.57e-2, while 3 sqrt(114) eps_g = 8.9e-5, and the run spends the budget (25590 minimization calls, abs_de 1.0e-6). With eps_f the median f32 floor stops in 63 calls, 18 of them the stencil, at abs_de 4.6e-6. The baseline gtol = 0.01 stops in 57 calls at abs_de 2.3e-7. The call count is within twice baseline. The energy is not: the test stops on the first search whose best decrease misses c * eps_f, and the following search still clears it. The geometry is not translation invariant. Its f32 stencil on the all-ones direction detects eps_g = 3.17e-7 and eps_f = 7.26e-8. Without eps_f the search first sticks at ||g|| = 7.98e-5 and spends the budget. With eps_f the median stop is 24 calls, 9 of them the stencil, at abs_de 1.3e-7. Baseline stops in 10 calls at abs_de 2.2e-6. The energy is closer. The calls are 2.4 times baseline.

LJ38 1e-2-floor seed 27 is the censored Gaussian row. The measured eps_g was 2.47e-3 against a true component σ of 1e-2, so 3 sqrt(114) eps_g = 0.079 stayed under ||g|| (minimum 0.082, final 0.101, fmax 0.023). The energy is exact. The minimization used 22029 calls and the stencil another 18, 22047 in total, abs_de 9.9e-6.

A fixed gtol set under the floor keeps iterating after the gradient is noise. On LJ38 (n = 114) and on Lennard-Jones at the eOn heptamer coordinates (lj-pt7-geom, 7 free atoms, n = 21) the scaled floor stops in tens of calls once sampling is counted. The tight tolerance 0.1 sqrt(n) σ is censored at 2000 iterations on every force-only Gaussian seed. At σ ≤ 3e-4 the two endpoints agree to about 1e-6 eV. At σ = 1e-2 the early stop sits higher, near 1e-3 eV, while the budget-exhausted run is closer to the noise-free minimum, because that run still sees an exact energy. The baseline gtol = 0.01 stops early at low noise and burns the same budget as the tight run once 0.01 is below the noisy gradient. The geometry run is Lennard-Jones on those coordinates. It is not a Pt heptamer. rgpot is not linked.