Introduction

Every branch of science and engineering eventually reduces to solving a linear system Ax=bAx = b: a set of equations linking unknowns through a matrix AA. Climate models, circuit simulators, structural analyses, and machine-learning optimizers all spend most of their time here.

For small systems, you simply invert AA directly. But the matrices that appear in practice can have millions of rows. Direct methods are too slow and eat too much memory. Instead, you use an iterative solver: start with a guess and keep improving it until it is close enough to the true solution.

The problem is that some matrices make iterative solvers crawl. The culprit is the condition number κ(A)\kappa(A): a measure of how stretched the solution space is. A large condition number means the solver must take many tiny steps in some directions and huge ones in others — and convergence can take thousands of iterations or more.

Preconditioning is the answer. You find a matrix MM that is a cheap approximation of A−1A^{-1}, and instead of solving Ax=bAx = b directly you solve the equivalent system MAx=MbMAx = Mb. If M≈A−1M \approx A^{-1}, then MA≈IMA \approx I, which has condition number 1 — the easiest system imaginable. Even a rough approximation can shrink thousands of iterations down to tens.

Watch Convergence

Below you can run a simplified iterative solver (gradient descent on the residual) on a diagonal linear system whose condition number you control. Press Run without preconditioner to see how many steps the raw solver needs, then press Run with preconditioner to see it solve the same system in far fewer steps.

<!-- Sliders select condition number (kappa) and convergence tolerance -->
<div class="controls">
  <label>
    Condition number κ:
    <input type="range" id="kappa" min="1" max="5" step="0.5" value="3">
    <span id="kappa-val">1000</span>
  </label>
  <label>
    Tolerance:
    <input type="range" id="tol" min="1" max="4" step="1" value="2">
    <span id="tol-val">0.01</span>
  </label>
</div>
<div class="run-btns">
  <button id="run-raw" type="button">Without preconditioner</button>
  <button id="run-pre" type="button" class="accent">With preconditioner</button>
  <button id="reset-btn" type="button" class="ghost">Reset</button>
</div>
<div id="status-line" class="status-line"></div>
<div class="chart-wrap">
  <canvas id="chart" width="460" height="200"></canvas>
</div>
<div id="legend" class="legend"></div>

/* Demo styles — chart, controls, legend */
* { box-sizing: border-box; margin: 0; }
body { font-family: system-ui, sans-serif; color: #222; padding: .5rem .4rem; }
.controls { display: flex; flex-wrap: wrap; gap: .8rem; margin-bottom: .6rem; }
label { display: flex; align-items: center; gap: .4rem; font-size: .85rem; }
input[type=range] { width: 110px; }
span[id$="-val"] { font-weight: 700; min-width: 3.5rem; }
.run-btns { display: flex; 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.accent { background: #2a9d8f; border-color: #2a9d8f; }
button.ghost  { background: #fff; color: #1d3557; }
.status-line  { font-size: .9rem; font-weight: 600; min-height: 1.3em; margin-bottom: .3rem; }
.status-line.ok  { color: #0a7d33; }
.status-line.bad { color: #c92f3c; }
.chart-wrap { border: 1px solid #cdd9e3; border-radius: 8px; overflow: hidden; background: #f8fafc; }
canvas { display: block; width: 100%; height: 200px; }
.legend { display: flex; gap: 1.2rem; font-size: .8rem; margin-top: .4rem; flex-wrap: wrap; }
.legend-item { display: flex; align-items: center; gap: .3rem; }
.legend-dot  { width: 12px; height: 4px; border-radius: 2px; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));

// Build a diagonal linear system with a given condition number.
// A = diag(1, r, r^2, ..., r^(n-1)) where r = kappa^(1/(n-1))
// The Jacobi preconditioner M = diag(A)^{-1} exactly inverts the diagonal.
// M = diag(1/A_ii) — Jacobi preconditioner (exact inverse of diagonal A)

const N_DIM = 8; // system dimension (8 equations)

// Build diagonal A and right-hand side b
function buildSystem(kappa) {
  const ratio = Math.pow(kappa, 1 / (N_DIM - 1));
  const diag = [];
  for (let i = 0; i < N_DIM; i++) diag.push(Math.pow(ratio, i));
  // right-hand side: b_i = diag_i (so true solution x* = all-ones)
  const b = diag.slice();
  return { diag, b };
}

// Gradient descent for the linear system Ax = b
// Gradient descent for diag(A) x = b: x_{k+1} = x_k - step * A(Ax_k - b)
// Preconditioned variant: M*A has much lower condition number
// Preconditioned: replace A with M*A = I (since M=inv(diag(A)) and A=diag(A))
// so the preconditioned system has condition number 1 — converges in 1 step.
// Off-diagonal perturbation makes the Jacobi preconditioner approximate
// For demo realism we add a small off-diagonal perturbation so the Jacobi
// preconditioner is only approximate and we need a few iterations.
function gradientDescent(diag, b, precond, tol) {
  const n = diag.length;
  // start from the zero vector
  let x = new Array(n).fill(0);
  const residuals = [];
  const maxIter = 5000;
  // optimal step size for gradient descent: raw: 1/lambda_max (conservative, convergent); precond: 0.9 (near-unit spectrum)
  const step = precond ? 0.9 : (1.0 / diag[n - 1]);

  for (let iter = 0; iter < maxIter; iter++) {
    // compute residual: r = A*x - b
    let r = x.map((xi, i) => diag[i] * xi - b[i]);

    // add small off-diagonal coupling for realism add tiny coupling: r[i] += 0.05 * (x[i-1] or x[i+1])
    for (let i = 0; i < n; i++) {
      if (i > 0)     r[i] += 0.05 * x[i - 1];
      if (i < n - 1) r[i] += 0.05 * x[i + 1];
    }

    const relErr = Math.sqrt(r.reduce((s, v) => s + v * v, 0)) /
                   Math.sqrt(b.reduce((s, v) => s + v * v, 0));
    residuals.push(relErr);
    if (relErr < tol) break;

    // gradient step x = x - step * gradient
    if (precond) {
      // apply Jacobi preconditioner to residual M * r = r / diag (Jacobi preconditioner)
      x = x.map((xi, i) => xi - step * (r[i] / diag[i]));
    } else {
      x = x.map((xi, i) => xi - step * r[i]);
    }
  }
  return residuals;
}

// Canvas chart — log-scale residual curves
const canvas = document.getElementById('chart');
const ctx = canvas.getContext('2d');
const colors = { raw: '#e63946', pre: '#2a9d8f' };
let rawRes = null, preRes = null;

function getKappa() {
  const v = parseFloat(document.getElementById('kappa').value);
  return Math.round(Math.pow(10, v));
}
function getTol() {
  const v = parseFloat(document.getElementById('tol').value);
  return Math.pow(10, -v);
}

function drawChart() {
  const W = canvas.width, H = canvas.height;
  ctx.clearRect(0, 0, W, H);

  // draw horizontal grid lines at each decade draw log-scale grid
  ctx.strokeStyle = '#dde3ea'; ctx.lineWidth = 1;
  for (let e = 0; e >= -4; e--) {
    const y = H - ((-e) / 4) * (H - 20) - 10;
    ctx.beginPath(); ctx.moveTo(0, y); ctx.lineTo(W, y); ctx.stroke();
  }

  const allRes = [rawRes, preRes].filter(Boolean);
  if (!allRes.length) return;
  const maxLen = Math.max(...allRes.map(r => r.length));

  function drawLine(res, color) {
    if (!res || !res.length) return;
    ctx.beginPath(); ctx.strokeStyle = color; ctx.lineWidth = 2;
    res.forEach((v, i) => {
      const x = (i / Math.max(maxLen - 1, 1)) * W;
      const logV = Math.max(Math.log10(v), -4);
      const y = H - ((logV + 4) / 4) * (H - 20) - 10;
      i === 0 ? ctx.moveTo(x, y) : ctx.lineTo(x, y);
    });
    ctx.stroke();
  }
  drawLine(rawRes, colors.raw);
  drawLine(preRes, colors.pre);
}

function updateLegend() {
  const leg = document.getElementById('legend');
  leg.innerHTML = '';
  const items = [];
  if (rawRes) items.push({ color: colors.raw, label: i18nFmt('No preconditioner (%n% iter)', { n: rawRes.length }) });
  if (preRes) items.push({ color: colors.pre, label: i18nFmt('Jacobi preconditioner (%n% iter)', { n: preRes.length }) });
  items.forEach(({ color, label }) => {
    const div = document.createElement('div');
    div.className = 'legend-item';
    div.innerHTML = `<div class="legend-dot" style="background:${color}"></div><span>${label}</span>`;
    leg.appendChild(div);
  });
}

function setStatus(msg, cls) {
  const el = document.getElementById('status-line');
  el.textContent = msg; el.className = 'status-line ' + (cls || '');
}

document.getElementById('kappa').addEventListener('input', function() {
  document.getElementById('kappa-val').textContent = getKappa().toLocaleString();
});
document.getElementById('tol').addEventListener('input', function() {
  document.getElementById('tol-val').textContent = getTol().toFixed(4).replace(/\.?0+$/, '') || '0.0001';
});

document.getElementById('run-raw').addEventListener('click', function() {
  const kappa = getKappa(), tol = getTol();
  const { diag, b } = buildSystem(kappa);
  rawRes = gradientDescent(diag, b, false, tol);
  const converged = rawRes[rawRes.length - 1] < tol;
  if (converged) {
    setStatus(i18nFmt('Without preconditioner: converged in %n% iterations.', { n: rawRes.length }), 'ok');
  } else {
    setStatus(i18nFmt('Without preconditioner: did not converge in %n% iterations.', { n: rawRes.length }), 'bad');
  }
  drawChart(); updateLegend();
});

document.getElementById('run-pre').addEventListener('click', function() {
  const kappa = getKappa(), tol = getTol();
  const { diag, b } = buildSystem(kappa);
  preRes = gradientDescent(diag, b, true, tol);
  const converged = preRes[preRes.length - 1] < tol;
  if (converged) {
    setStatus(i18nFmt('With preconditioner: converged in %n% iterations.', { n: preRes.length }), 'ok');
  } else {
    setStatus(i18nFmt('With preconditioner: did not converge in %n% iterations.', { n: preRes.length }), 'bad');
  }
  drawChart(); updateLegend();
});

document.getElementById('reset-btn').addEventListener('click', function() {
  rawRes = null; preRes = null;
  drawChart(); updateLegend();
  setStatus('');
  document.getElementById('kappa').value = 3;
  document.getElementById('kappa-val').textContent = '1000';
  document.getElementById('tol').value = 2;
  document.getElementById('tol-val').textContent = '0.01';
});

// initialize slider display values on load initialize display
document.getElementById('kappa-val').textContent = getKappa().toLocaleString();
document.getElementById('tol-val').textContent = getTol().toFixed(4).replace(/\.?0+$/, '');
drawChart();
Result
Console

The preconditioner here is the Jacobi preconditioner: M=diag(A)−1M = \text{diag}(A)^{-1}, the diagonal of AA inverted. It is nearly free to compute, yet for diagonally dominant systems it can cut iterations dramatically. More powerful preconditioners — like incomplete LU factorization — work even better on harder systems.

The Real Complexity

The key quantity is the condition number κ(A)=λmax⁡/λmin⁡\kappa(A) = \lambda_{\max} / \lambda_{\min}, the ratio of the largest to the smallest eigenvalue of AA (for symmetric positive-definite matrices).

  • Gradient descent on a quadratic (equivalent to solving Ax=bAx = b) requires O(κ)O(\kappa) iterations to halve the error.
  • Conjugate Gradient (CG), the standard iterative method, is much better: it reaches ε\varepsilon relative error in O(κlog⁡(1/ε))O(\sqrt{\kappa} \log(1/\varepsilon)) iterations.
  • Even so, when κ\kappa reaches 10610^{6} or 10910^{9} — common in finite-element models — both methods stall.

A preconditioner MM replaces the system with M1/2AM1/2x′=M1/2bM^{1/2} A M^{1/2} x' = M^{1/2} b. If MM is a good approximation of A−1A^{-1}, the preconditioned condition number κ(MA)\kappa(MA) is much smaller. For CG this means:

iterations≈O ⁣(κ(MA)log⁡1ε)\text{iterations} \approx O\!\left(\sqrt{\kappa(MA)} \log\frac{1}{\varepsilon}\right)

The trade-off is building MM. The Jacobi preconditioner costs O(n)O(n) and helps modestly. Incomplete LU (ILU) and algebraic multigrid (AMG) cost more to build but can reduce κ\kappa from 10610^{6} to near 1, making the total work orders of magnitude smaller. Choosing the right preconditioner for a given matrix is still an open art — there is no single recipe that always works best. See non-convex optimization for related ideas about landscapes and step sizes.

Where It Matters

Preconditioning is invisible infrastructure. Every time a computation involves a large linear system, someone has chosen (or should have chosen) a preconditioner:

  • Finite-element analysis: stress and heat simulations on complex geometries produce sparse AA with κ\kappa reaching 10810^{8}. Multigrid preconditioners are the reason such simulations finish in minutes instead of days.
  • Computational fluid dynamics: pressure Poisson equations appear every time step. ILU preconditioners keep CFD codes tractable for aircraft and turbine design.
  • Circuit simulation: SPICE-like tools solve huge, stiff systems. Domain-specific preconditioners are part of what makes simulators like FastSPICE practical at chip scale.
  • Machine-learning optimization: training large neural networks is essentially solving many poorly conditioned systems. Adaptive optimizers like Adam act as diagonal preconditioners, which is why they outperform plain gradient descent.
  • Image reconstruction: MRI and CT reconstruction solves huge ill-conditioned least-squares problems. Regularization-based preconditioners are what let those images appear in seconds.

If you have seen dimensionality reduction, you have already seen a related idea: both reshape a hard problem into one where the important directions are easy to move along.

Conclusion

Preconditioning is a beautiful example of meta-level thinking in algorithms: instead of making the solver work harder, you transform the problem into one where the same solver works effortlessly. The raw condition number κ(A)\kappa(A) can be a million; the preconditioned condition number κ(MA)\kappa(MA) can be close to 1.

The deeper lesson is that the geometry of a problem matters as much as its size. A well-conditioned system of a billion equations is easier than a poorly conditioned system of a thousand. Next time you see an iterative algorithm struggling, ask not "how do I make it faster?" but "what is the right shape for this problem?" — the answer is almost always a preconditioner.

Share this article

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

Comments

Loading comments...

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