Introduction

Somewhere inside every physics simulation, every 3-D renderer, every finite-element model, and many machine-learning pipelines sits a monster linear system: millions of equations, millions of unknowns, and the answers needed right now. Direct methods like Gaussian elimination are out — they'd chew through memory and time in cubic fashion. You need something smarter.

Conjugate gradient (CG) is that something. Invented by Magnus Hestenes and Eduard Stiefel in 1952, it solves a symmetric positive-definite system Ax = b by taking a sequence of steps, each perfectly crafted so it never undoes the progress of any earlier step. In theory it converges in at most n iterations for an n-variable problem. In practice, on the sparse systems that arise in real engineering, it often converges far sooner.

The trick is not just "follow the slope downhill." It is "follow a direction that is A-conjugate to everything before it" — a generalized orthogonality that keeps the search directions from interfering. Each step lands you on the exact minimum of the objective in the subspace explored so far, and those gains are never lost.

This makes CG one of the great algorithms of applied mathematics: provably finite, memory-light (you only ever need a handful of vectors), and fast enough to power daily computation at a scale naive methods cannot touch.

See It Converge

Below is a 2-D quadratic bowl — the simplest case of Ax = b where the solution is the bowl's minimum. Press Step (CG) to advance conjugate gradient one iteration at a time, or Run GD to watch plain gradient descent take the same journey. CG will reach the minimum in at most 2 steps; gradient descent zigzags for many more.

<p class="hint">A 2-D quadratic bowl: solving <b>Ax = b</b> is the same as finding its minimum. Watch <b>Conjugate Gradient</b> land there in exactly 2 steps, while plain <b>Gradient Descent</b> zigzags for many more.</p>
<canvas id="cvs" width="340" height="260"></canvas>
<div class="info" id="info">Press a button to start.</div>
<div class="btns">
  <button id="stepCG" type="button">Step (CG)</button>
  <button id="runGD" type="button">Run GD</button>
  <button id="reset" type="button" class="ghost">Reset</button>
</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; }
.hint { font-size: .88rem; color: #444; margin: 0 0 .6rem; line-height: 1.45; }
canvas { display: block; border-radius: 10px; border: 1px solid #cdd9e3; background: #f8fafc; max-width: 100%; }
.info { font-size: .92rem; font-weight: 600; min-height: 1.4em; margin: .5rem 0; color: #1d3557; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; }
button { font: 600 14px system-ui, sans-serif; padding: .45rem .9rem;
         border: 1px solid #1d3557; background: #1d3557; color: #fff;
         border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Conjugate Gradient on a 2-D quadratic f(x) = 0.5 x^T A x - b^T x.
// A = [[4,1],[1,3]], b = [1,2]  (SPD, solution x* = A^{-1}b).
// CG converges in EXACTLY 2 steps for a 2-D SPD system.
// Gradient descent is run for comparison (fixed step size alpha=0.1).

const A = [[4,1],[1,3]];
const b = [1,2];

// Exact solution: A x* = b  ->  x* = A^{-1} b
// det(A) = 12-1 = 11,  A^{-1} = (1/11)[[3,-1],[-1,4]]
// x* = (1/11)[[3*1+(-1)*2],[(-1)*1+4*2]] = (1/11)[1,7]
const xStar = [1/11, 7/11];

function mv(M, v) { return [M[0][0]*v[0]+M[0][1]*v[1], M[1][0]*v[0]+M[1][1]*v[1]]; }
function dot(u, v) { return u[0]*v[0]+u[1]*v[1]; }
function add(u, v, s=1) { return [u[0]+s*v[0], u[1]+s*v[1]]; }
function scale(s, v) { return [s*v[0], s*v[1]]; }

// Canvas mapping: world coords [-0.5, 1.5] -> canvas [20, 320] x [240, 20]
const W = 340, H = 260;
const wx = [-0.5, 1.5], wy = [-0.3, 1.3];
function toCanvas(x, y) {
  return [
    20 + (x - wx[0]) / (wx[1] - wx[0]) * (W - 40),
    H - 20 - (y - wy[0]) / (wy[1] - wy[0]) * (H - 40)
  ];
}

// f(x,y) = 0.5*(4x^2+2xy+3y^2) - x - 2y
function f(x, y) { return 0.5*(4*x*x+2*x*y+3*y*y) - x - 2*y; }

// Draw filled contour lines using brute-force pixel shading
const cvs = document.getElementById('cvs');
const ctx = cvs.getContext('2d');

function drawBackground() {
  const img = ctx.createImageData(W, H);
  // Compute range of f for normalizing
  const fMin = f(xStar[0], xStar[1]);
  const fMax = 3.5; // visually chosen to show ~5 contour bands
  for (let py = 0; py < H; py++) {
    for (let px = 0; px < W; px++) {
      const x = wx[0] + (px - 20) / (W - 40) * (wx[1] - wx[0]);
      const y = wy[0] + (H - 20 - py) / (H - 40) * (wy[1] - wy[0]);
      const v = f(x, y);
      const t = Math.max(0, Math.min(1, (v - fMin) / (fMax - fMin)));
      // Color: deep blue (close to min) -> light gray (far)
      const r = Math.round(215 + (232-215)*t);
      const g = Math.round(220 + (244-220)*t);
      const bl = Math.round(235 + (249-235)*t);
      const i = (py * W + px) * 4;
      img.data[i] = r; img.data[i+1] = g; img.data[i+2] = bl; img.data[i+3] = 255;
    }
  }
  ctx.putImageData(img, 0, 0);
  // Draw contour lines
  const levels = [0.05, 0.2, 0.5, 1.0, 1.8, 2.8];
  const fMinV = f(xStar[0], xStar[1]);
  ctx.strokeStyle = 'rgba(100,130,160,0.35)';
  ctx.lineWidth = 0.7;
  for (const lv of levels) {
    const target = fMinV + lv;
    // March over pixel borders to trace iso-level
    for (let py = 1; py < H-1; py++) {
      for (let px = 1; px < W-1; px++) {
        const x = wx[0] + (px - 20) / (W - 40) * (wx[1] - wx[0]);
        const y = wy[0] + (H - 20 - py) / (H - 40) * (wy[1] - wy[0]);
        const dx = (wx[1] - wx[0]) / (W - 40);
        const dy = (wy[1] - wy[0]) / (H - 40);
        const f00 = f(x, y) - target;
        const f10 = f(x+dx, y) - target;
        const f01 = f(x, y+dy) - target;
        if ((f00 < 0) !== (f10 < 0) || (f00 < 0) !== (f01 < 0)) {
          ctx.fillStyle = 'rgba(100,130,160,0.4)';
          ctx.fillRect(px, py, 1, 1);
        }
      }
    }
  }
  // Mark the minimum
  const [mx, my] = toCanvas(xStar[0], xStar[1]);
  ctx.beginPath(); ctx.arc(mx, my, 5, 0, 2*Math.PI);
  ctx.fillStyle = '#e63946'; ctx.fill();
  ctx.strokeStyle = '#fff'; ctx.lineWidth = 1.5; ctx.stroke();
}

// State
let cgPath = null, gdPath = null;

function initCG() {
  // Start at x0 = [1.2, 0.1]
  const x0 = [1.2, 0.1];
  const r0 = add(b, mv(A, x0), -1); // r = b - Ax
  const p0 = [...r0];
  return { x: x0, r: r0, p: p0, step: 0, done: false };
}

function initGD() {
  return { x: [1.2, 0.1], step: 0, done: false, alpha: 0.12, path: [[1.2, 0.1]] };
}

let cg = initCG();
let gdState = initGD();
cgPath = [cg.x.slice()];
gdPath = [gdState.x.slice()];

function drawPath(path, color, dashed=false) {
  if (path.length < 2) return;
  ctx.beginPath();
  ctx.setLineDash(dashed ? [5,3] : []);
  ctx.strokeStyle = color;
  ctx.lineWidth = 2;
  const [sx, sy] = toCanvas(path[0][0], path[0][1]);
  ctx.moveTo(sx, sy);
  for (let i = 1; i < path.length; i++) {
    const [cx2, cy2] = toCanvas(path[i][0], path[i][1]);
    ctx.lineTo(cx2, cy2);
  }
  ctx.stroke();
  ctx.setLineDash([]);
  // Draw dots
  for (const pt of path) {
    const [cx2, cy2] = toCanvas(pt[0], pt[1]);
    ctx.beginPath(); ctx.arc(cx2, cy2, 3.5, 0, 2*Math.PI);
    ctx.fillStyle = color; ctx.fill();
  }
}

function render() {
  drawBackground();
  drawPath(gdPath, '#f4a261', true);
  drawPath(cgPath, '#1d3557');
  // Start dot
  const [sx, sy] = toCanvas(1.2, 0.1);
  ctx.beginPath(); ctx.arc(sx, sy, 5, 0, 2*Math.PI);
  ctx.fillStyle = '#2a9d8f'; ctx.fill();
}

function stepCG() {
  if (cg.done) { setInfo(fmt('CG already converged in %v% step(s) — exact!', { v: cgPath.length-1 })); return; }
  const { x, r, p } = cg;
  const Ap = mv(A, p);
  const alpha = dot(r, r) / dot(p, Ap);
  const xNew = add(x, p, alpha);
  const rNew = add(r, Ap, -alpha);
  const beta = dot(rNew, rNew) / dot(r, r);
  const pNew = add(rNew, p, beta);
  cg = { x: xNew, r: rNew, p: pNew, step: cg.step + 1 };
  cgPath.push(xNew.slice());
  const err = Math.hypot(xNew[0]-xStar[0], xNew[1]-xStar[1]);
  if (err < 1e-9 || cg.step >= 2) {
    cg.done = true;
    setInfo(fmt('CG converged in %step% step(s). Distance to solution: %toExponential%', { step: cg.step, toExponential: err.toExponential(2) }));
    console.log('CG final x:', xNew, '| x*:', xStar, '| error:', err);
  } else {
    setInfo(fmt('CG step %step%: x = [%map%]', { step: cg.step, map: xNew.map(v=>v.toFixed(4)).join(', ') }));
  }
  render();
}

function runGD() {
  // Run gradient descent for up to 50 steps
  let { x, alpha } = gdState;
  gdPath = [x.slice()];
  for (let i = 0; i < 50; i++) {
    const grad = add(mv(A, x), b, -1); // gradient = Ax - b
    x = add(x, grad, -alpha);
    gdPath.push(x.slice());
    const err = Math.hypot(x[0]-xStar[0], x[1]-xStar[1]);
    if (err < 1e-4) break;
  }
  const steps = gdPath.length - 1;
  const finalErr = Math.hypot(x[0]-xStar[0], x[1]-xStar[1]);
  setInfo(fmt('GD took %steps% steps (err %toExponential%). CG needs at most 2.', { steps: steps, toExponential: finalErr.toExponential(2) }));
  console.log('GD steps:', steps, '| final x:', x, '| error:', finalErr);
  render();
}

function setInfo(t) { document.getElementById('info').textContent = t; }

document.getElementById('stepCG').addEventListener('click', stepCG);
document.getElementById('runGD').addEventListener('click', runGD);
document.getElementById('reset').addEventListener('click', () => {
  cg = initCG(); gdState = initGD();
  cgPath = [cg.x.slice()]; gdPath = [gdState.x.slice()];
  setInfo('Press a button to start.'); render();
});

render();
Result
Console

Notice that each CG step lands at the exact minimum along its search direction, and the next direction is chosen to be A-conjugate (non-interfering) with the previous one. Two independent directions in 2-D means two steps — done. Gradient descent, forced to always go straight downhill, bounces back and forth across the valley.

The Real Complexity

Conjugate gradient sits in P — the class of problems solvable in polynomial time. Specifically, for an n×nn \times n symmetric positive-definite (SPD) matrix:

  • Exact convergence in n steps. Each iteration produces a search direction A-conjugate to all previous ones. Because there are at most n such directions in Rn\mathbb{R}^n, CG reaches the exact answer in at most n steps. This is a proven upper bound, not a heuristic.
  • Cost per step: O(n)O(n) for sparse matrices. The dominant operation each iteration is a matrix-vector product Ap. When A is sparse (most entries zero, as in physics grids or graph problems), this is proportional to the number of non-zeros — far cheaper than the O(n2)O(n^{2}) dense case.
  • Total work: O(n⋅nnz)O(n \cdot nnz) in the worst case, vs O(n3)O(n^{3}) for Gaussian elimination. For large sparse systems this is the difference between feasible and impossible.
  • Eigenvalue clustering accelerates convergence. The actual convergence rate depends on the condition number κ(A)=λmax⁡/λmin⁡\kappa(A) = \lambda_{\max}/\lambda_{\min}. If eigenvalues cluster tightly (small κ\kappa), CG converges in far fewer than n steps. Preconditioners are designed to reshape A's spectrum toward clustering.
  • Status: proven polynomial, no open questions. Unlike P vs NP, CG's complexity is settled: it solves SPD linear systems in polynomial time and this has been known since Hestenes and Stiefel's 1952 paper.

The catch: CG only applies to symmetric positive-definite systems. For general linear systems, variants like GMRES or BiCGSTAB generalize the idea at higher memory cost.

Where It Matters

Whenever a large, sparse, symmetric positive-definite system appears — and this shape arises constantly — conjugate gradient is the default weapon of choice:

  • Finite-element analysis (FEA): simulating stress in a bridge, airflow over a wing, or heat through an engine produces SPD systems with millions of unknowns. CG with a good preconditioner makes these tractable on a laptop.
  • Image reconstruction (CT/MRI): tomographic reconstruction amounts to solving a large linear system from projections. Iterative solvers like CG let you trade reconstruction time against image quality.
  • Machine learning: L2-regularized least-squares (ridge regression) and certain neural-network sub-problems reduce to SPD linear systems. CG solves them without ever forming the full matrix.
  • Graph problems: the graph Laplacian is SPD, so electrical-flow calculations, semi-supervised learning, and spectral clustering all benefit from fast Laplacian solvers — many of which use CG internally.
  • Preconditioned CG (PCG) is so ubiquitous it appears as a standard black-box in libraries like MATLAB, SciPy, PETSc, and Eigen.

Understanding CG means understanding the heart of modern scientific computing — and seeing how a clever geometric idea (A-conjugate directions) turns an intractable brute-force sweep into an elegant, finite descent. Compare the brute-force exponential search in P vs NP with what happens when the problem structure lets you pick smart directions.

Conclusion

Conjugate gradient is a reminder that the right geometry can turn a seemingly endless search into a finite, elegant march. By choosing each new direction to be A-conjugate to all previous ones, the algorithm guarantees that no progress is ever wasted — and that for an n-variable system, n perfectly aimed steps are always enough.

In practice the power is even greater: sparse matrices, clustered eigenvalues, and smart preconditioners routinely bring that bound down from n to something far smaller, making problems with millions of variables routine. Not every hard problem enjoys such a lucky structure — but when it does, CG is the proof that polynomial time can be beautiful.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/conjugate-gradient/Content licensed under CC BY-NC 4.0.