Quasi-Newton Methods and the L-BFGS Two-Loop Recursion

Newton’s method converges in very few iterations and costs too much per iteration to use on a model with millions of parameters. Quasi-Newton methods keep the structure of the Newton step while replacing the Hessian with an approximation built from the gradients already computed along the trajectory. The full-memory version of this idea still needs a matrix, so the limited-memory variant stores a handful of vector pairs instead and reconstructs the effect of the matrix on a vector on demand. The result is the method of choice for smooth, moderate-dimensional optimization, and its failure modes are exactly the ways in which the secant information can be wrong.

What the Newton step buys, and its price

The Newton step solves H p = −∇f(x) for the search direction p, then takes a step of length controlled by a line search or trust region. Near a strict local minimum with a Lipschitz-continuous Hessian, the full step is accepted and the error contracts quadratically: the number of correct digits doubles per iteration. The price is the Hessian: forming it costs one second derivative per variable pair, storing it costs n² floating-point values, which is 800 MB at n = 10⁴ and 8 × 10¹² bytes at n = 10⁶, and factorising it costs O(n³). The scaling is what rules the method out for statistical models, where n is the parameter count.

Quasi-Newton methods never form the Hessian. They maintain a matrix B intended to approximate ∇²f(x), or its inverse H, and they update it using only the gradient differences that the iteration produces anyway. Convergence is superlinear rather than quadratic: the error ratio between successive steps tends to zero while the exponent stays at one, which in practice lands between the linear rate of gradient descent and the quadratic rate of Newton, at a cost per iteration of O(n²) for the matrix update in the dense case.

The secant equation and the curvature it encodes

The requirement that the model reproduce the observed change in gradient is the secant equation

B_{k+1} s_k = y_k, with s_k = x_{k+1} − x_k and y_k = ∇f(x_{k+1}) − ∇f(x_k).

One matrix must satisfy n scalar equations, so the equation alone does not determine the update, and the different quasi-Newton methods are the different choices of additional structure. The second requirement is that the update preserve positive definiteness, which a line search can enforce through the curvature condition sᵀy > 0: for a strictly convex objective the mean-value form of y guarantees it, and for a general objective a step satisfying the Wolfe conditions guarantees it. When sᵀy ≤ 0, the observed curvature contradicts the model’s positive-definiteness assumption. The two responses are to skip the update and keep the previous matrix, or to damp the difference y toward the multiple of s that makes the condition hold.

The BFGS update and its inverse form

The update named after Broyden, Fletcher, Goldfarb and Shanno is the rank-two correction that is closest to B_k among those satisfying the secant equation, in a weighted Frobenius sense:

B_{k+1} = B_k − (B_k s sᵀB_k)/(sᵀB_k s) + (y yᵀ)/(yᵀs).

The first correction subtracts the component of the old model along the new step direction, and the second adds the newly observed curvature. The inverse form, obtained by the Sherman–Morrison–Woodbury identity, replaces B with H and swaps the roles of s and y, so the search direction is a matrix–vector product rather than a linear solve. Because the update depends on the current point only through s and y, the pair (s, y) is a complete record of that iteration’s information, which is what makes the limited-memory version possible.

An important structural fact is that the update is invariant to a global rescaling of the objective: scaling f by a constant scales all gradients and hence all y by the same constant, and the resulting direction from the inverse form is unchanged. Gradient descent has no such invariance, which is one reason quasi-Newton methods need much less parameter tuning of the step size.

Limited memory: what is stored and the two-loop recursion

The limited-memory variant keeps the m most recent pairs (sᵢ, yᵢ), with m conventionally between 3 and 20, and discards the rest. Applying the implied inverse Hessian to a vector requires no stored matrix. The recursion starts from an initial inverse Hessian of the form H₀ = γI with γ = sᵀy / yᵀy, the scaling that makes the first pair consistent, then applies the rank-two corrections in reverse chronological order, the two-loop structure, and finishes by scaling with H₀.

def lbfgs_direction(g, S, Y, gamma):
    """Two-loop recursion: apply the implied inverse Hessian to g."""
    q = list(g)
    rho = [1.0 / sum(si * yi for si, yi in zip(S[i], Y[i])) for i in range(len(S))]
    alpha = [0.0] * len(S)
    for i in range(len(S) - 1, -1, -1):      # backward pass: the s,y pairs
        alpha[i] = rho[i] * sum(si * qi for si, qi in zip(S[i], q))
        q = [qi - alpha[i] * yi for qi, yi in zip(q, Y[i])]
    r = [gamma * qi for qi in q]             # initial Hessian scaling
    for i in range(len(S)):                  # forward pass
        beta = rho[i] * sum(yi * ri for yi, ri in zip(Y[i], r))
        r = [ri + (alpha[i] - beta) * si for ri, si in zip(r, S[i])]
    return r

Each inner product is O(n) and there are 2m of them, so one direction costs O(mn) and the storage is 2mn floating-point values. At n = 10⁶ and m = 10, that is 2 × 10⁷ doubles, about 160 MB, large but bounded and independent of the history length. The scaling γ = sᵀy/yᵀy is the only place where the objective’s magnitude enters, and setting m = 0 reduces the method to preconditioned steepest descent with that scale.

The line search has to earn its keep

The pairs are usable only if sᵀy > 0, and the line search is what ensures it. The Wolfe conditions require sufficient decrease, f(x + αp) ≤ f(x) + c₁α ∇f(x)ᵀp, with c₁ = 10⁻⁴ conventionally, together with the curvature condition ∇f(x + αp)ᵀp ≥ c₂ ∇f(x)ᵀp with c₂ = 0.9 conventionally. The first condition prevents steps so long that the model’s decline was imaginary, and the second prevents steps so short that sᵀy is dominated by noise. Both constants are conventions rather than theorems; the pair (10⁻⁴, 0.9) is the standard choice for quasi-Newton methods, while (10⁻⁴, 0.1) is common for nonlinear conjugate gradients where an accurate line search is less important.

The consequence is that the iteration and the line search are coupled in both directions. The line search needs a descent direction, which positive definiteness of B provides; the update needs sᵀy > 0, which the curvature condition provides. Breaking either link breaks the method in a way that shows up as a slow crawl rather than a failure: with sᵀy ≤ 0 and the update applied anyway, the matrix can become indefinite, the direction stops being a descent direction, and every subsequent line search terminates at its minimum step length.

When the curvature condition fails, skipping the update is the simple response and damping is the more effective one. Damping replaces the observed y with a combination of y and Bs chosen so that the modified pair satisfies the curvature condition with a margin, which keeps the direction usable and still extracts part of the observed curvature information. The price is that the update no longer solves the secant equation exactly, so the superlinear convergence argument does not cover that iteration. A related guard is the limited-memory habit of dropping the oldest pair when the buffer is full: an old pair describes curvature at a point the iteration has left, and retaining it biases the model toward the geometry of an earlier phase of the run.

Where it fails

The most consequential failure mode is the stochastic gradient. The secant pair y = g_{k+1} − g_k then contains the difference of two noise realisations in addition to the curvature term, so the implied curvature is dominated by noise when the step is small, exactly when the Hessian estimate should be most accurate. Methods that attempt the combination require variance-reduced gradient differences or increasing batch sizes, and the resulting algorithms are considerably more delicate than either ingredient alone. A second failure is ill-conditioning without curvature information: because B is updated from single steps, a nearly singular direction can persist for many iterations and produce enormous steps that the line search rejects, which manifests as many function evaluations per iteration and no progress.

Scale and memory limits are the practical bounds. Dense BFGS is unusable past a few thousand variables because of the n² storage, and the limited-memory form trades approximation quality for bounded cost: a too-small m loses the curvature information from longer-scale structure and degrades toward steepest descent, while a large m costs more per iteration for a marginally better direction. Finally, the positive-definiteness requirement is not an assumption that holds on non-convex problems but a restriction the method imposes on its own model: in a region of negative curvature, the true Hessian would suggest moving along the negative eigenvector, while the quasi-Newton model insists on a descent direction based on an approximation that cannot represent that structure. Trust-region variants that allow indefinite B handle this case at the cost of solving a subproblem, and the methods that explicitly compute negative curvature are covered in escaping saddle points in nonconvex optimization. For large models where neither dense nor limited-memory curvature is affordable, the per-coordinate approximations of adaptive step sizes are the surviving option, at the cost of ignoring all coupling between coordinates.