L-BFGS

Deep dive

How it works

In: objective f and exact gradient ∇f, start x₀, history size m, gradient tolerance · Out: a stationary point x with ‖∇f(x)‖ ≤ tol · State: the last m pairs (sᵢ, yᵢ) and the current x, g · Cost: about 4mn flops per direction plus one f+∇f evaluation per line-search trial

Two-loop recursion (Nocedal 1980), with ρᵢ = 1/(yᵢᵀsᵢ), newest pair last:

q = g
for i = newest → oldest:   αᵢ = ρᵢ sᵢᵀq;   q −= αᵢ yᵢ
r = γ q,  γ = sᵀy / yᵀy (newest pair)
for i = oldest → newest:   β = ρᵢ yᵢᵀr;   r += (αᵢ − β) sᵢ
d = −r

Line search along d from step t = 1, accepting t when

• sufficient decrease (Armijo): f(x + t·d) ≤ f(x) + c₁·t·gᵀd, c₁ = 10⁻⁴;
• curvature (strong Wolfe): |g(x + t·d)ᵀd| ≤ c₂·|gᵀd|, c₂ = 0.9.

The search brackets an acceptable interval (increase t while f decreases and the slope stays negative), then "zooms" by interpolation (bisection or cubic) until both hold (Nocedal-Wright Algorithms 3.5 and 3.6). After the step, store (s, y) if sᵀy > ε·‖s‖·‖y‖ and drop the oldest pair beyond m.

Code

use std::collections::VecDeque;

pub struct Lbfgs { pub m: usize, pub gtol: f64, pub max_iter: usize, pub c1: f64, pub c2: f64 }
impl Default for Lbfgs { fn default() -> Self { Lbfgs { m: 8, gtol: 1e-8, max_iter: 1000, c1: 1e-4, c2: 0.9 } } }

pub struct Report { pub iterations: usize, pub evals: usize, pub f: f64, pub grad_norm: f64, pub converged: bool }

fn dot(a: &[f64], b: &[f64]) -> f64 { a.iter().zip(b).map(|(x, y)| x * y).sum() }
fn axpy(t: f64, d: &[f64], x: &[f64]) -> Vec<f64> { x.iter().zip(d).map(|(a, b)| a + t * b).collect() }

impl Lbfgs {
    /// Minimize f; `fg(x, grad)` writes the gradient and returns f(x).
    pub fn minimize(&self, fg: &mut impl FnMut(&[f64], &mut [f64]) -> f64, x: &mut Vec<f64>) -> Report {
        let n = x.len();
        let mut g = vec![0.0; n];
        let mut f = fg(x, &mut g);
        let mut evals = 1;
        let mut hist: VecDeque<(Vec<f64>, Vec<f64>, f64)> = VecDeque::new(); // (s, y, rho)
        for it in 0..self.max_iter {
            let gn = dot(&g, &g).sqrt();
            if gn <= self.gtol * f.abs().max(1.0) { return Report { iterations: it, evals, f, grad_norm: gn, converged: true }; }
            // two-loop recursion
            let mut q = g.clone();
            let mut alpha = vec![0.0; hist.len()];
            for (i, (s, y, rho)) in hist.iter().enumerate().rev() {
                alpha[i] = rho * dot(s, &q);
                for k in 0..n { q[k] -= alpha[i] * y[k]; }
            }
            let gamma = hist.back().map_or(1.0 / gn.max(1.0), |(s, y, _)| dot(s, y) / dot(y, y));
            for v in q.iter_mut() { *v *= gamma; }
            for (i, (s, y, rho)) in hist.iter().enumerate() {
                let beta = rho * dot(y, &q);
                for k in 0..n { q[k] += (alpha[i] - beta) * s[k]; }
            }
            let mut d: Vec<f64> = q.iter().map(|v| -v).collect();
            let mut slope = dot(&g, &d);
            if slope >= 0.0 { d = g.iter().map(|v| -v).collect(); slope = -gn * gn; hist.clear(); } // safety reset
            // strong Wolfe line search
            let Some((t, f_new, g_new, e)) = self.wolfe(fg, x, f, slope, &d) else {
                return Report { iterations: it, evals, f, grad_norm: gn, converged: false };
            };
            evals += e;
            let s: Vec<f64> = d.iter().map(|v| t * v).collect();
            let y: Vec<f64> = g_new.iter().zip(&g).map(|(a, b)| a - b).collect();
            let sy = dot(&s, &y);
            *x = axpy(1.0, &s, x);
            f = f_new;
            g = g_new;
            if sy > 1e-10 * dot(&s, &s).sqrt() * dot(&y, &y).sqrt() {
                if hist.len() == self.m { hist.pop_front(); }
                hist.push_back((s, y, 1.0 / sy));
            }
        }
        let gn = dot(&g, &g).sqrt();
        Report { iterations: self.max_iter, evals, f, grad_norm: gn, converged: false }
    }

    /// Bracket + zoom (Nocedal-Wright 3.5/3.6) with bisection-safeguarded cubic interpolation.
    fn wolfe(&self, fg: &mut impl FnMut(&[f64], &mut [f64]) -> f64, x: &[f64], f0: f64, d0: f64, d: &[f64])
        -> Option<(f64, f64, Vec<f64>, usize)> {
        let mut g = vec![0.0; x.len()];
        let mut eval = |t: f64, g: &mut Vec<f64>| { let ft = fg(&axpy(t, d, x), g); (ft, dot(g, d)) };
        let (mut t_lo, mut f_lo, mut d_lo) = (0.0, f0, d0);
        let mut t = 1.0;
        let mut t_hi = f64::INFINITY;
        let (mut f_hi, mut d_hi) = (f64::NAN, f64::NAN);
        for e in 1..=50 {
            let (ft, dt) = eval(t, &mut g);
            if ft > f0 + self.c1 * t * d0 || (ft >= f_lo && e > 1) {
                t_hi = t; f_hi = ft; d_hi = dt;             // too far: minimum is in [t_lo, t]
            } else if dt.abs() <= -self.c2 * d0 {
                return Some((t, ft, g, e));                 // both Wolfe conditions hold
            } else if dt * (t_hi - t_lo) >= 0.0 {           // also when t_hi = inf and dt > 0
                t_hi = t_lo; f_hi = f_lo; d_hi = d_lo;      // slope sign says minimum is behind
                t_lo = t; f_lo = ft; d_lo = dt;
            } else {
                t_lo = t; f_lo = ft; d_lo = dt;
            }
            t = if t_hi.is_finite() {
                // cubic interpolation between (t_lo, f_lo, d_lo) and (t_hi, f_hi, d_hi), kept inside the bracket
                let d1 = d_lo + d_hi - 3.0 * (f_lo - f_hi) / (t_lo - t_hi);
                let rad = d1 * d1 - d_lo * d_hi;
                let (a, b) = (t_lo.min(t_hi), t_lo.max(t_hi));
                let c = if rad >= 0.0 {
                    let d2 = rad.sqrt() * (t_hi - t_lo).signum();
                    t_hi - (t_hi - t_lo) * (d_hi + d2 - d1) / (d_hi - d_lo + 2.0 * d2)
                } else { f64::NAN };
                if c.is_finite() && c > a + 0.1 * (b - a) && c < b - 0.1 * (b - a) { c } else { 0.5 * (a + b) }
            } else {
                2.0 * t // still descending steeply: extrapolate
            };
        }
        None
    }
}

/// Extended Rosenbrock: sum over pairs of 100 (x_{2i+1} - x_{2i}^2)^2 + (1 - x_{2i})^2.
pub fn rosenbrock(x: &[f64], g: &mut [f64]) -> f64 {
    let mut f = 0.0;
    g.iter_mut().for_each(|v| *v = 0.0);
    for i in (0..x.len()).step_by(2) {
        let (a, b) = (x[i], x[i + 1]);
        let t = b - a * a;
        f += 100.0 * t * t + (1.0 - a) * (1.0 - a);
        g[i] += -400.0 * a * t - 2.0 * (1.0 - a);
        g[i + 1] += 200.0 * t;
    }
    f
}

Complexity and constants

• Direction: 2 loops over m pairs, each a dot product and an axpy: 4mn flops; plus the n-vector operations of the step. For expensive objectives the cost is all in f and ∇f.
• Measured with the code above (m = 8, gtol 10⁻⁸) on extended Rosenbrock from (−1.2, 1, …): 38-41 iterations and 51-52 function+gradient evaluations for n = 2, 100 and 10⁴; n = 10⁶ takes 38 iterations, 54 evaluations, 0.86 s on one Apple M4 Pro core. A diagonal quadratic with eigenvalues 1…1000 (n = 1000) needs 367 iterations, about what CG needs for κ = 1000 at that tolerance.
• Memory: 2m vectors of n doubles plus a few work vectors; n = 10⁸ with m = 10 is 16 GB, so very large problems use m = 3-5.
• Larger m helps less than one expects: on Rosenbrock, m = 3, 8 and 20 all need 36-41 iterations, while the direction cost grows linearly in m.

Efficiency

• Cost: per iteration 4mn flops for the direction plus the cost of f and ∇f, which dominates for real objectives. Measured with the code above (m = 8) on extended Rosenbrock: 38-41 iterations and 51-54 function+gradient evaluations for n = 2 up to 10⁶; n = 10⁶ takes 0.86 s on one Apple M4 Pro core. Memory is 2m vectors of n doubles; n = 10⁸ with m = 10 is 16 GB. Bound: objective and gradient evaluation; memory bandwidth for the 2m vector passes.
• Levers:
  • small m (5-10) — m = 3, 8, 20 all need 36-41 iterations on Rosenbrock (measured above), direction cost linear in m — most problems
  • strong Wolfe with c₂ = 0.9 — step 1 mostly accepted, about 1.3 evaluations per iteration (51-54 over 38-41 iterations, measured above) — always
  • rescaling variables — H₀ = γI becomes a good start — mixed units
  • m = 3-5 — fits very large n (16 GB at n = 10⁸, m = 10) — memory-limited
• Naive → standard → state of the art:
  • steepest descent with line search · O(κ) iterations (367 for κ = 1000 quadratic, comparable to CG, above)
  • full BFGS · O(n²) memory and time per step, fine to a few thousand variables
  • L-BFGS (L-BFGS-B with bounds) · 4mn flops and 2mn doubles per step

Implementation notes

• c₂ = 0.9 for quasi-Newton (loose curvature condition, step 1 usually accepted); CG-type methods use c₂ = 0.1.
• First iteration has no curvature information: scale the steepest-descent step so ‖t·d‖ is reasonable (here 1/‖g‖ for large gradients), or the first line search wastes evaluations.
• Skip bad pairs: store (s, y) only when sᵀy is positive and not tiny; otherwise H loses positive definiteness.
• Convergence tests: ‖g‖ relative to max(|f|, 1), plus a stall test on relative decrease of f over several iterations.
• Bound constraints: L-BFGS-B (Byrd-Lu-Nocedal-Zhu) uses a generalized Cauchy point and subspace minimization; projecting plain L-BFGS steps onto bounds is not correct.
• Float precision: f32 gradients limit achievable ‖g‖ to about √ε_f32 relative; keep the history and dot products in f64.
• Gradient check: compare ∇f with central finite differences on a random direction before blaming the optimizer.

Pitfalls

• Using a line search that only checks Armijo (backtracking): sᵀy can be negative, the update becomes indefinite, directions go uphill.
• Noisy objectives (Monte Carlo, mini-batches, solver tolerances): curvature pairs are garbage; increase accuracy or switch method.
• Badly scaled variables (e.g. metres and micrometres in one vector): H₀ = γI is a poor start; rescale variables.
• Not resetting history after a non-descent direction or after the problem changes (warm starts in a sequence of problems).
• Expecting it to escape saddle points or handle non-convexity gracefully: it finds stationary points.

Variants and alternatives

In production

Statistics and ML: SciPy minimize(method="L-BFGS-B"), scikit-learn LogisticRegression(solver="lbfgs"), Stan (default optimizer), PyTorch torch.optim.LBFGS, large-scale log-linear NLP training (OWL-QN, Andrew and Gao). Numerical libraries: liblbfgs (C), Optim.jl, argmin in Rust. Robotics and vision: Ceres Solver (LBFGS line search). Chemistry: geometry optimizers in quantum chemistry codes.

Implementations

Sources

• Jorge Nocedal, Updating Quasi-Newton Matrices with Limited Storage, Mathematics of Computation 35(151), 1980. paper
• Dong C. Liu, Jorge Nocedal, On the Limited Memory BFGS Method for Large Scale Optimization, Mathematical Programming 45, 1989. paper
• Jorge Nocedal, Stephen J. Wright, Numerical Optimization, 2nd ed., Springer, 2006, chapters 3 and 7. book
• Richard H. Byrd, Peihuang Lu, Jorge Nocedal, Ciyou Zhu, A Limited Memory Algorithm for Bound Constrained Optimization, SIAM Journal on Scientific Computing 16(5), 1995. paper
• Galen Andrew, Jianfeng Gao, Scalable Training of L1-Regularized Log-Linear Models, ICML 2007. paper