QR Decomposition & Gram-Schmidt
Orthogonalizing vectors, and the numerically stable way to do least squares.
On this page
Beginner: orthogonalization means turning a set of vectors into a set that all point perpendicular to each other, without changing what they span (section 1.9). The classic recipe for doing this is the Gram-Schmidt process: take each new vector, subtract off whatever part of it points along the directions you already have, and keep only what's left over — which is automatically perpendicular to everything before it.
Intermediate: QR decomposition packages this idea into a matrix factorization: any matrix A can be written as A = QR, where Q has orthonormal columns (the Gram-Schmidt-ed, unit-length version of A's columns) and R is upper-triangular (recording exactly how much of each original column was "redundant" with the earlier ones).
Advanced: QR gives a numerically superior way to solve least-squares regression compared to the normal equation from section 1.6 (w = (XᵀX)⁻¹Xᵀy). Forming XᵀX squares the condition number of X, amplifying numerical error; solving via QR avoids that squaring entirely, which is why serious statistical software defaults to a QR-based solver rather than the textbook formula.
Gram-Schmidt, step by step, for vectors a₁, a₂: e₁ = a₁/‖a₁‖, then e₂ = (a₂ − (a₂·e₁)e₁), normalized. Each new vector only ever has the previous directions subtracted out.
Given a fixed direction e₁ (unit length) and a new vector a₂, define the projection coefficient c = a₂·e₁ and subtract it off:
Claim: u is exactly perpendicular to e₁. Proof: compute the dot product directly:
(using e₁·e₁ = ‖e₁‖² = 1, since e₁ is a unit vector). This holds for any choice of a₂ — c is defined precisely so this cancellation always happens, which is exactly the live check shown in the diagram's readout.
Where this is used: repeating this subtraction against every previously built direction, one at a time, is the entire Gram-Schmidt algorithm — each new vector is only ever guaranteed perpendicular to what came immediately before it, which is exactly why the process must go through the vectors in order.
Blue is the fixed reference direction. Drag orange — grey dashed is the projection being removed, green is what's left: always exactly perpendicular to blue.
np.linalg.lstsq is what you should reach for in practice — but knowing it's doing QR (not the textbook normal equation) explains why it's the more numerically trustworthy choice.
#include <cstdio>
#include <vector>
#include <cmath>
using Matrix = std::vector<std::vector<double>>;
using Vector = std::vector<double>;
// Modified Gram-Schmidt QR decomposition -- no linear-algebra library.
void qrDecompose(const Matrix& A, Matrix& Q, Matrix& R) {
int m = static_cast<int>(A.size());
int n = static_cast<int>(A[0].size());
std::vector<Vector> cols(n, Vector(m));
for (int j = 0; j < n; ++j)
for (int i = 0; i < m; ++i) cols[j][i] = A[i][j];
std::vector<Vector> qCols;
R.assign(n, Vector(n, 0.0));
for (int j = 0; j < n; ++j) {
Vector v = cols[j];
for (int i = 0; i < j; ++i) {
double dot = 0.0;
for (int k = 0; k < m; ++k) dot += qCols[i][k] * v[k];
R[i][j] = dot;
for (int k = 0; k < m; ++k) v[k] -= dot * qCols[i][k];
}
double norm = 0.0;
for (double vk : v) norm += vk * vk;
norm = std::sqrt(norm);
R[j][j] = norm;
for (double& vk : v) vk /= norm;
qCols.push_back(v);
}
Q.assign(m, Vector(n, 0.0));
for (int j = 0; j < n; ++j)
for (int i = 0; i < m; ++i) Q[i][j] = qCols[j][i];
}
Vector backSub(const Matrix& R, const Vector& y) {
int n = static_cast<int>(y.size());
Vector x(n, 0.0);
for (int i = n - 1; i >= 0; --i) {
double s = y[i];
for (int j = i + 1; j < n; ++j) s -= R[i][j] * x[j];
x[i] = s / R[i][i];
}
return x;
}
int main() {
Matrix X = {{1., 1.}, {1., 2.}, {1., 3.}, {1., 4.}};
Vector y = {2., 4., 5., 8.};
Matrix Q, R;
qrDecompose(X, Q, R);
// Q^T Q should be close to the identity.
double qtq[2][2] = {{0, 0}, {0, 0}};
for (size_t k = 0; k < Q.size(); ++k)
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j) qtq[i][j] += Q[k][i] * Q[k][j];
std::printf("Q^T Q = [[%.4f, %.4f], [%.4f, %.4f]]\n", qtq[0][0], qtq[0][1], qtq[1][0], qtq[1][1]);
Vector Qty(2, 0.0);
for (size_t k = 0; k < Q.size(); ++k)
for (int j = 0; j < 2; ++j) Qty[j] += Q[k][j] * y[k];
Vector w = backSub(R, Qty);
std::printf("w = [%.4f, %.4f]\n", w[0], w[1]);
return 0;
}- Robotics and computer graphics use QR (or the related Householder reflections) to keep a sequence of rotation matrices numerically "clean" — repeated multiplication of rotation matrices slowly accumulates floating-point drift away from true orthogonality, and re-orthogonalizing via QR fixes it.
- Every serious linear regression / least-squares solver (R's
lm(), scikit-learn'sLinearRegression, NumPy'slstsq) uses QR or SVD internally instead of the raw normal equation.
- Implementing "classical" Gram-Schmidt naively for many vectors — it's numerically unstable in practice; real libraries use modified Gram-Schmidt or Householder reflections, which are mathematically equivalent but far more stable in floating point.
- Solving least squares via
(XᵀX)⁻¹Xᵀydirectly in code that matters — prefernp.linalg.lstsqor an explicit QR/SVD-based solve.
Going deeper
QR decomposition isn't just for least squares — the QR algorithm (repeatedly factoring a matrix as QR, then multiplying the factors back together in reverse order, and iterating) is the actual method general-purpose numerical libraries use to compute eigenvalues (section 1.7) for matrices larger than 2×2 or 3×3. The "solve the characteristic polynomial" method taught in school is essentially never used in real software — it's numerically unreliable for anything but the smallest matrices.
At the master level: the columns of Q form an orthonormal basis for the same column space as A — meaning QR is simultaneously an orthogonalization procedure, a rank-revealing factorization, and (via the QR algorithm) an eigenvalue solver, which is a lot of mileage from one relatively simple idea.