Householder QR Factorization

A complete reference for the Householder QR factorization algorithm (Golub & Van Loan, Matrix Computations, Algorithm 5.2.1). Five parts covering the same algorithm from different angles — definitions, invariants, a question-and-answer walkthrough, a traditional narrative, and a concept map. They are not meant to be read in sequence; they are meant to be used as a reference you return to.

How to use this document

Part 2 (Definitions) and Part 3 (Invariants) are the core: every claim in the document is grounded there. Part 4 (Dialogue) is the recommended first read for a new reader — each question is an independent entry point. Part 5 (Narrative) is the synthesis, with a worked example and code. Part 6 (Concept Map) shows the dependency graph of ideas — which concepts are load-bearing, and which connections are non-obvious.

Part 1 — Context

The Householder QR factorization is Algorithm 5.2.1 in Golub and Van Loan, Matrix Computations (4th ed.). It is the standard method for computing the QR factorization of a dense matrix. The storage scheme it uses — packing Householder vectors into the zeros they create — is the format returned by the LAPACK routine DGEQRF. Anyone reading that output without knowing the storage scheme will be confused by the lower triangle.

Part 2 — Definitions

Every term used in Parts 3–5 is defined here. No term in the invariants requires knowledge from outside this section.

Matrix.

A rectangular array of numbers with $m$ rows and $n$ columns. Written $A$, size $m\times n$. Entry in row $i$, column $j$ is written $A_{ij}$.

Vector.

A column of numbers. Written $x$. The $i$-th entry is $x_i$. A vector of length $m$ is an $m\times 1$ matrix.

Length of a vector.

$\Vert x\Vert =\sqrt{x_1^2+x_2^2+\cdots +x_m^2}$. Also called the Euclidean norm.

Dot product.

$v^{T}x=v_1{x}_1+v_2{x}_2+\cdots +v_m{x}_m$. Always a single number.

Outer product.

$v{v}^T$ is an $m\times m$ matrix whose $(i,j)$ entry is $v_i{v}_j$.

Identity matrix.

$I$ is the square matrix with 1s on the diagonal and 0s elsewhere. $Ix=x$ for any vector $x$.

Transpose.

$A^T$ has rows and columns of $A$ swapped: $(A^T{)}_{ij}=A_{ji}$.

Orthogonal matrix.

A square matrix $Q$ satisfying $Q^{T}Q=I$, equivalently $Q^{-1}=Q^T$. Applying $Q$ to any vector preserves its length: $\Vert Qx\Vert =\Vert x\Vert$ for all $x$.

Upper triangular matrix.

A matrix $R$ where $R_{ij}=0$ whenever $i>j$. Everything below the main diagonal is zero.

$e_1$.

The vector $[1,0,0,\dots ,0{]}^T$. Points along the first coordinate axis.

$\mathrm{sign}(x_1)$.

Returns $+1$ if $x_1\ge 0$, returns $-1$ if $x_1<0$.

Submatrix $A(j:m,j:n)$.

The block of $A$ formed by rows $j$ through $m$ and columns $j$ through $n$.

Rank-1 update.

An operation of the form $B\leftarrow B-\tau v(v^{T}B)$, where $v$ is a vector and $\tau$ is a scalar. It modifies every column of $B$ by subtracting a scalar multiple of $v$. Cost: $O(mn)$.

Part 3 — Glossary of Invariants

Every statement below is necessarily true. Each uses only terms from Part 2 or from earlier invariants. Statements are ordered by dependency — no invariant uses a term introduced later.

The object

1. Every $m\times n$ matrix $A$ with $m\ge n$ can be written exactly as $A=QR$, where $Q$ is an $m\times m$ orthogonal matrix and $R$ is an $m\times n$ upper triangular matrix. This is not an approximation.
2. $R$ being upper triangular means every entry below the main diagonal is zero. The diagonal and above are generally nonzero.
3. $Q$ being orthogonal means $Q^{T}Q=I$: applying $Q$ changes the direction of vectors without changing their length.
4. Because $Q^{-1}=Q^T$, the operation $Q^{T}b$ — reorienting a vector $b$ — costs $O(mn)$ and destroys no information.
5. The factorization $A=QR$ always exists for any $m\times n$ matrix with $m\ge n$.
6. The factorization is not unique: negating the $j$-th column of $Q$ and the $j$-th row of $R$ simultaneously gives a different but equally valid factorization. The signs of the diagonal of $R$ are free.

Why it is useful

7. The system $Ax=b$ is equivalent to $Rx=Q^{T}b$, obtained by substituting $A=QR$ and multiplying both sides by $Q^T$.
8. Because $R$ is upper triangular, $Rx=Q^{T}b$ is solved from the bottom row upward, one unknown at a time. This is back-substitution; it costs $O(n^2)$ operations.
9. The least squares problem — find $x$ minimizing $\Vert Ax-b\Vert$ — has the solution obtained by: computing $Q^{T}b$, taking its first $n$ entries, then solving $Rx=(Q^{T}b{)}_{1:n}$ by back-substitution. No other steps are required.
10. Applying $Q^T$ to $b$ does not change $\Vert b\Vert$, because $Q^T$ is orthogonal. It reorients $b$ without distorting it.

The key constraint

11. Any transformation applied to $A$ during the factorization must be orthogonal — it must preserve lengths. Setting entries to zero directly changes the matrix itself, not just its coordinate representation, and is not permitted.
12. Zeroing all entries below position $j$ in a column $x$ is equivalent to finding an orthogonal transformation $H$ such that $Hx=[\pm \Vert x\Vert ,0,0,\dots ,0{]}^T$. The first entry must be $\pm \Vert x\Vert$ because $H$ preserves length and the right-hand side must have the same norm as $x$.

The Householder reflector

13. For any nonzero vector $v$, the scalar $\tau =2/(v^{T}v)$ satisfies $0<\tau \le 2$.
14. The matrix $H=I-\tau v{v}^T$, where $\tau =2/(v^{T}v)$, is called a Householder reflector.
15. $H$ is orthogonal: $H^{T}H=I$. Applying $H$ to any vector preserves its length.
16. $H$ is symmetric: $H^T=H$.
17. $H$ is its own inverse: $H^2=I$. Applying $H$ twice returns any vector exactly to where it started.
18. $H$ reflects space across the hyperplane $\{w:v^{T}w=0\}$ — the set of all vectors perpendicular to $v$. The component of any vector along $v$ is negated; components perpendicular to $v$ are unchanged.
19. $H$ is never stored or formed as an explicit $m\times m$ matrix. Applying $H$ to a vector $b$ uses: $Hb=b-\tau v(v^{T}b)$. This requires one dot product, one scaling, and one vector subtraction. Cost: $O(m)$.
20. Applying $H$ to an $m\times n$ matrix $B$ uses: $HB=B-\tau v(v^{T}B)$, where $v^{T}B$ is a row vector of $n$ numbers. Cost: $O(mn)$.

Constructing $v$ to zero a column

21. Given a vector $x$ of length $m$, define $v=x+\mathrm{sign}(x_1)\Vert x\Vert e_1$. The reflector $H=I-\tau v{v}^T$ with $\tau =2/(v^{T}v)$ then satisfies $Hx=-\mathrm{sign}(x_1)\Vert x\Vert e_1$: the first entry is $-\mathrm{sign}(x_1)\Vert x\Vert$ and all other entries are zero.
22. The vector $v$ from invariant 21 points from the target $-\mathrm{sign}(x_1)\Vert x\Vert e_1$ to the current position $x$. The hyperplane perpendicular to $v$ bisects $x$ and its target. Reflecting $x$ across this bisector maps it exactly onto the target.
23. The zeros in $Hx$ are not set directly. They appear because $x$ has been rotated to point along $e_1$. A vector pointing along $e_1$ has zero in every position except the first by definition.
24. The sign $\mathrm{sign}(x_1)$ in invariant 21 is chosen so that $v_1=x_1+\mathrm{sign}(x_1)\Vert x\Vert$ is always large in magnitude. The alternative sign gives $v_1=x_1-\Vert x\Vert$, which is nearly zero when $x_1\approx \Vert x\Vert$. Near-zero $v_1$ causes catastrophic cancellation when $\tau =2/(v^{T}v)$ is computed, amplifying floating point errors by factors of millions. Example: $x=[0.9999,0.01,0.01{]}^T$, $\Vert x\Vert \approx 1.0001$. Wrong sign: $v_1=0.9999-1.0001=-0.0002$ (dangerous). Correct sign: $v_1=0.9999+1.0001=2.0000$ (safe).

The algorithm

25. The algorithm performs $n$ steps. At step $j$ (for $j=1,2,\dots ,n$), it operates on the submatrix $A(j:m,j:n)$.
26. At step $j$: extract $x=A(j:m,j)$. Construct $v_j$ from $x$ using invariant 21. The reflector $H_j=I-{\tau }_j{v}_j{v}_j^T$ maps $x$ to $-\mathrm{sign}(x_1)\Vert x\Vert e_1$.
27. Apply $H_j$ to the entire submatrix $A(j:m,j:n)$ — not just column $j$. Applying it to one column but not the others would leave the matrix inconsistent.
28. $H_j$ acts only on rows $j$ through $m$. It therefore cannot disturb columns $1$ through $j-1$, which are already in their final form. This is what makes the column-by-column approach correct, not merely plausible.
29. After applying $H_j$: entry $A(j,j)$ becomes $r_{jj}=-\mathrm{sign}(x_1)\Vert x\Vert$. Entries $A(j+1:m,j)$ become exactly zero.
30. After all $n$ steps: the upper triangle of the working matrix (including diagonal) contains $R$. The product $Q=H_1{H}_2\cdots H_n$ is the orthogonal factor. $Q$ is never computed explicitly.
31. The total cost of the factorization is $2m{n}^2-\frac{2}{3}{n}^3$ floating point operations. For a square $m=n$ matrix this is $\frac{4}{3}{n}^3$ — twice the cost of LU factorization ($\frac{2}{3}{n}^3$). QR is preferred not for square systems but for least squares and rank-deficient problems, where its orthogonality guarantees are essential.

The storage scheme

32. After step $j$, entries $A(j+1:m,j)$ are zero. These slots held no information about the original $A$ and will not be needed going forward. They are free to reuse.
33. By normalizing $v_j$ so that $v_j(1)=1$ (dividing $v_j$ by its first component after construction), the first entry is always exactly 1. Only $v_j(2),\dots ,v_j(m-j+1)$ — a total of $m-j$ entries — need to be stored.
34. The number of free slots created at step $j$ (from invariant 32) is $m-j$. The number of entries of $v_j$ that need storing (from invariant 33) is also $m-j$. These counts are equal. This exact coincidence is what enables the storage trick.
35. The free slots in the lower triangle of column $j$ store $v_j(2)$ through $v_j(m-j+1)$. The diagonal entry stores $r_{jj}$. The upper triangle stores $R$.
36. One scalar ${\tau }_j=2/(v_j^T{v}_j)$ is stored per step in a separate length-$n$ vector $\tau$. It cannot be recovered from the stored $v_j$ entries alone without recomputation.
37. The output is: (a) the overwritten array, upper triangle = $R$, strictly lower triangle = essential parts of $v_1,\dots ,v_n$; and (b) the vector $\tau$ of length $n$. Together these are a complete lossless encoding of both $Q$ and $R$.
38. The lower triangle of the output is not part of $R$, not zero, and not padding. It is the compressed representation of $Q$. This is the LAPACK DGEQRF format.

Recovering and applying $Q$

39. To apply $Q^T$ to a vector $b$: for $j=1,2,\dots ,n$ (forward), reconstruct $v_j$ by prepending the implied 1 to the stored subdiagonal of column $j$, then compute $b(j:m)\leftarrow b(j:m)-{\tau }_j\cdot v_j\cdot (v_j^{T}b(j:m))$.
40. To apply $Q$ to a vector $b$: run the same loop in reverse, $j=n,n-1,\dots ,1$. This is correct because $(H_1{H}_2\cdots H_n{)}^T=H_n^T\cdots H_1^T=H_n\cdots H_1$, since each $H_j$ is symmetric by invariant 16.
41. Forming $Q$ as an explicit $m\times m$ matrix costs $O(m^{2}n)$. Applying $Q$ or $Q^T$ to a single vector using the stored representation costs $O(mn)$. Explicit $Q$ is almost never needed.

What the algorithm is not doing

42. It is not setting entries to zero directly. Zeros appear because columns are rotated onto coordinate axes.
43. It is not building $Q$ first and then computing $Q^{T}A$. It applies each reflector directly to the working matrix as it goes.
44. It is not storing Householder vectors in a separate data structure. They live in the zeros they created.
45. The lower triangle of the output is not a mistake. It is the compressed recipe for $Q$.

Part 4 — Dialogue

Each question is an independent entry point. Answers are self-contained.

for j = 1 to n:
    reconstruct v_j: prepend 1 to stored subdiagonal of column j
    b[j:] ← b[j:] − τⱼ · vⱼ · (vⱼᵀ b[j:])

Part 5 — Traditional Narrative

What QR factorization is

Given an $m\times n$ matrix $A$ with $m\ge n$, QR factorization produces $Q$ and $R$ such that $A=QR$. $Q$ is $m\times m$ orthogonal ($Q^{T}Q=I$); $R$ is $m\times n$ upper triangular. The factorization always exists and is unique up to sign choices on the diagonal of $R$.

The system $Ax=b$ becomes $Rx=Q^{T}b$, solved by back-substitution. The least squares minimizer of $\Vert Ax-b\Vert$ satisfies the same triangular system. In both cases, $Q^T$ is applied without forming $Q$ explicitly.

The Householder reflector

For any nonzero vector $v$, the matrix

is symmetric and orthogonal. It reflects space across the hyperplane $\{w:v^{T}w=0\}$. It is its own inverse: $H^2=I$.

The key property: given any vector $x$, choose

Then $Hx=-\mathrm{sign}(x_1)\Vert x\Vert e_1$: first entry $-\mathrm{sign}(x_1)\Vert x\Vert$, all other entries zero.

The sign convention ensures $v_1=x_1+\mathrm{sign}(x_1)\Vert x\Vert$ is large, avoiding catastrophic cancellation. The wrong sign gives $v_1\approx 0$ when $x$ nearly points along $e_1$, causing $\tau =2/(v^{T}v)$ to explode.

$H$ is never formed explicitly. Applied to matrix $B$:

Cost: $O(mn)$ rather than $O(m^{2}n)$.

The algorithm

At step $j$, extract $x=A(j:m,j)$. Construct $v_j$ and ${\tau }_j$. Apply:

After this update: $A(j,j)=r_{jj}$, entries $A(j+1:m,j)=0$. Earlier columns are untouched because $H_j$ acts only on rows $j$ through $m$.

After all $n$ steps: $A$ has been overwritten with $R$, and $Q=H_1{H}_2\cdots H_n$.

The compact storage scheme

After step $j$, the entries $A(j+1:m,j)$ are zero and free. Normalizing so $v_j(1)=1$, only $v_j(2:end)$ needs storing — exactly $m-j$ entries, matching the $m-j$ free slots. The diagonal stores $r_{jj}$.

The output contains two things simultaneously: upper triangle = $R$, strictly lower triangle = essential Householder vectors. A separate length-$n$ vector $\tau$ completes the encoding. This is the LAPACK DGEQRF format.

Applying $Q$ without forming it

To apply $Q^T$ to $b$:

To apply $Q$: same loop, reversed. Cost: $O(mn)$. Forming $Q$ explicitly costs $O(m^{2}n)$ and is almost never needed.

Complete numeric example

Step 1. $x=[1,1,0{]}^T$, $\Vert x\Vert =\sqrt{2}$.

Apply $H_1$ to all three columns via $A\leftarrow A-{\tau }_1{v}_1(v_1^{T}A)$.

$v_1^{T}A$: compute row by row.

• Column 1: $v_1^{T}x=(1+\sqrt{2})(1)+(1)(1)+0=2+\sqrt{2}$
• Column 2: $(1+\sqrt{2})(1)+(1)(0)+0=1+\sqrt{2}$
• Column 3: $(1+\sqrt{2})(0)+(1)(1)+0=1$

So ${\tau }_1(v_1^{T}A)=\frac{1}{2+\sqrt{2}}[2+\sqrt{2},1+\sqrt{2},1]$.

After update:

where $v_1(2)=1/(1+\sqrt{2})$ and $v_1(3)=0$ are packed into the lower triangle, and $r_{11}=-\sqrt{2}$.

Step 2. Subvector: $x=[-\frac{1}{2},1{]}^T$, $\Vert x\Vert =\frac{\sqrt{5}}{2}$.

Since $x_1=-\frac{1}{2}<0$, $\mathrm{sign}(x_1)=-1$:

Apply $H_2$ to the $2\times 2$ subblock. After this step: $A(3,2)=0$, $r_{22}=\frac{\sqrt{5}}{2}$. $v_2(2)=1$ packed into position $(3,2)$.

Step 3. Single entry: $A(3,3)$ is already in upper triangular position. $r_{33}$ is read directly. No reflection needed.

Final packed array:

Upper triangle: the answer ($R$). Lower triangle: the recipe for $Q$. These two objects share one array.

Implementation

import numpy as np

def householder_qr(A):
    A = A.astype(float).copy()
    m, n = A.shape
    tau = np.zeros(n)
    for j in range(n):
        x = A[j:, j]
        v = x.copy()
        # sign convention: avoid cancellation when x[0] ≈ ||x||
        v[0] += (1.0 if v[0] >= 0 else -1.0) * np.linalg.norm(x)
        v /= v[0]                                      # normalize so v[0] = 1
        tau[j] = 2.0 / (v @ v)
        A[j:, j:] -= tau[j] * np.outer(v, v @ A[j:, j:])  # rank-1 update
        A[j+1:, j] = v[1:]                            # pack v into lower triangle
    return A, tau                                      # upper = R, lower = v's, tau separate

Line 7 is the sign convention fix (np.sign(0) = 0 in NumPy; this version defaults to $+1$ when $v[0]=0$). Line 10 is the entire reflector application. Lines 11–12 are the storage trick.

To apply $Q^T$ to a vector $b$ given the packed output:

def apply_Qt(QR, tau, b):
    b = b.copy().astype(float)
    m, n = QR.shape
    for j in range(n):
        v = np.concatenate([[1.0], QR[j+1:, j]])      # reconstruct v with implied 1
        b[j:] -= tau[j] * v * (v @ b[j:])
    return b

Part 6 — Concept Map

A diagram of the dependency graph of ideas across three clusters — geometry, mechanism, and storage & use — with labelled cross-cluster connections (sign rule → numerical stability; rank-1 update → packed array; column-by-column → zeros are free; the exact coincidence zero count = v entry count). SVG is pending.

Known remaining issues

The following items are identified for completion:

1. Geometric diagram. A 2D SVG illustration showing: vector $x$ as an arrow from origin; target $\pm \Vert x\Vert e_1$ on the horizontal axis; $v$ as the vector from target to $x$; the hyperplane as the perpendicular bisector; $Hx$ landing on the axis. Should make visually obvious that zeros appear from alignment, not force.
2. Numeric example step 2. The $2\times 2$ subblock update after constructing $v_2$ should be carried through to produce exact values of $r_{22}$, $r_{23}$, $r_{33}$, and $v_2(2)$. Currently only $r_{22}=\sqrt{5}/2$ is stated.
3. Concept map SVG. The SVG described in Part 6 needs to be built and embedded.
4. Invariant ordering audit. Confirm that no invariant uses a term introduced in a later invariant. Current known issue: check that the ordering of 13–24 is strict (definitions before use).