QR Decomposition

Why QR Decomposition Matters

Suppose three measurements describe a model with only two unknown coefficients. Noise can make the equations incompatible: there may be no vector $x$ for which $Ax=b$ exactly. We still want the coefficients that minimize the Euclidean error $\lVert Ax-b\rVert_2$. QR decomposition separates the geometry of the model's column space from a triangular system that is easy to solve.

The columns of $Q$ form perpendicular unit directions. The upper triangular matrix $R$ describes the original columns in those directions. This lesson connects three ways to obtain them: Gram–Schmidt removes projections, Householder reflects whole column tails, and Givens rotates two coordinates at a time. Our complete implementations use Householder reflections with explicit column pivoting, then solve a least-squares problem whose optimal residual is deliberately nonzero.

You will learn to distinguish a valid nonzero fitting error from an implementation error, trace the transformations, undo a column permutation correctly, and check both reconstruction and orthogonality. QR also appears in polynomial fitting, numerical optimization, orthogonal basis construction, and eigenvalue algorithms.

Core Concepts and Terminology

For a real matrix $A\in\mathbb{R}^{m\times n}$ with $m\ge n$ and independent columns, the thin QR factorization is $A=QR$, where $Q\in\mathbb{R}^{m\times n}$, $Q^TQ=I_n$, and $R\in\mathbb{R}^{n\times n}$ is upper triangular. Thin $Q$ has orthonormal columns but is usually not square, so $QQ^T$ is a projection rather than the identity.

Column pivoting reorders the variables: write $AP=QR$. Our array $p$ means column $j$ of $AP$ is original column $p_j$ of $A$, or $P_{p_j,j}=1$. To solve, set $z=P^Tx$, calculate $c=Q^Tb$, solve $Rz=c$, and scatter the result back with $x_{p_j}=z_j$. Pivoting swaps columns, not equations or entries of $b$.

A full QR factorization instead uses an $m\times m$ orthogonal matrix $Q_{\mathrm{full}}$ and an $m\times n$ upper trapezoidal factor, whose last $m-n$ rows are zero. Its first $n$ columns are thin $Q$. The remaining columns span the perpendicular complement of the model's column space.

In exact arithmetic, the least-squares solution satisfies

$$ Ax=QQ^Tb, \qquad A^T(Ax-b)=0. $$

Thus the residual is perpendicular to every column of $A$. Its norm need not be zero. Full column rank makes the minimizer unique; this lesson's code rejects a column tail that is too small at the chosen tolerance instead of returning a rank-deficient or minimum-norm solution.

The LAPACK DGEQP3 reference documents the same $AP=QR$ convention and a representation of $Q$ as stored reflectors. Our teaching implementation explicitly constructs thin $Q$ and uses a simple rank-rejection policy.

How QR Decomposition Works, Step by Step

Gram–Schmidt: build directions from projections. For column $a_j$, classical Gram–Schmidt computes $r_{ij}=q_i^Ta_j$ for each earlier $i$, subtracts their combined projections, and normalizes the remainder:

$$ w=a_j-\sum_{i=0}^{j-1}r_{ij}q_i, \qquad r_{jj}=\lVert w\rVert_2, \qquad q_j=w/r_{jj}. $$

Modified Gram–Schmidt starts with $w=a_j$ and repeatedly computes $r_{ij}=q_i^Tw$ from the current remainder, followed immediately by $w\leftarrow w-r_{ij}q_i$. They agree in exact arithmetic but differ in rounding behavior; modified Gram–Schmidt usually preserves orthogonality better. Reorthogonalization makes an additional projection pass when needed. A zero or rejected remainder must be handled before normalization.

For the sample's original columns $(1,1,0)^T$ and $(2,0,2)^T$, unpivoted Gram–Schmidt gives $q_0=(1,1,0)^T/\sqrt2$, $r_{01}=\sqrt2$, remainder $(1,-1,2)^T$, and $q_1=(1,-1,2)^T/\sqrt6$, with $r_{11}=\sqrt6$. This is a valid unpivoted factorization; the pivoted Householder factors below use a different column order and signs.

Householder: reflect tails into coordinate axes. This is the algorithm implemented in the three panels:

1. Validate a finite rectangular input with $m\ge n\ge1$. Set $s=\max_{i,j}\lvert A_{ij}\rvert$ and reject $s=0$. Work on $W=A/s$, preserving the input. Let $\rho$ be the largest initial column norm of $W$; keep it fixed.
2. At stage $k$, recompute $\lVert W_{k:m-1,j}\rVert_2$ for each remaining column $j\ge k$. Select the largest, retaining the first on an exact tie. Reject if the largest norm divided by $\rho$ is at most $\tau$.
3. Swap the chosen entire columns of $W$, including entries above row $k$, and update $p$. Previously computed rows of $R$ must follow their variables too.
4. Extract the active column tail $w$. Choose $\alpha=-\lVert w\rVert_2$ if $w_0\ge0$, otherwise $\alpha=\lVert w\rVert_2$. Set $u=w-\alpha e_0$ and $v=u/\lVert u\rVert_2$. The opposite sign avoids cancellation in the first component of $u$.
5. Apply $H_k=I-2vv^T$ to the active rows and remaining columns: for each column tail $t$, calculate $d=v^Tt$ and replace $t$ by $t-2vd$. Store the reflector, explicitly set the eliminated entries to zero, and set the pivot to $\alpha$.
6. After all stages, extract the first $n$ rows and rescale them by $s$ to obtain $R$. Form thin $Q=H_0H_1\cdots H_{n-1}E$, where $E$ is the first $n$ columns of $I_m$. Applying the stored reflectors to $E$ requires processing them in reverse order, since each update multiplies on the left.
7. For each right-hand side, calculate $c=Q^Tb$, back-substitute for $z$, then assign $x_{p_j}=z_j$. Verify the fit using the original $A$ and $b$.

Givens: rotate two coordinates. For entries $a$ and $b$ in a target column, set $r=\operatorname{hypot}(a,b)$, $c=a/r$, and $s_g=b/r$. When $r\ne0$,

$$ \begin{pmatrix}c&s_g\\-s_g&c\end{pmatrix} \begin{pmatrix}a\\b\end{pmatrix} =\begin{pmatrix}r\\0\end{pmatrix}. $$

For $(a,b)=(3,4)$, take $r=5$, $c=3/5$, and $s_g=4/5$. Apply that rotation to the same two rows of all affected columns, then continue clearing entries below the diagonal. If both entries are zero, use the identity rotation. The product of the transposed rotations builds $Q$. Givens is useful when changes or existing zeros make localized two-row work attractive.

flowchart TD
    A["Validate A; scale; initialize column order"] --> B{"More columns?"}
    B -->|Yes| C["Recompute active norms; choose column pivot"]
    C --> D{"Norm passes tolerance?"}
    D -->|No| E["Report rejected column rank"]
    D -->|Yes| F["Swap whole columns and order entries"]
    F --> G["Build stable reflector; clear tail; store reflector"]
    G --> H["Advance column"]
    H --> B
    B -->|No| I["Rescale R; form thin Q using reverse reflector order"]
    I --> J["Compute Q transpose b; back solve R z"]
    J --> K["Scatter z into original variable order"]
    K --> L["Check reconstruction, orthogonality, and fit stationarity"]

Worked Example with a Full Trace

Use

$$ A=\begin{pmatrix}1&2\\1&0\\0&2\end{pmatrix}, \qquad b=\begin{pmatrix}6\\0\\3\end{pmatrix}. $$

We show matrix states in original units for readability. The code scales by $s=2$ internally, so its working entries and pivot values are half of those displayed; normalized reflectors and pivot ratios are identical. Its fixed reference norm is $\rho=\sqrt2$, corresponding to the original maximum column norm $2\sqrt2$.

Stage $k=0$: choose a column. Original column norms are $\sqrt2$ and $2\sqrt2$. Choose column $1$ and swap it with column $0$, giving $p=(1,0)$ and

$$ AP=\begin{pmatrix}2&1\\0&1\\2&0\end{pmatrix}. $$

The active tail is $w=(2,0,2)^T$, with norm $2\sqrt2$. Its pivot ratio is $1$. Choose $\alpha_0=-2\sqrt2$. A normalized reflector vector and the resulting full reflector are

$$ v_0=\frac{(1+\sqrt2,0,1)^T}{\sqrt{4+2\sqrt2}}, \qquad H_0=\begin{pmatrix} -1/\sqrt2&0&-1/\sqrt2\\ 0&1&0\\ -1/\sqrt2&0&1/\sqrt2 \end{pmatrix}. $$

Applying it to both columns clears the first tail:

$$ W_1=H_0AP= \begin{pmatrix}-2\sqrt2&-1/\sqrt2\\0&1\\0&-1/\sqrt2\end{pmatrix}. $$

Stage $k=1$: clear the remaining tail. Only one column remains. Its active tail is $(1,-1/\sqrt2)^T$, with norm $\sqrt{3/2}=\sqrt6/2$. Relative to $2\sqrt2$, its pivot ratio is $\sqrt3/4\gt10^{-12}$. No further column swap occurs. Choose $\alpha_1=-\sqrt{3/2}$ and normalized active reflector vector

$$ v_1=\frac{(1+\sqrt{3/2},-1/\sqrt2)^T}{\sqrt{3+\sqrt6}}. $$

Embedded in three dimensions, its reflector is

$$ H_1=\begin{pmatrix} 1&0&0\\ 0&-\sqrt{2/3}&1/\sqrt3\\ 0&1/\sqrt3&\sqrt{2/3} \end{pmatrix}, \qquad W_2=H_1W_1= \begin{pmatrix}-2\sqrt2&-1/\sqrt2\\0&-\sqrt6/2\\0&0\end{pmatrix}. $$

Recover the thin factors. The top two rows of $W_2$ are $R$. The first two columns of $H_0H_1$ are $Q$:

$$ Q=\begin{pmatrix} -1/\sqrt2&-1/\sqrt6\\ 0&-2/\sqrt6\\ -1/\sqrt2&1/\sqrt6 \end{pmatrix}, \qquad R=\begin{pmatrix}-2\sqrt2&-1/\sqrt2\\0&-\sqrt6/2\end{pmatrix}. $$

The column norms of $Q$ are $1$, and their dot product is $1/\sqrt{12}-1/\sqrt{12}=0$. Multiplication gives $QR=AP$, including the original column $(1,1,0)^T$ as the second product column.

Transform the measurements. The first reflection gives $H_0b=(-9/\sqrt2,0,-3/\sqrt2)^T$. The second gives $H_1H_0b=(-9/\sqrt2,-3/\sqrt6,-\sqrt3)^T$. Its first two entries are $c=Q^Tb$; the last entry is the irreducible perpendicular component.

Back-substitute and restore variables. Solve $Rz=c$. The last equation gives $z_1=1$. The first then gives $z_0=2$. These are in pivoted variable order: $x_{p_0}=x_1=z_0=2$ and $x_{p_1}=x_0=z_1=1$. Thus $x=(1,2)^T$.

Check the fit. We obtain $Ax=(5,1,4)^T$ and $e=Ax-b=(-1,1,1)^T$. Therefore $\lVert e\rVert_2=\sqrt3$, and $A^Te=(0,0)^T$. The nonzero residual is correct: it is perpendicular to the model's column space and cannot be reduced by changing $x$.

Implementation in Python, C++, and Java

Each panel is a complete program with the same Householder and column-pivoting logic. Use Python 3, C++14 or newer, or save the Java program as Main.java. Inputs are nested real arrays with $m\ge n\ge1$; the code preserves $A$ and $b$. The factorization returns explicit thin $Q$, square $R$, and a zero-based final column order $p$. The solve routine expects an unchanged successful factorization and can be reused for another $b$.

All norms use repeated hypot calls rather than a direct sum of squares. Uniform scaling keeps reflector construction in a manageable range; the final $R$ is rescaled and checked for finite entries and representable nonzero pivots. The default $\tau=10^{-12}$ compares each active norm with the largest initial normalized column norm. This rejects deficient or numerically small directions, but it is not a singular-value calculation or a guarantee of accurate rank estimation. Extreme dynamic range can still cause underflow or loss of information; this educational code does not replace a carefully scaled numerical library.

The diagnostics check three different properties. The relative reconstruction error is $\max_{i,j}\lvert(QR)_{ij}-A_{i,p_j}\rvert/s$. Orthogonality error is the largest entry magnitude of $Q^TQ-I$. The stationarity error is $\lVert A^T(Ax-b)\rVert_\infty$; it should be small even though the fit residual norm is $\sqrt3$. The sample's absolute stationarity threshold of $10^{-10}$ assumes these modest data scales. Diagnostic dot products are simple floating-point sums intended for moderate inputs, not overflow-safe error estimators for every finite input.

import math


def norm(v):
    result = 0.0
    for value in v:
        result = math.hypot(result, value)
    return result


def qr_factor(a, tol=1e-12):
    if not math.isfinite(tol) or not 0.0 < tol < 1.0:
        raise ValueError("tol must be finite and between 0 and 1")
    m = len(a)
    n = len(a[0]) if m else 0
    if n == 0 or m < n or any(len(row) != n for row in a):
        raise ValueError("A must be rectangular with m >= n >= 1")
    if any(not math.isfinite(v) for row in a for v in row):
        raise ValueError("A must contain finite values")
    scale = max(abs(v) for row in a for v in row)
    if scale == 0.0:
        raise ValueError("zero matrix has deficient column rank")
    w = [[v / scale for v in row] for row in a]
    p = list(range(n))
    reference = max(norm([w[i][j] for i in range(m)]) for j in range(n))
    reflectors = []
    for k in range(n):
        pivot, best = k, -1.0
        for j in range(k, n):
            candidate = norm([w[i][j] for i in range(k, m)])
            if candidate > best:
                pivot, best = j, candidate
        if best / reference <= tol:
            raise ValueError("column rank rejected at the chosen tolerance")
        if pivot != k:
            for row in w:
                row[k], row[pivot] = row[pivot], row[k]
            p[k], p[pivot] = p[pivot], p[k]
        v = [w[i][k] for i in range(k, m)]
        alpha = -best if v[0] >= 0.0 else best
        v[0] -= alpha
        length = norm(v)
        v = [value / length for value in v]
        for j in range(k, n):
            dot = sum(v[t] * w[k + t][j] for t in range(len(v)))
            for t in range(len(v)):
                w[k + t][j] -= 2.0 * v[t] * dot
        w[k][k] = alpha
        for i in range(k + 1, m):
            w[i][k] = 0.0
        reflectors.append(v)
    # Form only the first n columns of H0 H1 ... H(n-1).
    q = [[float(i == j) for j in range(n)] for i in range(m)]
    for k in range(n - 1, -1, -1):
        v = reflectors[k]
        for j in range(n):
            dot = sum(v[t] * q[k + t][j] for t in range(len(v)))
            for t in range(len(v)):
                q[k + t][j] -= 2.0 * v[t] * dot
    r = [[w[i][j] * scale for j in range(n)] for i in range(n)]
    if any(not math.isfinite(v) for row in r for v in row):
        raise ValueError("R is not representable as finite values")
    if any(r[i][i] == 0.0 for i in range(n)):
        raise ValueError("R has an unrepresentable zero pivot")
    return q, r, p


def least_squares(factors, b):
    # factors must be an unchanged successful result of qr_factor.
    q, r, p = factors
    m, n = len(q), len(p)
    if len(b) != m or any(not math.isfinite(v) for v in b):
        raise ValueError("b must have m finite entries")
    c = [sum(q[i][j] * b[i] for i in range(m)) for j in range(n)]
    if any(not math.isfinite(v) for v in c):
        raise ValueError("nonfinite transformed right-hand side")
    z = [0.0] * n
    for i in range(n - 1, -1, -1):
        value = c[i]
        for j in range(i + 1, n):
            value -= r[i][j] * z[j]
        z[i] = value / r[i][i]
        if not math.isfinite(z[i]):
            raise ValueError("nonfinite triangular solve")
    x = [0.0] * n
    for j in range(n):
        x[p[j]] = z[j]
    return x


def diagnostics(a, factors, x, b):
    q, r, p = factors
    m, n = len(a), len(p)
    scale = max(abs(v) for row in a for v in row)
    fe = max(abs(sum(q[i][k] * r[k][j] for k in range(n)) - a[i][p[j]])
             for i in range(m) for j in range(n)) / scale
    oe = max(abs(sum(q[k][i] * q[k][j] for k in range(m)) - float(i == j))
             for i in range(n) for j in range(n))
    residual = [sum(a[i][j] * x[j] for j in range(n)) - b[i] for i in range(m)]
    stationarity = max(abs(sum(a[i][j] * residual[i] for i in range(m)))
                       for j in range(n))
    return fe, oe, norm(residual), stationarity


def main():
    a = [[1.0, 2.0], [1.0, 0.0], [0.0, 2.0]]
    b = [6.0, 0.0, 3.0]
    factors = qr_factor(a)
    x = least_squares(factors, b)
    fe, oe, residual, stationarity = diagnostics(a, factors, x, b)
    print("column order = " + str(factors[2]))
    print("x = [" + ", ".join(f"{v:.6f}" for v in x) + "]")
    print(f"fit residual norm = {residual:.6f}")
    print(f"relative factor error = {fe:.6f}")
    print(f"orthogonality error = {oe:.6f}")
    print(f"stationarity error = {stationarity:.6f}")
    good = max(fe, oe, stationarity, abs(residual - math.sqrt(3.0))) < 1e-10
    print("check = " + ("PASS" if good else "FAIL"))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <stdexcept>
#include <vector>

using Matrix = std::vector<std::vector<double>>;
struct Factors { Matrix q, r; std::vector<std::size_t> p; };

double norm(const std::vector<double>& v) {
    double result = 0.0;
    for (double value : v) result = std::hypot(result, value);
    return result;
}

double column_norm(const Matrix& w, std::size_t k, std::size_t j) {
    double result = 0.0;
    for (std::size_t i = k; i < w.size(); ++i) result = std::hypot(result, w[i][j]);
    return result;
}

Factors qr_factor(const Matrix& a, double tol = 1e-12) {
    if (!std::isfinite(tol) || !(tol > 0.0 && tol < 1.0))
        throw std::invalid_argument("tol must be finite and between 0 and 1");
    const std::size_t m = a.size(), n = m ? a[0].size() : 0;
    if (n == 0 || m < n) throw std::invalid_argument("require m >= n >= 1");
    double scale = 0.0;
    for (const auto& row : a) {
        if (row.size() != n) throw std::invalid_argument("A must be rectangular");
        for (double v : row) {
            if (!std::isfinite(v)) throw std::invalid_argument("A must be finite");
            scale = std::max(scale, std::abs(v));
        }
    }
    if (scale == 0.0) throw std::runtime_error("zero matrix has deficient column rank");
    Matrix w = a, reflectors;
    for (auto& row : w) for (auto& v : row) v /= scale;
    std::vector<std::size_t> p(n);
    std::iota(p.begin(), p.end(), 0);
    double reference = 0.0;
    for (std::size_t j = 0; j < n; ++j) reference = std::max(reference, column_norm(w, 0, j));
    for (std::size_t k = 0; k < n; ++k) {
        std::size_t pivot = k;
        double best = -1.0;
        for (std::size_t j = k; j < n; ++j) {
            const double candidate = column_norm(w, k, j);
            if (candidate > best) { pivot = j; best = candidate; }
        }
        if (best / reference <= tol)
            throw std::runtime_error("column rank rejected at the chosen tolerance");
        if (pivot != k) {
            for (auto& row : w) std::swap(row[k], row[pivot]);
            std::swap(p[k], p[pivot]);
        }
        std::vector<double> v(m - k);
        for (std::size_t t = 0; t < v.size(); ++t) v[t] = w[k + t][k];
        const double alpha = v[0] >= 0.0 ? -best : best;
        v[0] -= alpha;
        const double length = norm(v);
        for (auto& value : v) value /= length;
        for (std::size_t j = k; j < n; ++j) {
            double dot = 0.0;
            for (std::size_t t = 0; t < v.size(); ++t) dot += v[t] * w[k + t][j];
            for (std::size_t t = 0; t < v.size(); ++t) w[k + t][j] -= 2.0 * v[t] * dot;
        }
        w[k][k] = alpha;
        for (std::size_t i = k + 1; i < m; ++i) w[i][k] = 0.0;
        reflectors.push_back(v);
    }
    Matrix q(m, std::vector<double>(n, 0.0)), r(n, std::vector<double>(n));
    for (std::size_t i = 0; i < n; ++i) q[i][i] = 1.0;
    for (std::size_t k = n; k-- > 0;) {
        const auto& v = reflectors[k];
        for (std::size_t j = 0; j < n; ++j) {
            double dot = 0.0;
            for (std::size_t t = 0; t < v.size(); ++t) dot += v[t] * q[k + t][j];
            for (std::size_t t = 0; t < v.size(); ++t) q[k + t][j] -= 2.0 * v[t] * dot;
        }
    }
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) {
            r[i][j] = w[i][j] * scale;
            if (!std::isfinite(r[i][j])) throw std::runtime_error("R is not finite");
        }
    for (std::size_t i = 0; i < n; ++i)
        if (r[i][i] == 0.0) throw std::runtime_error("R has an unrepresentable zero pivot");
    return {q, r, p};
}

std::vector<double> least_squares(const Factors& f, const std::vector<double>& b) {
    // f must be an unchanged successful result of qr_factor.
    const std::size_t m = f.q.size(), n = f.p.size();
    if (b.size() != m) throw std::invalid_argument("b must have m entries");
    for (double v : b) if (!std::isfinite(v)) throw std::invalid_argument("b must be finite");
    std::vector<double> c(n, 0.0), z(n, 0.0), x(n);
    for (std::size_t j = 0; j < n; ++j) {
        for (std::size_t i = 0; i < m; ++i) c[j] += f.q[i][j] * b[i];
        if (!std::isfinite(c[j])) throw std::runtime_error("nonfinite transformed right-hand side");
    }
    for (std::size_t i = n; i-- > 0;) {
        double value = c[i];
        for (std::size_t j = i + 1; j < n; ++j) value -= f.r[i][j] * z[j];
        z[i] = value / f.r[i][i];
        if (!std::isfinite(z[i])) throw std::runtime_error("nonfinite triangular solve");
    }
    for (std::size_t j = 0; j < n; ++j) x[f.p[j]] = z[j];
    return x;
}

std::vector<double> diagnostics(const Matrix& a, const Factors& f,
                                const std::vector<double>& x, const std::vector<double>& b) {
    const std::size_t m = a.size(), n = f.p.size();
    double scale = 0.0, fe = 0.0, oe = 0.0, stationarity = 0.0;
    for (const auto& row : a) for (double v : row) scale = std::max(scale, std::abs(v));
    for (std::size_t i = 0; i < m; ++i)
        for (std::size_t j = 0; j < n; ++j) {
            double value = 0.0;
            for (std::size_t k = 0; k < n; ++k) value += f.q[i][k] * f.r[k][j];
            fe = std::max(fe, std::abs(value - a[i][f.p[j]]));
        }
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) {
            double value = 0.0;
            for (std::size_t k = 0; k < m; ++k) value += f.q[k][i] * f.q[k][j];
            oe = std::max(oe, std::abs(value - (i == j ? 1.0 : 0.0)));
        }
    std::vector<double> residual(m);
    for (std::size_t i = 0; i < m; ++i) {
        double value = 0.0;
        for (std::size_t j = 0; j < n; ++j) value += a[i][j] * x[j];
        residual[i] = value - b[i];
    }
    for (std::size_t j = 0; j < n; ++j) {
        double value = 0.0;
        for (std::size_t i = 0; i < m; ++i) value += a[i][j] * residual[i];
        stationarity = std::max(stationarity, std::abs(value));
    }
    return {fe / scale, oe, norm(residual), stationarity};
}

int main() {
    const Matrix a = { {1.0, 2.0}, {1.0, 0.0}, {0.0, 2.0} };
    const std::vector<double> b = {6.0, 0.0, 3.0};
    const auto f = qr_factor(a);
    const auto x = least_squares(f, b), d = diagnostics(a, f, x, b);
    std::cout << "column order = [";
    for (std::size_t i = 0; i < f.p.size(); ++i) {
        if (i) std::cout << ", ";
        std::cout << f.p[i];
    }
    std::cout << "]\nx = [" << std::fixed << std::setprecision(6);
    for (std::size_t i = 0; i < x.size(); ++i) {
        if (i) std::cout << ", ";
        std::cout << x[i];
    }
    std::cout << "]\nfit residual norm = " << d[2] << "\nrelative factor error = " << d[0]
              << "\northogonality error = " << d[1] << "\nstationarity error = " << d[3];
    const bool good = std::max({d[0], d[1], d[3], std::abs(d[2] - std::sqrt(3.0))}) < 1e-10;
    std::cout << "\ncheck = " << (good ? "PASS" : "FAIL") << '\n';
}

import java.util.Arrays;
import java.util.Locale;

public class Main {
    static class Factors {
        final double[][] q, r;
        final int[] p;
        Factors(double[][] q, double[][] r, int[] p) { this.q = q; this.r = r; this.p = p; }
    }

    static double norm(double[] v) {
        double result = 0.0;
        for (double value : v) result = Math.hypot(result, value);
        return result;
    }

    static double columnNorm(double[][] w, int k, int j) {
        double result = 0.0;
        for (int i = k; i < w.length; ++i) result = Math.hypot(result, w[i][j]);
        return result;
    }

    static Factors qrFactor(double[][] a, double tol) {
        if (!Double.isFinite(tol) || !(tol > 0.0 && tol < 1.0))
            throw new IllegalArgumentException("tol must be finite and between 0 and 1");
        int m = a.length, n = m > 0 && a[0] != null ? a[0].length : 0;
        if (n == 0 || m < n) throw new IllegalArgumentException("require m >= n >= 1");
        double scale = 0.0;
        for (double[] row : a) {
            if (row == null || row.length != n) throw new IllegalArgumentException("A must be rectangular");
            for (double v : row) {
                if (!Double.isFinite(v)) throw new IllegalArgumentException("A must be finite");
                scale = Math.max(scale, Math.abs(v));
            }
        }
        if (scale == 0.0) throw new IllegalArgumentException("zero matrix has deficient column rank");
        double[][] w = new double[m][n], reflectors = new double[n][];
        for (int i = 0; i < m; ++i) for (int j = 0; j < n; ++j) w[i][j] = a[i][j] / scale;
        int[] p = new int[n];
        double reference = 0.0;
        for (int j = 0; j < n; ++j) { p[j] = j; reference = Math.max(reference, columnNorm(w, 0, j)); }
        for (int k = 0; k < n; ++k) {
            int pivot = k;
            double best = -1.0;
            for (int j = k; j < n; ++j) {
                double candidate = columnNorm(w, k, j);
                if (candidate > best) { pivot = j; best = candidate; }
            }
            if (best / reference <= tol)
                throw new IllegalArgumentException("column rank rejected at the chosen tolerance");
            if (pivot != k) {
                for (double[] row : w) { double value = row[k]; row[k] = row[pivot]; row[pivot] = value; }
                int index = p[k]; p[k] = p[pivot]; p[pivot] = index;
            }
            double[] v = new double[m - k];
            for (int t = 0; t < v.length; ++t) v[t] = w[k + t][k];
            double alpha = v[0] >= 0.0 ? -best : best;
            v[0] -= alpha;
            double length = norm(v);
            for (int t = 0; t < v.length; ++t) v[t] /= length;
            for (int j = k; j < n; ++j) {
                double dot = 0.0;
                for (int t = 0; t < v.length; ++t) dot += v[t] * w[k + t][j];
                for (int t = 0; t < v.length; ++t) w[k + t][j] -= 2.0 * v[t] * dot;
            }
            w[k][k] = alpha;
            for (int i = k + 1; i < m; ++i) w[i][k] = 0.0;
            reflectors[k] = v;
        }
        double[][] q = new double[m][n], r = new double[n][n];
        for (int i = 0; i < n; ++i) q[i][i] = 1.0;
        for (int k = n - 1; k >= 0; --k) {
            double[] v = reflectors[k];
            for (int j = 0; j < n; ++j) {
                double dot = 0.0;
                for (int t = 0; t < v.length; ++t) dot += v[t] * q[k + t][j];
                for (int t = 0; t < v.length; ++t) q[k + t][j] -= 2.0 * v[t] * dot;
            }
        }
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j) {
                r[i][j] = w[i][j] * scale;
                if (!Double.isFinite(r[i][j])) throw new IllegalArgumentException("R is not finite");
            }
        for (int i = 0; i < n; ++i)
            if (r[i][i] == 0.0) throw new IllegalArgumentException("R has an unrepresentable zero pivot");
        return new Factors(q, r, p);
    }

    static double[] leastSquares(Factors f, double[] b) {
        // f must be an unchanged successful result of qrFactor.
        int m = f.q.length, n = f.p.length;
        if (b.length != m) throw new IllegalArgumentException("b must have m entries");
        for (double v : b) if (!Double.isFinite(v)) throw new IllegalArgumentException("b must be finite");
        double[] c = new double[n], z = new double[n], x = new double[n];
        for (int j = 0; j < n; ++j) {
            for (int i = 0; i < m; ++i) c[j] += f.q[i][j] * b[i];
            if (!Double.isFinite(c[j])) throw new IllegalArgumentException("nonfinite transformed right-hand side");
        }
        for (int i = n - 1; i >= 0; --i) {
            double value = c[i];
            for (int j = i + 1; j < n; ++j) value -= f.r[i][j] * z[j];
            z[i] = value / f.r[i][i];
            if (!Double.isFinite(z[i])) throw new IllegalArgumentException("nonfinite triangular solve");
        }
        for (int j = 0; j < n; ++j) x[f.p[j]] = z[j];
        return x;
    }

    static double[] diagnostics(double[][] a, Factors f, double[] x, double[] b) {
        int m = a.length, n = f.p.length;
        double scale = 0.0, fe = 0.0, oe = 0.0, stationarity = 0.0;
        for (double[] row : a) for (double v : row) scale = Math.max(scale, Math.abs(v));
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < n; ++j) {
                double value = 0.0;
                for (int k = 0; k < n; ++k) value += f.q[i][k] * f.r[k][j];
                fe = Math.max(fe, Math.abs(value - a[i][f.p[j]]));
            }
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j) {
                double value = 0.0;
                for (int k = 0; k < m; ++k) value += f.q[k][i] * f.q[k][j];
                oe = Math.max(oe, Math.abs(value - (i == j ? 1.0 : 0.0)));
            }
        double[] residual = new double[m];
        for (int i = 0; i < m; ++i) {
            double value = 0.0;
            for (int j = 0; j < n; ++j) value += a[i][j] * x[j];
            residual[i] = value - b[i];
        }
        for (int j = 0; j < n; ++j) {
            double value = 0.0;
            for (int i = 0; i < m; ++i) value += a[i][j] * residual[i];
            stationarity = Math.max(stationarity, Math.abs(value));
        }
        return new double[]{fe / scale, oe, norm(residual), stationarity};
    }

    public static void main(String[] args) {
        double[][] a = { {1.0, 2.0}, {1.0, 0.0}, {0.0, 2.0} };
        double[] b = {6.0, 0.0, 3.0};
        Factors f = qrFactor(a, 1e-12);
        double[] x = leastSquares(f, b), d = diagnostics(a, f, x, b);
        System.out.println("column order = " + Arrays.toString(f.p));
        System.out.printf(Locale.ROOT, "x = [%.6f, %.6f]%n", x[0], x[1]);
        System.out.printf(Locale.ROOT, "fit residual norm = %.6f%nrelative factor error = %.6f%northogonality error = %.6f%nstationarity error = %.6f%n",
                          d[2], d[0], d[1], d[3]);
        boolean good = Math.max(Math.max(d[0], d[1]), Math.max(d[3], Math.abs(d[2] - Math.sqrt(3.0)))) < 1e-10;
        System.out.println("check = " + (good ? "PASS" : "FAIL"));
    }
}

All three programs print the same seven lines. Six decimal places hide the small nonzero roundoff in the diagnostic errors:

• column order = [1, 0]
• x = [1.000000, 2.000000]
• fit residual norm = 1.732051
• relative factor error = 0.000000
• orthogonality error = 0.000000
• stationarity error = 0.000000
• check = PASS

Complexity, Correctness, and Invariants

For a successful dense $m\times n$ input with $m\ge n$, the implementations take $\Theta(mn^2)$ time in the best, average, and worst cases. They recompute the remaining column norms, update the trailing matrix, and explicitly form thin $Q$ using the same loop structure regardless of existing zeros. A rejected direction can terminate early. Unpivoted Householder reduction alone has a leading arithmetic count near $2mn^2-2n^3/3$; column-norm recomputation and explicit construction of $Q$ add work, so that expression is not the total count for these programs.

The working matrix, reflector vectors, thin $Q$, and $R$ together use $\Theta(mn)$ storage because $n\le m$. A later solve needs $\Theta(mn+n^2)=\Theta(mn)$ time and $\Theta(n)$ additional workspace, excluding its input vector. Reusing the factors for $t$ right-hand sides costs $\Theta(mn^2+tmn)$ overall. Forming a full $m\times m$ orthogonal matrix would require more storage and is deliberately avoided.

In exact arithmetic, each normalized reflector satisfies $H_k^T=H_k$ and $H_k^TH_k=I$. After $k$ stages, the evolving relation is

$$ W_k=H_{k-1}\cdots H_0(A/s)P_k. $$

The first $k$ columns have zero entries below their diagonals. A new pivot swaps only remaining columns, preserving those zeros. The next reflector acts only on rows $k,\ldots,m-1$, clearing its chosen column tail without disturbing previously completed columns. After the final stage, the transformed matrix is a top triangular block followed by zero rows.

Reversing the orthogonal transformations gives $AP=Q_{\mathrm{full}}\begin{pmatrix}R\\0\end{pmatrix}=QR$. Because thin $Q$ consists of columns of an orthogonal matrix, $Q^TQ=I_n$. The code's reverse loop forms this product by left-multiplying the rectangular identity $E$ in the required order.

For any candidate $x$, write $x=Pz$. Orthogonal projection splits the squared error into two perpendicular parts:

$$ \lVert Ax-b\rVert_2^2 =\lVert Rz-Q^Tb\rVert_2^2 +\lVert(I_m-QQ^T)b\rVert_2^2. $$

The second term is independent of $z$. Nonzero accepted diagonal entries make $R$ invertible, so solving $Rz=Q^Tb$ uniquely makes the first term zero. Scattering $z$ back through $P$ yields the unique minimizer in original variable order. Floating-point arithmetic makes these identities approximate; reconstruction, orthogonality, and stationarity checks test different parts of that argument.

Variants, Pitfalls, and Optimizations

Do not lose the variable order. For this sample, the triangular result is $z=(2,1)^T$, but the answer in original coordinates is $x=(1,2)^T$. If you incorrectly return $x=z$, the fitted vector becomes $(4,2,2)^T$, its residual norm is $3$ instead of $\sqrt3$, and its stationarity error is $6$. For a general permutation, assign x[p[j]] = z[j]; guessing the inverse direction can silently fail. At later pivots, swap complete working columns so earlier rows of $R$ move with their variables.

Preserve orthogonality and use the stable sign. Nearly dependent columns produce a small remainder in Gram–Schmidt, where cancellation can magnify rounding error. Modified Gram–Schmidt and a second projection pass improve this behavior, while Householder is a usual dense choice. In a reflector, choosing the same sign as the first component can subtract nearly equal values: for $w=(1,\varepsilon)^T$ with tiny $\varepsilon$, $1-\lVert w\rVert_2$ can round to zero. The opposite sign uses a first component near $2$ instead.

Factor signs are not unique. Flipping one column of $Q$ and the corresponding row of $R$ leaves the product unchanged. Negative diagonal entries of $R$ in our example are valid. Comparing different implementations requires allowing for sign conventions, and pivot ties can also produce different valid column orders. Check mathematical identities before comparing raw factor entries.

Rank deficiency needs a separate solution policy. A tiny diagonal suggests a weak direction but does not by itself give an accurate condition number or a universal numerical rank. The LAPACK discussion of pivoted QR distinguishes factorization from rank determination. For rank-deficient or underdetermined problems, a basic solution is not necessarily the minimum-norm solution; use an appropriate complete orthogonal factorization or SVD when that is the requirement. Reducing $\tau$ blindly can amplify a poorly determined coefficient.

Choose the transformation for the workload. Givens rotations can update a triangular factor after adding a row, making them useful in online least squares. Sparse QR needs fill-reducing ordering because localized operations can still create nonzeros. Modified Gram–Schmidt appears naturally when building Krylov bases, often with reorthogonalization. Tall-skinny QR reduces communication on distributed tall matrices. Production dense implementations block reflectors into matrix operations, as described in the LAPACK Householder QR discussion.

Avoid unnecessary matrices. A solver can store reflectors and apply them directly to $b$ without explicitly forming $Q$ or $QQ^T$. Our explicit thin $Q$ makes the geometry and orthogonality check visible. Forming normal equations can square the condition number, so their lower storage footprint comes with a numerical tradeoff. Column rescaling can improve badly scaled data but changes coefficient coordinates and must be undone; row weighting changes the least-squares objective itself.

Practice Problems and Real-World Use

Begin by computing two-column orthogonalizations by hand, then compare methods on the same data. A full-rank random test can start from a known $x$ and set $b=Ax$. To test a nonzero least-squares residual independently, choose a vector $e$ perpendicular to the columns and use $b=Ax+e$; the same $x$ should remain the minimizer. Include column permutations with three or more entries so an incorrect scatter direction cannot hide behind a single swap.

For numerical tests, check reconstruction relative to matrix scale, orthogonality separately, and a scale-aware stationarity measure when data magnitudes change. Reuse a factorization for several right-hand sides, verify the original arrays remain intact, and include dependent columns, near-dependence, a zero leading component, and a one-column input.