L-BFGS Two-Loop Recursion

Newton's method converges beautifully and needs the inverse Hessian — an n × n matrix, impossible to store when n is in the millions. L-BFGS approximates multiplying by that inverse while storing nothing but the last few steps the optimiser took.

It keeps two short histories, oldest first:

s_list[i] — how much the parameters moved at that step,
y_list[i] — how much the gradient changed at that step.

From those it reconstructs the product $H \nabla f$ in two passes, which is where the name comes from.

First, define $\rho_i = \dfrac{1}{y_i \cdot s_i}$ for each remembered step.

Loop one — newest to oldest. Start with q = grad, and for each i going backwards:

\alpha_i = \rho_i \,(s_i \cdot q), \qquad q \leftarrow q - \alpha_i y_i

Save every $\alpha_i$ ; the second loop needs them.

Scale. Using only the newest pair:

\gamma = \frac{s_{\text{last}} \cdot y_{\text{last}}}{y_{\text{last}} \cdot y_{\text{last}}}, \qquad r = \gamma q

Loop two — oldest to newest. For each i going forwards:

\beta = \rho_i \,(y_i \cdot r), \qquad r \leftarrow r + s_i (\alpha_i - \beta)

Task: write two_loop_recursion(grad, s_list, y_list) returning the final r, each entry rounded to 4 decimal places.

grad, and every vector in the two histories, have the same length. Both histories are the same length and hold at least one pair, ordered oldest first.
Every $y_i \cdot s_i$ is non-zero, so no $\rho$ blows up.
Return r itself. It approximates $H \nabla f$ , so the descent step the optimiser actually takes is its negative.

Nothing here ever builds a matrix — every line is a dot product or a vector add. That's the trick: the curvature information lives entirely in the history, and the cost scales with how many steps you remember rather than with how many parameters you have.

L-BFGS Two-Loop Recursion

Examples

Hints