Introduction

Every square matrix hides a set of special numbers called eigenvalues. They govern how the matrix stretches or rotates space, and they appear throughout science: the resonant frequencies of a structure, the energy levels of a quantum system, the ranking scores behind PageRank, the principal components in data analysis.

For decades, computing even one eigenvalue of a large matrix was a painful chore. Then in 1961, independently and almost simultaneously, J.H. Wilkinson, John Francis (in England) and Vera Kublanovskaya (in the USSR) discovered an iteration so elegant it seems like magic: factor the matrix, swap the two pieces, multiply them back together, and repeat.

After enough repetitions the matrix converges to a triangular form, and the eigenvalues simply read off the diagonal. This is the QR algorithm โ€” named after the two matrices in its core factorization โ€” and it was later voted one of the top ten algorithms of the 20th century by the journal Computing in Science and Engineering.

Try It: Watch Eigenvalues Emerge

Below is a live 3ร—3 matrix undergoing QR iteration. At each step the matrix A is factored into Q (orthogonal) and R (upper-triangular), then the pieces are recombined as RยทQ to form the next iterate.

<div class="demo-wrap">
  <div class="controls">
    <button id="stepBtn" type="button">Step</button>
    <button id="runBtn" type="button">Run</button>
    <button id="resetBtn" type="button" class="ghost">Reset</button>
    <span class="iter-label">Iteration: <span id="iterCount">0</span></span>
  </div>
  <div class="panels">
    <div class="panel">
      <div class="panel-title">Current matrix A<sub>k</sub></div>
      <div id="matA" class="matrix"></div>
    </div>
    <div class="panel">
      <div class="panel-title">Off-diagonal magnitude (log scale)</div>
      <div id="chart" class="chart">
        <canvas id="cvs" width="256" height="90"></canvas>
      </div>
    </div>
  </div>
  <div class="eigen-row">
    <span class="eigen-label">Diagonal (โ†’ eigenvalues):</span>
    <span id="diagVals" class="diag-vals"></span>
  </div>
  <div class="eigen-row">
    <span class="eigen-label">True eigenvalues:</span>
    <span id="trueVals" class="diag-vals true"></span>
  </div>
  <div id="status" class="status"></div>
</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; font-size: 14px; }
.demo-wrap { padding: 4px 0; }
.controls { display: flex; align-items: center; gap: 8px; flex-wrap: wrap; margin-bottom: 10px; }
button { font: 600 13px system-ui; padding: 6px 14px; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 7px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
.iter-label { font-size: 13px; color: #555; margin-left: 4px; }
.panels { display: flex; gap: 14px; flex-wrap: wrap; align-items: flex-start; margin-bottom: 10px; }
.panel { background: #f4f7fa; border: 1px solid #d0d9e4; border-radius: 8px; padding: 10px 12px; }
.panel-title { font-size: 12px; font-weight: 600; color: #4a6080; margin-bottom: 8px; }
.matrix { display: grid; grid-template-columns: repeat(3, 68px); gap: 3px; }
.cell { width: 68px; height: 34px; display: flex; align-items: center; justify-content: center;
        font: 500 12px ui-monospace, monospace; border-radius: 5px; transition: background 0.3s; }
.cell.diag { background: #dbeafe; border: 1px solid #93c5fd; color: #1e3a5f; }
.cell.off  { background: #fef3c7; border: 1px solid #fcd34d; color: #78350f; }
.cell.off.tiny { background: #f0fdf4; border: 1px solid #86efac; color: #166534; }
.chart { background: #fff; border: 1px solid #d0d9e4; border-radius: 6px; padding: 4px; }
canvas { display: block; }
.eigen-row { display: flex; align-items: baseline; gap: 8px; margin: 3px 0; flex-wrap: wrap; }
.eigen-label { font-size: 12px; font-weight: 600; color: #4a6080; white-space: nowrap; }
.diag-vals { font: 500 12px ui-monospace, monospace; color: #1d3557; }
.diag-vals.true { color: #15803d; }
.status { font-size: 13px; font-weight: 600; min-height: 1.3em; margin-top: 6px; color: #15803d; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// QR Algorithm Demo
// Uses Gram-Schmidt QR decomposition, then A_{k+1} = R * Q.
// The symmetric matrix A = [[4,2,1],[2,5,3],[1,3,6]] has eigenvalues โ‰ˆ 1.9213, 3.7302, 9.3485 (verified via bisection on the characteristic poly).
// As iteration proceeds, the off-diagonal entries decay to 0 and the diagonal converges to the true eigenvalues.

// โ”€โ”€ Matrix helpers (flat row-major, n=3) โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
const N = 3;
const zeros = () => Array(N*N).fill(0);
const idx = (r,c) => r*N+c;
const get = (M,r,c) => M[idx(r,c)];
const set = (M,r,c,v) => { M[idx(r,c)] = v; };

function matMul(A,B){
  const C = zeros();
  for(let r=0;r<N;r++) for(let c=0;c<N;c++){
    let s=0; for(let k=0;k<N;k++) s += get(A,r,k)*get(B,k,c);
    set(C,r,c,s);
  }
  return C;
}

// Gram-Schmidt QR decomposition: A = Q * R
function qrDecomp(A){
  const cols = [];
  for(let c=0;c<N;c++){
    let v = Array.from({length:N}, (_,r) => get(A,r,c));
    for(let j=0;j<c;j++){
      const dot = v.reduce((s,x,i) => s + x*cols[j][i], 0);
      v = v.map((x,i) => x - dot*cols[j][i]);
    }
    const norm = Math.sqrt(v.reduce((s,x) => s+x*x, 0));
    cols.push(norm < 1e-15 ? v : v.map(x => x/norm));
  }
  const Q = zeros();
  for(let r=0;r<N;r++) for(let c=0;c<N;c++) set(Q,r,c, cols[c][r]);
  // R = Q^T A  (upper-triangular since Q is orthogonal)
  const QT = zeros();
  for(let r=0;r<N;r++) for(let c=0;c<N;c++) set(QT,r,c, get(Q,c,r));
  return { Q, R: matMul(QT, A) };
}

// Off-diagonal Frobenius norm โ€” drives to 0 as matrix becomes triangular
function offDiagNorm(A){
  let s=0;
  for(let r=0;r<N;r++) for(let c=0;c<N;c++) if(r!==c) s += get(A,r,c)**2;
  return Math.sqrt(s);
}

// โ”€โ”€ True eigenvalues via bisection on the characteristic polynomial โ”€โ”€โ”€โ”€โ”€โ”€
// For A = [[4,2,1],[2,5,3],[1,3,6]]: char poly = det(A - ฮปI) = -ฮปยณ + 15ฮปยฒ - 60ฮป + 67
// Roots bracket from scanning: near 1.92, 3.73, 9.35
function charPoly(l){
  const a=4-l, b=2, c=1, d=2, e=5-l, f=3, g=1, h=3, ii=6-l;
  return a*(e*ii - f*h) - b*(d*ii - f*g) + c*(d*h - e*g);
}
function bisect(lo, hi){
  for(let i=0;i<80;i++){
    const mid=(lo+hi)/2;
    if(charPoly(lo)*charPoly(mid)<=0) hi=mid; else lo=mid;
  }
  return (lo+hi)/2;
}
// Brackets determined by sign-change scan of charPoly
const TRUE_EV = [
  bisect(1.5, 2.5),   // โ‰ˆ 1.9213
  bisect(3.0, 4.5),   // โ‰ˆ 3.7302
  bisect(8.5, 10.0)   // โ‰ˆ 9.3485
];

// โ”€โ”€ State โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
const A0 = [4,2,1, 2,5,3, 1,3,6];
let A = [...A0];
let iter = 0;
let history = [offDiagNorm(A)];
let running = false;
let runTimer = null;

// โ”€โ”€ DOM โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
const matDiv   = document.getElementById('matA');
const diagDiv  = document.getElementById('diagVals');
const trueDiv  = document.getElementById('trueVals');
const iterEl   = document.getElementById('iterCount');
const statusEl = document.getElementById('status');
const cvs      = document.getElementById('cvs');
const ctx      = cvs.getContext('2d');

function renderMatrix(){
  matDiv.innerHTML = '';
  for(let r=0;r<N;r++) for(let c=0;c<N;c++){
    const d = document.createElement('div');
    const v = get(A,r,c);
    if(r===c){
      d.className = 'cell diag';
    } else {
      d.className = Math.abs(v) < 1e-3 ? 'cell off tiny' : 'cell off';
    }
    d.textContent = v.toFixed(3);
    matDiv.appendChild(d);
  }
}

function renderDiag(){
  const diag = [get(A,0,0), get(A,1,1), get(A,2,2)];
  const sorted = [...diag].sort((a,b)=>a-b);
  diagDiv.textContent = sorted.map(v=>v.toFixed(4)).join(',  ');
  trueDiv.textContent = TRUE_EV.map(v=>v.toFixed(4)).join(',  ');
}

function renderChart(){
  const W=256, H=90, pad=6;
  ctx.clearRect(0,0,W,H);
  if(history.length < 2){
    ctx.fillStyle='#aaa'; ctx.font='11px system-ui';
    ctx.fillText('Step to see convergenceโ€ฆ',8,H/2+4); return;
  }
  const posH = history.filter(v=>v>1e-15);
  const maxLog = Math.log10(Math.max(...posH)) + 0.5;
  const minLog = -10;
  const xScale = (W-2*pad) / Math.max(history.length-1,1);
  const yScale = (H-2*pad) / (maxLog - minLog);

  // Convergence line
  ctx.strokeStyle = '#86efac'; ctx.lineWidth=1; ctx.setLineDash([3,3]);
  const y0 = H - pad - (Math.log10(1e-6) - minLog)*yScale;
  if(y0>pad && y0<H-pad){
    ctx.beginPath(); ctx.moveTo(pad,y0); ctx.lineTo(W-pad,y0); ctx.stroke();
    ctx.fillStyle='#15803d'; ctx.font='10px system-ui';
    ctx.fillText('converged',W-pad-56,y0-2);
  }
  ctx.setLineDash([]);

  // Main curve
  ctx.strokeStyle='#1d6fb8'; ctx.lineWidth=2;
  ctx.beginPath();
  history.forEach((v,i)=>{
    const lv = Math.log10(Math.max(v,1e-15));
    const x = pad + i*xScale;
    const y = H - pad - (lv - minLog)*yScale;
    i===0 ? ctx.moveTo(x,y) : ctx.lineTo(x,y);
  });
  ctx.stroke();
}

function step(){
  if(iter >= 80){ statusEl.textContent = fmt('Converged after %iter% iterations.', { iter: iter }); return; }
  const {Q,R} = qrDecomp(A);
  A = matMul(R, Q);  // A_{k+1} = R * Q  (similarity transform, preserves eigenvalues)
  iter++;
  history.push(offDiagNorm(A));
  iterEl.textContent = iter;
  renderMatrix(); renderDiag(); renderChart();
  const od = offDiagNorm(A);
  if(od < 1e-6){
    statusEl.textContent = fmt('Converged! Off-diagonal norm = %toExponential% after %iter% steps.', { toExponential: od.toExponential(2), iter: iter });
    stopRun();
  } else {
    statusEl.textContent = fmt('Off-diagonal norm: %toExponential%', { toExponential: od.toExponential(3) });
  }
}

function startRun(){
  if(running) return;
  running = true;
  document.getElementById('runBtn').textContent = 'Pause';
  runTimer = setInterval(()=>{ step(); if(offDiagNorm(A)<1e-6) stopRun(); }, 180);
}
function stopRun(){
  running = false; clearInterval(runTimer);
  document.getElementById('runBtn').textContent = 'Run';
}
function reset(){
  stopRun();
  A = [...A0]; iter=0; history=[offDiagNorm(A)];
  iterEl.textContent=0; statusEl.textContent='';
  renderMatrix(); renderDiag(); renderChart();
}

document.getElementById('stepBtn').addEventListener('click', step);
document.getElementById('runBtn').addEventListener('click', ()=>{ running ? stopRun() : startRun(); });
document.getElementById('resetBtn').addEventListener('click', reset);

trueDiv.textContent = TRUE_EV.map(v=>v.toFixed(4)).join(',  ');
renderMatrix(); renderDiag(); renderChart();

console.log('True eigenvalues (bisection):', TRUE_EV.map(v=>v.toFixed(6)));
Result
Console

Press Step to advance one iteration, or Run to animate. Watch the off-diagonal entries shrink toward zero โ€” when they are small enough, the numbers on the main diagonal are the eigenvalues. The true eigenvalues (computed exactly) are shown for comparison.

Convergence and Cost

Status: solved (1961, Francis & Kublanovskaya) โ€” efficient and numerically stable in practice.

The QR algorithm is beautifully analyzable:

  • Why it works. Each QR step is an orthogonal similarity transformation: Ak+1=RkQk=QkTAkQkA_{k+1} = R_k Q_k = Q_k^T A_k Q_k. Similarity transformations preserve eigenvalues, so every iterate has the same eigenvalues. As the iterates accumulate, the off-diagonal entries are driven toward zero by the geometry of repeated matrix multiplication.
  • Convergence rate. Without modifications, convergence can be slow โ€” linear at best. Add a shift (subtract a scalar before factoring, add it back after) and convergence becomes cubic for simple eigenvalues. The Wilkinson shift chooses the shift adaptively, giving outstanding reliability in practice.
  • Cost per step. A naive QR step on an nร—n matrix costs O(n3)O(n^{3}) arithmetic operations. The practical algorithm first reduces A to upper Hessenberg form (O(n3)O(n^{3}) once), then each QR step costs only O(n2)O(n^{2}). Total cost for all n eigenvalues: O(n3)O(n^{3}), which is optimal to within constants for dense matrices.
  • Symmetric vs general. For symmetric (or Hermitian) matrices, the algorithm reduces to the symmetric QR algorithm with guaranteed real eigenvalues and globally quadratic convergence. For general matrices, convergence is guaranteed under mild conditions, but pathological cases can arise โ€” the Francis double-step and implicit-shift variants handle these in practice.

The QR algorithm is considered solved in the sense that the standard implementations (LAPACK's dsyev, dgeev) are numerically stable, efficient, and trusted in production software worldwide. The theoretical question of an optimal O(n2)O(n^{2}) eigenvalue algorithm remains open, connecting to deep questions in matrix multiplication.

Where It Matters

Eigenvalue computation is one of the most broadly useful operations in all of applied mathematics. The QR algorithm sits at the center of:

  • Structural engineering: the natural frequencies of a building or bridge are eigenvalues of its stiffness matrix. Compute them wrong and resonance can be catastrophic.
  • Quantum physics: the energy levels of an atom or molecule are eigenvalues of its Hamiltonian operator, discretized into a matrix and solved by QR-based methods.
  • Machine learning (PCA): principal component analysis finds the directions of greatest variance by computing eigenvalues (and eigenvectors) of the data covariance matrix โ€” the QR algorithm under the hood.
  • Graph analysis and PageRank: the ranking score in PageRank is the dominant eigenvector of a huge sparse matrix; power iteration (a cousin of QR) drives it.
  • Control systems: the stability of a controller depends on whether all eigenvalues of a system matrix lie inside the unit disk or left half-plane โ€” a question answered by the QR algorithm.
  • Scientific simulation: differential equations, fluid dynamics, and population models are all studied via the eigenstructure of linearized operators.

Without fast, reliable eigenvalue computation, vast swaths of modern science and engineering would grind to a halt.

Conclusion

The QR algorithm is a rare thing in computing: an idea that is both simple to state and provably optimal in practice. Factor the matrix into Q and R, multiply in reverse order, repeat โ€” and the eigenvalues crystallize on the diagonal as if by magic.

The magic is not accidental. Each step is a similarity transformation, so eigenvalues are preserved while the off-diagonal entries are steadily squeezed. Add shifts for speed, Hessenberg reduction to cut cost, and implicit steps for numerical stability, and you have the algorithm that powers LAPACK, MATLAB, NumPy, and virtually every piece of scientific software on earth.

If you have ever computed a principal component, simulated a vibrating structure, or trusted a control system, you have quietly relied on Francis and Kublanovskaya's 1961 discovery โ€” one of the most consequential loops in the history of algorithms.

Share this article

Pick a channel โ€” or use your device's native share sheet.

Comments

Loading comments...

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