Householder QR

Introduction

Somewhere deep inside your spreadsheet solver, your least-squares fit, and your eigenvalue routine lives a single workhorse: QR factorization. Every matrix $A$ can be written as $A = QR$ , where $Q$ is orthogonal ( $Q^T Q = I$ ) and $R$ is upper triangular. That decomposition turns messy simultaneous equations into a clean back-substitution.

The challenge is computing $Q$ and $R$ without accumulating the floating-point errors that plague naive approaches like Gram–Schmidt. The answer, developed by Alston Scott Householder in 1958, is disarmingly geometric: take the first column of $A$ and reflect it onto the first coordinate axis, killing every entry below the diagonal in one shot. Repeat column by column and you are done.

Each reflection is an orthogonal transformation — it preserves lengths and angles, so it cannot amplify rounding errors. That one insight made Householder QR the algorithm of choice for dense matrices in numerical software for the next seven decades.

Watch the Reflections

The demo below shows a random $4 \times 4$ matrix. At each step one Householder reflection zeroes out everything below the diagonal in the current column. Click Next step to apply each reflection and watch the matrix march toward upper-triangular form.

<!-- Householder QR demo: matrix display + step button -->
<p class="hint">Each click applies one Householder reflection, zeroing out the column below the diagonal. After N-1 steps the matrix is upper-triangular (R).</p>
<div class="matrix-area">
  <div class="mat-label">Matrix A (current state)</div>
  <table id="matA" class="mat"></table>
</div>
<div class="status" id="status">Press "Next step" to apply the first reflection.</div>
<div class="btns">
  <button id="btnStep" type="button">Next step</button>
  <button id="btnReset" type="button" class="ghost">New matrix</button>
</div>
<div class="step-info" id="stepInfo"></div>

/* Styles for the Householder QR demo */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; font-size: 14px; }
.hint { font-size: .9rem; color: #444; margin: 0 0 .7rem; line-height: 1.45; }
.matrix-area { margin: .5rem 0; }
.mat-label { font-size: .8rem; font-weight: 700; color: #1d3557; margin-bottom: .2rem; letter-spacing: .04em; }
.mat { border-collapse: separate; border-spacing: 3px; }
.mat td {
  width: 54px; height: 36px; text-align: right; padding: 0 6px;
  font: 600 12px ui-monospace, monospace; border-radius: 5px;
  transition: background .35s, color .35s;
}
.cell-normal { background: #e8eef3; color: #1d3557; border: 1px solid #cdd9e3; }
.cell-zero   { background: #d4edda; color: #0a7d33; border: 1px solid #b8dbc3; }
.cell-pivot  { background: #fff3cd; color: #856404; border: 1px solid #ffe08a; }
.cell-active { background: #f8d7da; color: #842029; border: 1px solid #f5c2c7; }
.status { font-size: 1rem; font-weight: 600; margin: .5rem 0; min-height: 1.4em; }
.status.ok   { color: #0a7d33; }
.status.info { color: #1d3557; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 14px system-ui; padding: .45rem .9rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .4; cursor: not-allowed; }
.step-info { font-size: .82rem; color: #555; min-height: 1.2em; line-height: 1.5; }

// Householder QR factorization — step-by-step demo
// At each step k we build a Householder reflector H_k that zeroes out A[k+1..n, k].
const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));

// Matrix dimension
const N = 4;

// Generate a random matrix with entries in [-5, 5]
function randMatrix() {
  const A = [];
  for (let i = 0; i < N; i++) {
    A.push([]);
    for (let j = 0; j < N; j++) {
      A[i].push(Math.round((Math.random() * 10 - 5) * 10) / 10);
    }
  }
  return A;
}

// Deep-copy a matrix
function copyM(A) { return A.map(r => r.slice()); }

// Euclidean norm of a vector
function norm2(v) { return Math.sqrt(v.reduce((s, x) => s + x * x, 0)); }

// Apply one Householder step to column k of matrix A
// Returns the updated matrix with zeros below the diagonal in column k.
function householderStep(A, k) {
  const m = A.length;
  // Extract the subcolumn from row k downward
  const x = [];
  for (let i = k; i < m; i++) x.push(A[i][k]);
  const nrm = norm2(x);
  if (nrm < 1e-12) return copyM(A); // column already zero, nothing to do
  // Build the Householder vector u = x + sign(x[0])*||x||*e1
  const u = x.slice();
  u[0] += (x[0] >= 0 ? 1 : -1) * nrm;
  const nu = norm2(u);
  if (nu < 1e-15) return copyM(A);
  const un = u.map(v => v / nu);
  // Apply H = I - 2*u*u^T to the submatrix from column k onward
  const B = copyM(A);
  for (let j = k; j < m; j++) {
    // dot product of u with each column of the submatrix
    let dot = 0;
    for (let i = 0; i < x.length; i++) dot += un[i] * B[i + k][j];
    for (let i = 0; i < x.length; i++) B[i + k][j] -= 2 * un[i] * dot;
  }
  // Force exact zeros below diagonal to avoid floating-point drift
  for (let i = k + 1; i < m; i++) B[i][k] = 0;
  return B;
}

let origA, curA, step;

function init() {
  origA = randMatrix();
  curA = copyM(origA);
  step = 0;
  document.getElementById('btnStep').disabled = false;
  render(curA, -1);
  setStatus('Press "Next step" to apply the first reflection.', 'info');
  document.getElementById('stepInfo').textContent = '';
}

// Format a number for display: treat tiny values as zero
function fmtNum(v) {
  if (Math.abs(v) < 5e-10) return '0';
  return v.toFixed(2);
}

function render(A, activeCol) {
  const tbl = document.getElementById('matA');
  tbl.innerHTML = '';
  for (let i = 0; i < N; i++) {
    const tr = document.createElement('tr');
    for (let j = 0; j < N; j++) {
      const td = document.createElement('td');
      td.textContent = fmtNum(A[i][j]);
      if (j === activeCol && i > activeCol) {
        td.className = 'cell-active'; // entries being zeroed this step
      } else if (j < activeCol && i > j) {
        td.className = 'cell-zero'; // already-zeroed entries from previous steps
      } else if (i === j && j === activeCol) {
        td.className = 'cell-pivot'; // current pivot diagonal entry
      } else {
        td.className = 'cell-normal';
      }
      tr.appendChild(td);
    }
    tbl.appendChild(tr);
  }
}

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

document.getElementById('btnStep').addEventListener('click', () => {
  if (step >= N - 1) return;
  const next = householderStep(curA, step);
  curA = next;
  render(curA, step);
  step++;
  if (step < N - 1) {
    setStatus(i18nFmt('Step %k% complete. Column %k% is now upper-triangular.', { k: step }), 'info');
    document.getElementById('stepInfo').textContent =
      i18nFmt('Next: reflection will zero out column %k%.', { k: step + 1 });
  } else {
    setStatus('Done! The matrix is now upper-triangular (R).', 'ok');
    document.getElementById('btnStep').disabled = true;
    document.getElementById('stepInfo').textContent = 'All reflections applied. The accumulated reflectors give Q^T.';
    render(curA, -1);
  }
});

document.getElementById('btnReset').addEventListener('click', init);

init();

<!-- Householder QR demo: matrix display + step button -->
<p class="hint">Each click applies one Householder reflection, zeroing out the column below the diagonal. After N-1 steps the matrix is upper-triangular (R).</p>
<div class="matrix-area">
  <div class="mat-label">Matrix A (current state)</div>
  <table id="matA" class="mat"></table>
</div>
<div class="status" id="status">Press "Next step" to apply the first reflection.</div>
<div class="btns">
  <button id="btnStep" type="button">Next step</button>
  <button id="btnReset" type="button" class="ghost">New matrix</button>
</div>
<div class="step-info" id="stepInfo"></div>

/* Styles for the Householder QR demo */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; font-size: 14px; }
.hint { font-size: .9rem; color: #444; margin: 0 0 .7rem; line-height: 1.45; }
.matrix-area { margin: .5rem 0; }
.mat-label { font-size: .8rem; font-weight: 700; color: #1d3557; margin-bottom: .2rem; letter-spacing: .04em; }
.mat { border-collapse: separate; border-spacing: 3px; }
.mat td {
  width: 54px; height: 36px; text-align: right; padding: 0 6px;
  font: 600 12px ui-monospace, monospace; border-radius: 5px;
  transition: background .35s, color .35s;
}
.cell-normal { background: #e8eef3; color: #1d3557; border: 1px solid #cdd9e3; }
.cell-zero   { background: #d4edda; color: #0a7d33; border: 1px solid #b8dbc3; }
.cell-pivot  { background: #fff3cd; color: #856404; border: 1px solid #ffe08a; }
.cell-active { background: #f8d7da; color: #842029; border: 1px solid #f5c2c7; }
.status { font-size: 1rem; font-weight: 600; margin: .5rem 0; min-height: 1.4em; }
.status.ok   { color: #0a7d33; }
.status.info { color: #1d3557; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 14px system-ui; padding: .45rem .9rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .4; cursor: not-allowed; }
.step-info { font-size: .82rem; color: #555; min-height: 1.2em; line-height: 1.5; }

// Householder QR factorization — step-by-step demo
// At each step k we build a Householder reflector H_k that zeroes out A[k+1..n, k].
const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));

// Matrix dimension
const N = 4;

// Generate a random matrix with entries in [-5, 5]
function randMatrix() {
  const A = [];
  for (let i = 0; i < N; i++) {
    A.push([]);
    for (let j = 0; j < N; j++) {
      A[i].push(Math.round((Math.random() * 10 - 5) * 10) / 10);
    }
  }
  return A;
}

// Deep-copy a matrix
function copyM(A) { return A.map(r => r.slice()); }

// Euclidean norm of a vector
function norm2(v) { return Math.sqrt(v.reduce((s, x) => s + x * x, 0)); }

// Apply one Householder step to column k of matrix A
// Returns the updated matrix with zeros below the diagonal in column k.
function householderStep(A, k) {
  const m = A.length;
  // Extract the subcolumn from row k downward
  const x = [];
  for (let i = k; i < m; i++) x.push(A[i][k]);
  const nrm = norm2(x);
  if (nrm < 1e-12) return copyM(A); // column already zero, nothing to do
  // Build the Householder vector u = x + sign(x[0])*||x||*e1
  const u = x.slice();
  u[0] += (x[0] >= 0 ? 1 : -1) * nrm;
  const nu = norm2(u);
  if (nu < 1e-15) return copyM(A);
  const un = u.map(v => v / nu);
  // Apply H = I - 2*u*u^T to the submatrix from column k onward
  const B = copyM(A);
  for (let j = k; j < m; j++) {
    // dot product of u with each column of the submatrix
    let dot = 0;
    for (let i = 0; i < x.length; i++) dot += un[i] * B[i + k][j];
    for (let i = 0; i < x.length; i++) B[i + k][j] -= 2 * un[i] * dot;
  }
  // Force exact zeros below diagonal to avoid floating-point drift
  for (let i = k + 1; i < m; i++) B[i][k] = 0;
  return B;
}

let origA, curA, step;

function init() {
  origA = randMatrix();
  curA = copyM(origA);
  step = 0;
  document.getElementById('btnStep').disabled = false;
  render(curA, -1);
  setStatus('Press "Next step" to apply the first reflection.', 'info');
  document.getElementById('stepInfo').textContent = '';
}

// Format a number for display: treat tiny values as zero
function fmtNum(v) {
  if (Math.abs(v) < 5e-10) return '0';
  return v.toFixed(2);
}

function render(A, activeCol) {
  const tbl = document.getElementById('matA');
  tbl.innerHTML = '';
  for (let i = 0; i < N; i++) {
    const tr = document.createElement('tr');
    for (let j = 0; j < N; j++) {
      const td = document.createElement('td');
      td.textContent = fmtNum(A[i][j]);
      if (j === activeCol && i > activeCol) {
        td.className = 'cell-active'; // entries being zeroed this step
      } else if (j < activeCol && i > j) {
        td.className = 'cell-zero'; // already-zeroed entries from previous steps
      } else if (i === j && j === activeCol) {
        td.className = 'cell-pivot'; // current pivot diagonal entry
      } else {
        td.className = 'cell-normal';
      }
      tr.appendChild(td);
    }
    tbl.appendChild(tr);
  }
}

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

document.getElementById('btnStep').addEventListener('click', () => {
  if (step >= N - 1) return;
  const next = householderStep(curA, step);
  curA = next;
  render(curA, step);
  step++;
  if (step < N - 1) {
    setStatus(i18nFmt('Step %k% complete. Column %k% is now upper-triangular.', { k: step }), 'info');
    document.getElementById('stepInfo').textContent =
      i18nFmt('Next: reflection will zero out column %k%.', { k: step + 1 });
  } else {
    setStatus('Done! The matrix is now upper-triangular (R).', 'ok');
    document.getElementById('btnStep').disabled = true;
    document.getElementById('stepInfo').textContent = 'All reflections applied. The accumulated reflectors give Q^T.';
    render(curA, -1);
  }
});

document.getElementById('btnReset').addEventListener('click', init);

init();

Result

Console

Notice that after step $k$ the entries below the diagonal in column $k$ are exactly zero (up to floating-point rounding). The matrix never needs to store $Q$ explicitly — the reflection vectors are enough to reconstruct it later. By the final step the matrix has become $R$ , and the accumulated reflections give $Q^T$ .

Cost and Stability

How expensive is Householder QR, and how much should you trust the answer?

Cost. Factoring an $m \times n$ matrix ( $m \ge n$ ) requires roughly $2n^2(m - n/3)$ floating-point operations — $O(n^3)$ for square matrices. That is the same order as Gaussian elimination, but with a constant factor about twice as large. The extra work buys something valuable.
Backward stability. Householder reflections are orthogonal transformations: they rotate and reflect but never stretch. A theorem due to Wilkinson shows that the computed $\hat{Q}$ and $\hat{R}$ satisfy $\hat{Q}\hat{R} = A + \delta A$ where $\|\delta A\| / \|A\| = O(\varepsilon_{\text{machine}})$ . In plain terms: the result is the exact QR of a matrix that differs from $A$ by a tiny rounding-error perturbation. No other QR method offers this guarantee as cleanly.
Gram–Schmidt comparison. Classical Gram–Schmidt is cheaper per column but can lose orthogonality catastrophically when columns are nearly dependent. Modified Gram–Schmidt is better but still only conditionally stable. Householder wins on stability for any well-posed problem.
LAPACK's choice. The routine dgeqrf in LAPACK — the reference for nearly all scientific computing — implements Householder QR with blocked updates for cache efficiency.

The algorithm is solved: it has the optimal stability class for this problem, and its cost is within a constant of the information-theoretic minimum.

Where It Matters

QR factorization by Householder reflections is one of the most-called subroutines in all of scientific software:

Least-squares regression. To fit a model $Ax \approx b$ in the least-squares sense, form the QR factorization of $A$ and solve $Rx = Q^T b$ by back-substitution. This is numerically far superior to forming the normal equations $A^T A x = A^T b$ , which square the condition number.
Eigenvalue algorithms. The QR algorithm — the standard method for computing all eigenvalues of a dense matrix — applies Householder QR at every iteration. Without a stable QR step the iteration would diverge in finite precision.
Signal and image processing. Orthogonal projections computed via QR appear in beamforming, noise removal, and the discrete cosine transform used in JPEG compression.
Compressed sensing and sparse recovery. Many algorithms for recovering sparse signals (LASSO, OMP) rely on QR to maintain numerical orthogonality as columns are selected greedily.
Quantum chemistry. Self-consistent field methods (Hartree–Fock, DFT) iterate on large dense matrices; QR stabilizes each step.

Householder's 1958 paper solved the stability problem so thoroughly that its approach has barely changed. The blocked variant in LAPACK is a tuning detail; the mathematical core is exactly what Householder wrote.

Conclusion

Householder QR is a rare case in numerical computing where a single geometric idea — reflect a vector onto an axis — solved a problem so completely that the field has not had to revisit it. Chain $n$ such reflections, accumulate them implicitly, and you have an orthogonal $Q$ and upper-triangular $R$ whose product is your original matrix, accurate to machine precision.

The next time you call a least-squares solver, run a dimensionality reduction, or compute eigenvalues, a Householder reflection is almost certainly doing the heavy lifting — quietly, stably, and exactly as Alston Scott Householder intended in 1958.

Comments

Loading comments...

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

Cookies & privacy

Reflect columns onto axes, zero out the mess