The secant family: BFGS, its limited memory, and the session hooks

Newton’s method prices a step at one Hessian solve; the secant family refuses to pay it. This page derives what the crate’s quasi-Newton solvers (src/qn.rs, src/lbfgs.rs) actually compute, and why the warm-restart and recognition hooks on the session are correctness- preserving.

1 The secant condition

After a step s = x_new - x_old with gradient change y = g_new - g_old, any model Hessian B worth keeping should reproduce the observed curvature: B s = y. In inverse form, with H = B^-1, the condition reads H y = s. One equation cannot determine a matrix, so each method adds a minimality principle. BFGS chooses the inverse update closest to the old H in a weighted Frobenius norm subject to the secant condition and symmetry, which lands on

H_new= (I - rho s y^T) H (I - rho y s^T) + rho s s^T, rho = 1 / (y . s).

The update preserves positive definiteness exactly when the curvature pair condition y . s > 0 holds, which the strong Wolfe curvature condition guarantees at an accepted line-search point. That is the architectural reason the crate pairs Lbfgs with LineSearch::Wolfe by default: the line search is not a convenience, it is the hypothesis of the update’s soundness.

2 The two-loop recursion

Storing H costs n²; limited-memory BFGS stores the last m pairs (s_i, y_i) and applies the update implicitly. Unrolling the BFGS formula m times around a seed matrix H0 gives the two-loop recursion (Nocedal and Wright 7.4):

q = g for i = m-1 .. 0: alpha_i= rho_i(s_i. q); q = q - alpha_iy_ir = H0 q for i = 0 .. m-1: beta = rho_i(y_i. r); r = r + (alpha_i- beta) s_idirection = -r

with the standard scalar seed H0 = (s . y) / (y . y) I from the most recent pair, a one-parameter fit of the average curvature along the last step; positivity of that scalar is again the pair condition. Every dot and axpy in the crate’s two-loop routes through the vecops seam, so the recursion is written once and the storage backend (serial, par, and by design a later device backend) is a detail behind four functions.

3 Stale pairs and the cautious update

A pair with y . s <= 0 would corrupt the inverse model, so the update skips it; the cautious knob (Li-Fukushima) additionally requires y . s >= eps norm(s)^alphabefore a pair enters the history, which keeps far-from-quadratic steps from polluting the memory. Skipping is always sound: the two-loop over any subset of valid pairs is the BFGS inverse built from exactly those pairs.

4 Session hooks: forget, watched, recognized

Solver (and anneal’s WarmLbfgs face over it) exposes three lifecycle hooks whose safety arguments are one sentence each:

• forget() clears the pair history. The next step is plain steepest descent scaled by the seed; every stored invariant is vacuous.

• minimize_watched runs a caller predicate per accepted iterate and stops early when asked. Stopping early returns an accepted iterate, so every property of the trajectory up to that point stands.

• minimize_recognized additionally lets the caller substitute a known result (energy and coordinates) when it recognizes the basin mid-descent, returning a flag that says the substitution happened. The optimizer never manufactures the substitute; it only carries the caller’s answer out, so the correctness burden sits exactly where the knowledge is. The identity and refund tests in tests/lbfgs_state.rs pin both behaviors.

5 Non-finite input is not convergence

A gradient with a NaN component has no norm worth trusting; the crate maps any non-finite gradient to an infinite norm in both the Euclidean and infinity norms (vecops::nrminf documents the rule as: a broken vector is never a small one). The test a_nan_gradient_is_not_convergence holds the L-BFGS driver to it, and the C ABI returns NaN-filled gradients on callback failure for the same reason: a consumer must never read half-written storage as a small residual.