BFGS Demo: Curvature Estimate Speeds Descent

Introduction

You want to find the minimum of a function. Gradient descent is the obvious first tool: measure the slope, step downhill, repeat. But the slope alone is a blunt instrument. If the landscape curves steeply in one direction and gently in another, a step that works in one direction overshoots in the other. You waste iterations zigzagging when you could be gliding.

The fix is to model the curvature — the matrix of second derivatives known as the Hessian. If you know the Hessian, you can correct each step so that it lands precisely at the bottom of the local quadratic bowl. This is Newton's method, and it converges in very few iterations.

The catch is the Hessian itself. For a function of $n$ variables it is an $n \times n$ matrix, which costs $O(n^{2})$ to store and $O(n^{3})$ to invert on every step — completely impractical for modern machine learning where $n$ might be in the millions.

Quasi-Newton methods split the difference. Instead of computing the exact Hessian, they estimate it from the gradient differences that accumulate as the algorithm moves. They ask: given that the gradient changed by $\mathbf{y}_{k}$ when the position changed by $\mathbf{s}_{k}$ , what does that tell us about curvature? The answer shapes a rank-2 update to the current Hessian approximation.

The most successful version of this idea is BFGS, named after Broyden, Fletcher, Goldfarb and Shanno, who independently derived it in 1970. It is the default optimizer in scipy, the backbone of many statistical fitting routines, and the intellectual parent of L-BFGS — the optimizer behind much of modern deep learning.

Try It

Below, both gradient descent (blue) and BFGS (orange) start at the same point on a curved 2-D surface and race to the minimum. The surface is deliberately ill-conditioned — stretched much further in one direction than the other — so gradient descent zigzags while BFGS adapts.

<!-- Controls and canvas for the BFGS vs gradient-descent comparison -->
<div class="controls">
  <button id="btn-step" type="button">Step</button>
  <button id="btn-run"  type="button">Run</button>
  <button id="btn-reset" type="button" class="ghost">Reset</button>
  <span class="legend"><span class="dot gd"></span> Gradient Descent &nbsp;<span class="dot bfgs"></span> BFGS</span>
</div>
<canvas id="cv" width="420" height="300" title="Optimization paths on an ill-conditioned quadratic surface"></canvas>
<div id="info" class="info">Press Step or Run to start the race.</div>

/* Layout and styling for the optimization demo */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; background: transparent; }
.controls { display: flex; align-items: center; gap: .5rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .45; cursor: default; }
.legend { font-size: .82rem; display: flex; align-items: center; gap: .3rem; margin-left: auto; }
.dot { display: inline-block; width: 10px; height: 10px; border-radius: 50%; }
.dot.gd   { background: #3a86ff; }
.dot.bfgs { background: #fb8500; }
canvas { display: block; border: 1px solid #dde; border-radius: 8px; max-width: 100%; }
.info { font-size: .88rem; font-weight: 600; margin-top: .4rem; min-height: 1.3em; }
.info.done-gd   { color: #3a86ff; }
.info.done-bfgs { color: #e76f00; }
.info.done-both { color: #0a7d33; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// BFGS vs Gradient Descent on f(x,y) = 0.5*(A*x^2 + B*y^2).
// Condition number = B/A = 20; gradient descent zigzags, BFGS adapts quickly.

const canvas = document.getElementById('cv');
const ctx    = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;

// Surface definition
// f(x,y) = 0.5*(a*x^2 + b*y^2)  — a simple quadratic with condition number b/a
const A = 1, B = 20; // condition number = B/A = 20

function f(x, y)  { return 0.5 * (A * x * x + B * y * y); }
function gf(x, y) { return [A * x, B * y]; } // gradient of f

// Coordinate mapping: domain to canvas pixels
// Canvas: (0,0) top-left; domain: x in [-3,3], y in [-1.5,1.5]
const XLO = -3, XHI = 3, YLO = -1.5, YHI = 1.5;
function toCanvas(x, y) {
  return [
    (x - XLO) / (XHI - XLO) * W,
    (1 - (y - YLO) / (YHI - YLO)) * H
  ];
}

// Draw contour lines of the objective function
function drawContours() {
  const levels = [0.05, 0.2, 0.6, 1.4, 3.0, 6.0, 11.0];
  ctx.save();
  ctx.lineWidth = 1;
  for (const lv of levels) {
    // simple marching-squares iso-line — simple marching-squares iso-line
    const NX = 80, NY = 60;
    const dx = (XHI - XLO) / NX, dy = (YHI - YLO) / NY;
    ctx.beginPath();
    ctx.strokeStyle = `rgba(100,130,160,${0.18 + lv * 0.015})`;
    for (let i = 0; i < NX; i++) {
      for (let j = 0; j < NY; j++) {
        const x0 = XLO + i * dx, y0 = YLO + j * dy;
        const x1 = x0 + dx,      y1 = y0 + dy;
        const vals = [f(x0,y0), f(x1,y0), f(x1,y1), f(x0,y1)];
        const above = vals.map(v => v >= lv);
        if (above.every(Boolean) || above.every(v => !v)) continue;
        // linear interpolation along edges — linear interpolation along edges
        const pts = [];
        const edges = [[0,1],[1,2],[2,3],[3,0]];
        const xs = [x0,x1,x1,x0], ys = [y0,y0,y1,y1];
        for (const [a, b] of edges) {
          if (above[a] !== above[b]) {
            const t = (lv - vals[a]) / (vals[b] - vals[a]);
            pts.push(toCanvas(xs[a]+(xs[b]-xs[a])*t, ys[a]+(ys[b]-ys[a])*t));
          }
        }
        if (pts.length >= 2) {
          ctx.moveTo(...pts[0]);
          ctx.lineTo(...pts[1]);
        }
      }
    }
    ctx.stroke();
  }
  ctx.restore();
}

// Draw an optimization path with the given color
function drawPath(pts, color) {
  if (pts.length < 2) return;
  ctx.save();
  ctx.strokeStyle = color;
  ctx.lineWidth = 2;
  ctx.beginPath();
  const [cx0, cy0] = toCanvas(...pts[0]);
  ctx.moveTo(cx0, cy0);
  for (let i = 1; i < pts.length; i++) {
    const [cx, cy] = toCanvas(...pts[i]);
    ctx.lineTo(cx, cy);
  }
  ctx.stroke();
  // draw a dot at the current position — draw a dot at the current position
  const [lx, ly] = toCanvas(...pts[pts.length - 1]);
  ctx.fillStyle = color;
  ctx.beginPath();
  ctx.arc(lx, ly, 5, 0, 2 * Math.PI);
  ctx.fill();
  ctx.restore();
}

function drawMinMarker() {
  const [mx, my] = toCanvas(0, 0);
  ctx.save();
  ctx.strokeStyle = '#0a7d33';
  ctx.lineWidth = 1.5;
  const r = 6;
  ctx.beginPath();
  ctx.moveTo(mx - r, my); ctx.lineTo(mx + r, my);
  ctx.moveTo(mx, my - r); ctx.lineTo(mx, my + r);
  ctx.stroke();
  ctx.restore();
}

function redraw() {
  ctx.clearRect(0, 0, W, H);
  drawContours();
  drawMinMarker();
  drawPath(gdPath,   '#3a86ff');
  drawPath(bfgsPath, '#fb8500');
}

// Gradient Descent optimizer state
// Gradient Descent: x_{k+1} = x_k - lr * grad f(x_k)
const GD_LR = 0.045; // learning rate tuned for visual clarity
let gdPath = [];
let gdDone = false;

function gdStep() {
  if (gdDone) return;
  const [x, y] = gdPath[gdPath.length - 1];
  const [gx, gy] = gf(x, y);
  const nx = x - GD_LR * gx;
  const ny = y - GD_LR * gy;
  gdPath.push([nx, ny]);
  if (Math.sqrt(gx * gx + gy * gy) < 1e-4 || gdPath.length > 300) gdDone = true;
}

// BFGS optimizer state (inverse Hessian approximation)
// BFGS in 2-D: maintains H (inverse Hessian approximation, 2x2 matrix)
// stored as [h00, h01, h10, h11]
let bfgsPath = [];
let bfgsH    = [1, 0, 0, 1]; // start with identity (= steepest descent) — start with identity (= steepest descent)
let bfgsDone = false;
let bfgsPrevG = null;

function matVec(H, v) {
  // multiply 2x2 matrix H by vector v — multiply 2x2 matrix H by vector v
  return [
    H[0] * v[0] + H[1] * v[1],
    H[2] * v[0] + H[3] * v[1]
  ];
}

function bfgsUpdate(H, s, y) {
  // BFGS inverse-Hessian update (H is the *inverse* approx) — BFGS inverse-Hessian update (H is the *inverse* approx)
  // H_{k+1} = (I - rho*s*y^T) H_k (I - rho*y*s^T) + rho*s*s^T
  // where rho = 1 / (y^T s)
  const sy = s[0] * y[0] + s[1] * y[1]; // y^T s — inner product
  if (sy <= 0) return H; // skip update if curvature condition fails
  const rho = 1 / sy;
  // compute H*y for rank-2 update
  const Hy = matVec(H, y);
  // H_new = H - rho*(s outer Hy) - rho*(Hy outer s) + rho*(yHy + rho^-1)*s outer s
  // simplified rank-2 update:
  const yHy = y[0] * Hy[0] + y[1] * Hy[1];
  const H00 = H[0] - rho * (s[0] * Hy[0] + Hy[0] * s[0]) + rho * (rho * yHy + 1) * s[0] * s[0];
  const H01 = H[1] - rho * (s[0] * Hy[1] + Hy[0] * s[1]) + rho * (rho * yHy + 1) * s[0] * s[1];
  const H10 = H[2] - rho * (s[1] * Hy[0] + Hy[1] * s[0]) + rho * (rho * yHy + 1) * s[1] * s[0];
  const H11 = H[3] - rho * (s[1] * Hy[1] + Hy[1] * s[1]) + rho * (rho * yHy + 1) * s[1] * s[1];
  return [H00, H01, H10, H11];
}

const BFGS_LR = 1.0; // BFGS step size (controlled by backtracking line search)

function bfgsStep() {
  if (bfgsDone) return;
  const [x, y] = bfgsPath[bfgsPath.length - 1];
  const g = gf(x, y);
  // search direction: d = -H * g — search direction: d = -H * g
  const d = matVec(bfgsH, g);
  let lr = BFGS_LR;
  // simple backtracking Armijo line search — simple backtracking Armijo line search
  const f0 = f(x, y);
  const dg = d[0] * g[0] + d[1] * g[1]; // directional derivative (should be positive for descent)
  for (let i = 0; i < 20; i++) {
    const nx = x - lr * d[0], ny = y - lr * d[1];
    if (f(nx, ny) <= f0 - 1e-4 * lr * dg) break;
    lr *= 0.5;
  }
  const nx = x - lr * d[0], ny = y - lr * d[1];
  const ng = gf(nx, ny);
  const s = [nx - x, ny - y];
  const yv = [ng[0] - g[0], ng[1] - g[1]];
  bfgsH = bfgsUpdate(bfgsH, s, yv);
  bfgsPath.push([nx, ny]);
  bfgsPrevG = ng;
  if (Math.sqrt(ng[0] * ng[0] + ng[1] * ng[1]) < 1e-4 || bfgsPath.length > 100) bfgsDone = true;
}

// shared starting point for both optimizers
const START = [2.5, 1.2];

function reset() {
  gdPath   = [START.slice()];
  bfgsPath = [START.slice()];
  bfgsH    = [1, 0, 0, 1];
  bfgsPrevG = null;
  gdDone = bfgsDone = false;
  running = false;
  clearInterval(timer);
  updateInfo('Press Step or Run to start the race.', '');
  redraw();
  updateButtons();
}

function setInfo(txt, cls) {
  const el = document.getElementById('info');
  el.textContent = txt;
  el.className = 'info ' + cls;
}

function updateInfo(fallback, cls) {
  if (gdDone && bfgsDone) {
    setInfo(
      fmt('Done! GD took %v% steps, BFGS took %v2% steps.', { v: gdPath.length - 1, v2: bfgsPath.length - 1 }),
      'done-both'
    );
  } else if (gdDone) {
    setInfo(fmt('Gradient Descent converged in %v% steps. BFGS still running…', { v: gdPath.length - 1 }), 'done-gd');
  } else if (bfgsDone) {
    setInfo(fmt('BFGS converged in %v% steps. Gradient Descent still running…', { v: bfgsPath.length - 1 }), 'done-bfgs');
  } else {
    setInfo(fallback, cls);
  }
}

function updateButtons() {
  const allDone = gdDone && bfgsDone;
  document.getElementById('btn-step').disabled = allDone;
  document.getElementById('btn-run').disabled  = allDone || running;
}

function step() {
  if (!gdDone)   gdStep();
  if (!bfgsDone) bfgsStep();
  redraw();
  updateInfo(fmt('GD step %v% | BFGS step %v2%', { v: gdPath.length - 1, v2: bfgsPath.length - 1 }), '');
  updateButtons();
}

let running = false;
let timer   = null;
function run() {
  if (running) return;
  running = true;
  updateButtons();
  timer = setInterval(() => {
    if (gdDone && bfgsDone) {
      clearInterval(timer);
      running = false;
      updateInfo('', '');
      updateButtons();
      return;
    }
    step();
  }, 60);
}

document.getElementById('btn-step').addEventListener('click', step);
document.getElementById('btn-run').addEventListener('click', run);
document.getElementById('btn-reset').addEventListener('click', reset);

reset();

<!-- Controls and canvas for the BFGS vs gradient-descent comparison -->
<div class="controls">
  <button id="btn-step" type="button">Step</button>
  <button id="btn-run"  type="button">Run</button>
  <button id="btn-reset" type="button" class="ghost">Reset</button>
  <span class="legend"><span class="dot gd"></span> Gradient Descent &nbsp;<span class="dot bfgs"></span> BFGS</span>
</div>
<canvas id="cv" width="420" height="300" title="Optimization paths on an ill-conditioned quadratic surface"></canvas>
<div id="info" class="info">Press Step or Run to start the race.</div>

/* Layout and styling for the optimization demo */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; background: transparent; }
.controls { display: flex; align-items: center; gap: .5rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .45; cursor: default; }
.legend { font-size: .82rem; display: flex; align-items: center; gap: .3rem; margin-left: auto; }
.dot { display: inline-block; width: 10px; height: 10px; border-radius: 50%; }
.dot.gd   { background: #3a86ff; }
.dot.bfgs { background: #fb8500; }
canvas { display: block; border: 1px solid #dde; border-radius: 8px; max-width: 100%; }
.info { font-size: .88rem; font-weight: 600; margin-top: .4rem; min-height: 1.3em; }
.info.done-gd   { color: #3a86ff; }
.info.done-bfgs { color: #e76f00; }
.info.done-both { color: #0a7d33; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// BFGS vs Gradient Descent on f(x,y) = 0.5*(A*x^2 + B*y^2).
// Condition number = B/A = 20; gradient descent zigzags, BFGS adapts quickly.

const canvas = document.getElementById('cv');
const ctx    = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;

// Surface definition
// f(x,y) = 0.5*(a*x^2 + b*y^2)  — a simple quadratic with condition number b/a
const A = 1, B = 20; // condition number = B/A = 20

function f(x, y)  { return 0.5 * (A * x * x + B * y * y); }
function gf(x, y) { return [A * x, B * y]; } // gradient of f

// Coordinate mapping: domain to canvas pixels
// Canvas: (0,0) top-left; domain: x in [-3,3], y in [-1.5,1.5]
const XLO = -3, XHI = 3, YLO = -1.5, YHI = 1.5;
function toCanvas(x, y) {
  return [
    (x - XLO) / (XHI - XLO) * W,
    (1 - (y - YLO) / (YHI - YLO)) * H
  ];
}

// Draw contour lines of the objective function
function drawContours() {
  const levels = [0.05, 0.2, 0.6, 1.4, 3.0, 6.0, 11.0];
  ctx.save();
  ctx.lineWidth = 1;
  for (const lv of levels) {
    // simple marching-squares iso-line — simple marching-squares iso-line
    const NX = 80, NY = 60;
    const dx = (XHI - XLO) / NX, dy = (YHI - YLO) / NY;
    ctx.beginPath();
    ctx.strokeStyle = `rgba(100,130,160,${0.18 + lv * 0.015})`;
    for (let i = 0; i < NX; i++) {
      for (let j = 0; j < NY; j++) {
        const x0 = XLO + i * dx, y0 = YLO + j * dy;
        const x1 = x0 + dx,      y1 = y0 + dy;
        const vals = [f(x0,y0), f(x1,y0), f(x1,y1), f(x0,y1)];
        const above = vals.map(v => v >= lv);
        if (above.every(Boolean) || above.every(v => !v)) continue;
        // linear interpolation along edges — linear interpolation along edges
        const pts = [];
        const edges = [[0,1],[1,2],[2,3],[3,0]];
        const xs = [x0,x1,x1,x0], ys = [y0,y0,y1,y1];
        for (const [a, b] of edges) {
          if (above[a] !== above[b]) {
            const t = (lv - vals[a]) / (vals[b] - vals[a]);
            pts.push(toCanvas(xs[a]+(xs[b]-xs[a])*t, ys[a]+(ys[b]-ys[a])*t));
          }
        }
        if (pts.length >= 2) {
          ctx.moveTo(...pts[0]);
          ctx.lineTo(...pts[1]);
        }
      }
    }
    ctx.stroke();
  }
  ctx.restore();
}

// Draw an optimization path with the given color
function drawPath(pts, color) {
  if (pts.length < 2) return;
  ctx.save();
  ctx.strokeStyle = color;
  ctx.lineWidth = 2;
  ctx.beginPath();
  const [cx0, cy0] = toCanvas(...pts[0]);
  ctx.moveTo(cx0, cy0);
  for (let i = 1; i < pts.length; i++) {
    const [cx, cy] = toCanvas(...pts[i]);
    ctx.lineTo(cx, cy);
  }
  ctx.stroke();
  // draw a dot at the current position — draw a dot at the current position
  const [lx, ly] = toCanvas(...pts[pts.length - 1]);
  ctx.fillStyle = color;
  ctx.beginPath();
  ctx.arc(lx, ly, 5, 0, 2 * Math.PI);
  ctx.fill();
  ctx.restore();
}

function drawMinMarker() {
  const [mx, my] = toCanvas(0, 0);
  ctx.save();
  ctx.strokeStyle = '#0a7d33';
  ctx.lineWidth = 1.5;
  const r = 6;
  ctx.beginPath();
  ctx.moveTo(mx - r, my); ctx.lineTo(mx + r, my);
  ctx.moveTo(mx, my - r); ctx.lineTo(mx, my + r);
  ctx.stroke();
  ctx.restore();
}

function redraw() {
  ctx.clearRect(0, 0, W, H);
  drawContours();
  drawMinMarker();
  drawPath(gdPath,   '#3a86ff');
  drawPath(bfgsPath, '#fb8500');
}

// Gradient Descent optimizer state
// Gradient Descent: x_{k+1} = x_k - lr * grad f(x_k)
const GD_LR = 0.045; // learning rate tuned for visual clarity
let gdPath = [];
let gdDone = false;

function gdStep() {
  if (gdDone) return;
  const [x, y] = gdPath[gdPath.length - 1];
  const [gx, gy] = gf(x, y);
  const nx = x - GD_LR * gx;
  const ny = y - GD_LR * gy;
  gdPath.push([nx, ny]);
  if (Math.sqrt(gx * gx + gy * gy) < 1e-4 || gdPath.length > 300) gdDone = true;
}

// BFGS optimizer state (inverse Hessian approximation)
// BFGS in 2-D: maintains H (inverse Hessian approximation, 2x2 matrix)
// stored as [h00, h01, h10, h11]
let bfgsPath = [];
let bfgsH    = [1, 0, 0, 1]; // start with identity (= steepest descent) — start with identity (= steepest descent)
let bfgsDone = false;
let bfgsPrevG = null;

function matVec(H, v) {
  // multiply 2x2 matrix H by vector v — multiply 2x2 matrix H by vector v
  return [
    H[0] * v[0] + H[1] * v[1],
    H[2] * v[0] + H[3] * v[1]
  ];
}

function bfgsUpdate(H, s, y) {
  // BFGS inverse-Hessian update (H is the *inverse* approx) — BFGS inverse-Hessian update (H is the *inverse* approx)
  // H_{k+1} = (I - rho*s*y^T) H_k (I - rho*y*s^T) + rho*s*s^T
  // where rho = 1 / (y^T s)
  const sy = s[0] * y[0] + s[1] * y[1]; // y^T s — inner product
  if (sy <= 0) return H; // skip update if curvature condition fails
  const rho = 1 / sy;
  // compute H*y for rank-2 update
  const Hy = matVec(H, y);
  // H_new = H - rho*(s outer Hy) - rho*(Hy outer s) + rho*(yHy + rho^-1)*s outer s
  // simplified rank-2 update:
  const yHy = y[0] * Hy[0] + y[1] * Hy[1];
  const H00 = H[0] - rho * (s[0] * Hy[0] + Hy[0] * s[0]) + rho * (rho * yHy + 1) * s[0] * s[0];
  const H01 = H[1] - rho * (s[0] * Hy[1] + Hy[0] * s[1]) + rho * (rho * yHy + 1) * s[0] * s[1];
  const H10 = H[2] - rho * (s[1] * Hy[0] + Hy[1] * s[0]) + rho * (rho * yHy + 1) * s[1] * s[0];
  const H11 = H[3] - rho * (s[1] * Hy[1] + Hy[1] * s[1]) + rho * (rho * yHy + 1) * s[1] * s[1];
  return [H00, H01, H10, H11];
}

const BFGS_LR = 1.0; // BFGS step size (controlled by backtracking line search)

function bfgsStep() {
  if (bfgsDone) return;
  const [x, y] = bfgsPath[bfgsPath.length - 1];
  const g = gf(x, y);
  // search direction: d = -H * g — search direction: d = -H * g
  const d = matVec(bfgsH, g);
  let lr = BFGS_LR;
  // simple backtracking Armijo line search — simple backtracking Armijo line search
  const f0 = f(x, y);
  const dg = d[0] * g[0] + d[1] * g[1]; // directional derivative (should be positive for descent)
  for (let i = 0; i < 20; i++) {
    const nx = x - lr * d[0], ny = y - lr * d[1];
    if (f(nx, ny) <= f0 - 1e-4 * lr * dg) break;
    lr *= 0.5;
  }
  const nx = x - lr * d[0], ny = y - lr * d[1];
  const ng = gf(nx, ny);
  const s = [nx - x, ny - y];
  const yv = [ng[0] - g[0], ng[1] - g[1]];
  bfgsH = bfgsUpdate(bfgsH, s, yv);
  bfgsPath.push([nx, ny]);
  bfgsPrevG = ng;
  if (Math.sqrt(ng[0] * ng[0] + ng[1] * ng[1]) < 1e-4 || bfgsPath.length > 100) bfgsDone = true;
}

// shared starting point for both optimizers
const START = [2.5, 1.2];

function reset() {
  gdPath   = [START.slice()];
  bfgsPath = [START.slice()];
  bfgsH    = [1, 0, 0, 1];
  bfgsPrevG = null;
  gdDone = bfgsDone = false;
  running = false;
  clearInterval(timer);
  updateInfo('Press Step or Run to start the race.', '');
  redraw();
  updateButtons();
}

function setInfo(txt, cls) {
  const el = document.getElementById('info');
  el.textContent = txt;
  el.className = 'info ' + cls;
}

function updateInfo(fallback, cls) {
  if (gdDone && bfgsDone) {
    setInfo(
      fmt('Done! GD took %v% steps, BFGS took %v2% steps.', { v: gdPath.length - 1, v2: bfgsPath.length - 1 }),
      'done-both'
    );
  } else if (gdDone) {
    setInfo(fmt('Gradient Descent converged in %v% steps. BFGS still running…', { v: gdPath.length - 1 }), 'done-gd');
  } else if (bfgsDone) {
    setInfo(fmt('BFGS converged in %v% steps. Gradient Descent still running…', { v: bfgsPath.length - 1 }), 'done-bfgs');
  } else {
    setInfo(fallback, cls);
  }
}

function updateButtons() {
  const allDone = gdDone && bfgsDone;
  document.getElementById('btn-step').disabled = allDone;
  document.getElementById('btn-run').disabled  = allDone || running;
}

function step() {
  if (!gdDone)   gdStep();
  if (!bfgsDone) bfgsStep();
  redraw();
  updateInfo(fmt('GD step %v% | BFGS step %v2%', { v: gdPath.length - 1, v2: bfgsPath.length - 1 }), '');
  updateButtons();
}

let running = false;
let timer   = null;
function run() {
  if (running) return;
  running = true;
  updateButtons();
  timer = setInterval(() => {
    if (gdDone && bfgsDone) {
      clearInterval(timer);
      running = false;
      updateInfo('', '');
      updateButtons();
      return;
    }
    step();
  }, 60);
}

document.getElementById('btn-step').addEventListener('click', step);
document.getElementById('btn-run').addEventListener('click', run);
document.getElementById('btn-reset').addEventListener('click', reset);

reset();

Result

Console

Press Step to advance one iteration at a time, or Run to watch the full race. Notice how BFGS builds a better model of the local shape after just a few steps and then takes a nearly direct path, while gradient descent keeps bouncing off the steep walls. The step count shown in the corner makes the difference vivid.

The Real Complexity

What exactly does BFGS do on each step, and why does it work?

The secant condition. After moving from $\mathbf{x}_{k}$ to $\mathbf{x}_{k+1}$ , we have two gradient vectors. Their difference $\mathbf{y}_{k} = \nabla f(\mathbf{x}_{k+1}) - \nabla f(\mathbf{x}_{k})$ and the step $\mathbf{s}_{k} = \mathbf{x}_{k+1} - \mathbf{x}_{k}$ must satisfy the secant condition: the new Hessian approximation $B_{k+1}$ should predict the observed gradient change, i.e. $B_{k+1} \mathbf{s}_{k} = \mathbf{y}_{k}$ .

The minimal update. There are infinitely many matrices satisfying the secant condition. BFGS picks the one closest to $B_{k}$ in a natural matrix norm — a symmetric rank-2 correction:

B_{k+1} = B_{k} - \frac{B_{k}\mathbf{s}_{k}\mathbf{s}_{k}^{\top}B_{k}}{\mathbf{s}_{k}^{\top}B_{k}\mathbf{s}_{k}} + \frac{\mathbf{y}_{k}\mathbf{y}_{k}^{\top}}{\mathbf{y}_{k}^{\top}\mathbf{s}_{k}}

If $B_{k}$ is positive-definite and the curvature condition $\mathbf{y}_{k}^{\top}\mathbf{s}_{k} > 0$ holds (guaranteed by a Wolfe-condition line search), then $B_{k+1}$ stays positive-definite — meaning the search direction always points downhill.

Convergence. On smooth, strongly convex functions, BFGS achieves superlinear convergence: the error shrinks faster than any geometric rate, approaching quadratic convergence near the minimum. Gradient descent, by contrast, converges only linearly, at a rate governed by the condition number $\kappa$ of the Hessian.

L-BFGS. For large-scale problems, storing the full $n \times n$ matrix is too costly. Limited-memory BFGS stores only the last $m$ pairs $(\mathbf{s}_{k}, \mathbf{y}_{k})$ — typically $m \in \{5, \dots, 20\}$ — and reconstructs the product $B_{k}^{-1}\mathbf{g}$ implicitly in $O(mn)$ time. This is the optimizer used to train many large language models and neural networks, and it is also related to algorithms used in non-convex optimization.

Where It Matters

Quasi-Newton methods are the workhorse of smooth, medium-to-large-scale optimization:

Machine learning: L-BFGS trains logistic regression, SVMs, and small-to-medium neural networks. PyTorch's torch.optim.LBFGS is a standard option for full-batch settings. Some large language model fine-tuning recipes use L-BFGS variants for their fast convergence.
Statistical fitting: maximum-likelihood estimation in mixed models, survival analysis, and structural equation models almost always uses a quasi-Newton solver internally — the user types fit() and BFGS is doing the work.
Scientific computing: scipy.optimize.minimize defaults to L-BFGS-B (a bound-constrained variant) for unconstrained smooth problems. Molecular dynamics force-field minimization, geophysical inversion, and optimal control all rely on it.
Engineering design: aerodynamic shape optimization, antenna design, and PDE-constrained optimization use adjoint methods to compute gradients cheaply, then feed those gradients to L-BFGS.

The key insight connecting all of these is the same: whenever you can compute gradients but not second derivatives, and the problem is smooth enough for a local quadratic model to be useful, quasi-Newton is the right tool. See also gradient descent's failure modes for the cases where it breaks down.

Conclusion

The core insight of BFGS is almost disarmingly simple: every time you take a step, the change in gradient tells you something about curvature. Collect those clues, keep the cheapest update that respects them, and you get a Hessian approximation that improves with every iteration — without ever computing a single second derivative.

The result is an algorithm that achieves superlinear convergence on smooth problems, costs only $O(n^{2})$ per step in its full form, and shrinks to $O(mn)$ per step in L-BFGS. It is the default optimizer in scipy, the go-to for statistical fitting, and the backbone of neural network training when full-batch optimization is viable.

The next time a model trains in seconds instead of hours, there is a good chance a rank-2 matrix update is quietly doing the heavy lifting behind the scenes.

Comments

Loading comments...

https://www.kipuhub.com/en/article/bfgs-quasi-newton/Content licensed under CC BY-NC 4.0.

Cookies & privacy

BFGS & Quasi-Newton Methods

Approximating curvature, one gradient at a time