Householder vs. Gram-Schmidt QR

Both factor A = Q R, but Householder reflections stay orthogonal to machine precision even for ill-conditioned A, while classical Gram-Schmidt can lose orthogonality badly – modified Gram-Schmidt is a partial (but not complete) remedy. This is the textbook motivation for preferring Householder in production numerical software.

import numpy as np

from mathematicskit.linalg import gram_schmidt_qr, householder_qr
from mathematicskit.linalg.systems.qr import orthogonality_error

A well-conditioned matrix: all three methods agree

rng = np.random.default_rng(0)
A = rng.uniform(-2, 2, size=(6, 4))
for label, result in (
    ("householder", householder_qr(A)),
    ("gram_schmidt (modified)", gram_schmidt_qr(A, modified=True)),
    ("gram_schmidt (classical)", gram_schmidt_qr(A, modified=False)),
):
    print(f"{label:>28s}: ||QR - A|| = {np.linalg.norm(result.Q @ result.R - A):.2e}, orthogonality error = {orthogonality_error(result.Q):.2e}")

             householder: ||QR - A|| = 2.36e-15, orthogonality error = 4.44e-16
 gram_schmidt (modified): ||QR - A|| = 3.38e-16, orthogonality error = 2.22e-16
gram_schmidt (classical): ||QR - A|| = 5.12e-16, orthogonality error = 2.22e-16

An ill-conditioned matrix: orthogonality breaks down differently

Nearly-collinear columns (Trefethen & Bau’s classic example) expose the gap between the three methods.

eps = 1e-8
A_ill = np.array([[1.0, 1.0, 1.0], [eps, 0.0, 0.0], [0.0, eps, 0.0], [0.0, 0.0, eps]])
for label, result in (
    ("householder", householder_qr(A_ill)),
    ("gram_schmidt (modified)", gram_schmidt_qr(A_ill, modified=True)),
    ("gram_schmidt (classical)", gram_schmidt_qr(A_ill, modified=False)),
):
    print(f"{label:>28s}: orthogonality error = {orthogonality_error(result.Q):.3e}")

             householder: orthogonality error = 6.661e-16
 gram_schmidt (modified): orthogonality error = 7.071e-09
gram_schmidt (classical): orthogonality error = 5.000e-01