Add QR decomposition

7 files changed, 730 insertions, 11 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 87cd6a3..1d2374c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -18,6 +18,7 @@ add_library(linear_algebra
     src/norms.cpp
     src/triangular_solve.cpp
     src/lu.cpp
+    src/qr.cpp
 )
 
 add_library(linear_algebra::core ALIAS linear_algebra)
@@ -65,6 +66,9 @@ option(LINEAR_ALGEBRA_BUILD_EXPERIMENTS "Build experiment programs" ON)
 if(LINEAR_ALGEBRA_BUILD_EXPERIMENTS)
     add_executable(pivoting_vs_no_pivoting experiments/pivoting_vs_no_pivoting.cpp)
     target_link_libraries(pivoting_vs_no_pivoting PRIVATE linear_algebra::core)
+
+    add_executable(hilbert_qr experiments/hilbert_qr.cpp)
+    target_link_libraries(hilbert_qr PRIVATE linear_algebra::core)
 endif()
 
 if(LINEAR_ALGEBRA_BUILD_TESTS)
@@ -88,6 +92,7 @@ if(LINEAR_ALGEBRA_BUILD_TESTS)
         tests/test_matrix.cpp
         tests/test_triangular_solve.cpp
         tests/test_lu.cpp
+        tests/test_qr.cpp
     )
 
     target_link_libraries(linear_algebra_tests
diff --git a/README.md b/README.md
index e95b936..172f1a7 100644
--- a/README.md
+++ b/README.md
@@ -25,9 +25,23 @@ Valid values are `AUTO`, `NONE`, `AVX`, `AVX2`, and `AVX512`. `AUTO` uses the co
 ctest --test-dir build --output-on-failure
 ```
 
+## What's implemented
+
+- Matrix / Vector core with SIMD matmul
+- Triangular solvers (forward / backward substitution)
+- LU factorization with partial pivoting (`lu_factor`, `lu_solve`)
+- QR factorization — classical GS, modified GS, and Householder (`qr_classical_gs`, `qr_modified_gs`, `qr_householder`)
+
 ## Run examples
 
 ```bash
 ./build/linear_system
 ./build/matmul
-```
\ No newline at end of file
+```
+
+## Run experiments
+
+```bash
+./build/pivoting_vs_no_pivoting   
+./build/hilbert_qr              
+```
diff --git a/experiments/hilbert_qr.cpp b/experiments/hilbert_qr.cpp
new file mode 100644
index 0000000..9209a30
--- /dev/null
+++ b/experiments/hilbert_qr.cpp
@@ -0,0 +1,167 @@
+// Experiment: QR methods on Hilbert matrices
+//
+// The Hilbert matrix H[i][j] = 1/(i+j+1) is the canonical ill-conditioned
+// dense matrix.  Its condition number grows roughly as (3.5 * e)^n / sqrt(n),
+// reaching ~10^13 at n=10 and ~10^18 at n=14.
+//
+// We compare classical GS, modified GS, and Householder QR on:
+//   - reconstruction error   ||A - QR||_F
+//   - orthogonality error    ||Q^T Q - I||_F
+//   - wall-clock time (minimum over several trials)
+
+#include "matrix.hpp"
+#include "qr.hpp"
+
+#include <chrono>
+#include <cmath>
+#include <cstddef>
+#include <functional>
+#include <iomanip>
+#include <iostream>
+#include <optional>
+#include <stdexcept>
+#include <string>
+
+using linalg::Matrix;
+using linalg::QRResult;
+
+// ---------------------------------------------------------------------------
+// Matrix construction
+// ---------------------------------------------------------------------------
+
+Matrix hilbert(std::size_t n) {
+    Matrix H(n, n);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            H(i, j) = 1.0 / static_cast<double>(i + j + 1);
+    return H;
+}
+
+// ---------------------------------------------------------------------------
+// Metrics
+// ---------------------------------------------------------------------------
+
+double reconstruction_error(const Matrix& A, const QRResult& qr) {
+    const std::size_t m = A.rows();
+    const std::size_t n = A.cols();
+    const Matrix QR = qr.Q * qr.R;
+    double err = 0.0;
+    for (std::size_t i = 0; i < m; ++i)
+        for (std::size_t j = 0; j < n; ++j) {
+            const double d = A(i, j) - QR(i, j);
+            err += d * d;
+        }
+    return std::sqrt(err);
+}
+
+double orthogonality_error(const QRResult& qr) {
+    const Matrix& Q = qr.Q;
+    const std::size_t n = Q.cols();
+    double err = 0.0;
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j) {
+            double s = 0.0;
+            for (std::size_t k = 0; k < Q.rows(); ++k) s += Q(k, i) * Q(k, j);
+            const double d = s - (i == j ? 1.0 : 0.0);
+            err += d * d;
+        }
+    return std::sqrt(err);
+}
+
+// ---------------------------------------------------------------------------
+// Timing
+// ---------------------------------------------------------------------------
+
+using Clock = std::chrono::high_resolution_clock;
+using Seconds = std::chrono::duration<double>;
+
+// Run fn() `trials` times, return minimum elapsed seconds.
+template<typename Fn>
+double min_time(Fn fn, int trials = 5) {
+    double best = 1e18;
+    for (int t = 0; t < trials; ++t) {
+        const auto t0 = Clock::now();
+        fn();
+        const auto t1 = Clock::now();
+        best = std::min(best, Seconds(t1 - t0).count());
+    }
+    return best;
+}
+
+// ---------------------------------------------------------------------------
+// Run one method on one size, return {recon, ortho, time} or nullopt on failure
+// ---------------------------------------------------------------------------
+
+using QRFn = std::function<QRResult(const Matrix&)>;
+
+struct Result { double recon, ortho, time_s; };
+
+std::optional<Result> measure(const Matrix& A, QRFn fn) {
+    try {
+        // Run once to get metrics.
+        const QRResult qr = fn(A);
+        const double re = reconstruction_error(A, qr);
+        const double oe = orthogonality_error(qr);
+        // Time over multiple trials.
+        const double t = min_time([&] { fn(A); });
+        return Result{re, oe, t};
+    } catch (const std::exception&) {
+        return std::nullopt;
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Pretty printing
+// ---------------------------------------------------------------------------
+
+void print_row(const std::string& method, std::optional<Result> r) {
+    std::cout << std::left << std::setw(16) << method;
+    if (!r) {
+        std::cout << "  FAILED (linearly dependent columns)\n";
+        return;
+    }
+    std::cout << std::scientific << std::setprecision(2)
+              << std::setw(14) << r->recon
+              << std::setw(14) << r->ortho
+              << std::fixed << std::setprecision(3)
+              << std::setw(10) << r->time_s * 1e6 << " µs\n";
+}
+
+// ---------------------------------------------------------------------------
+// main
+// ---------------------------------------------------------------------------
+
+int main() {
+    std::cout << std::string(70, '*') << "\n";
+    std::cout << "  Hilbert QR Experiment — comparing GS variants and Householder\n";
+    std::cout << std::string(70, '*') << "\n\n";
+    std::cout <<
+        "H[i][j] = 1/(i+j+1).  Condition number grows ~exponentially with n.\n"
+        "Orthogonality loss in classical GS tracks condition number directly.\n"
+        "Modified GS recovers ~half the lost digits.  Householder is unaffected.\n\n";
+
+    const std::size_t sizes[] = {2, 3, 4, 5, 6, 7, 8, 10, 12};
+
+    for (std::size_t n : sizes) {
+        const Matrix H = hilbert(n);
+
+        std::cout << std::string(70, '-') << "\n";
+        std::cout << "  n = " << n << "\n";
+        std::cout << std::string(70, '-') << "\n";
+        std::cout << std::left
+                  << std::setw(16) << "Method"
+                  << std::setw(14) << "||A-QR||_F"
+                  << std::setw(14) << "||QtQ-I||_F"
+                  << std::setw(10) << "Time\n";
+        std::cout << std::string(70, ' ') << "\n";
+
+        print_row("classical_gs",
+            measure(H, [](const Matrix& A) { return linalg::qr_classical_gs(A); }));
+        print_row("modified_gs",
+            measure(H, [](const Matrix& A) { return linalg::qr_modified_gs(A); }));
+        print_row("householder",
+            measure(H, [](const Matrix& A) { return linalg::qr_householder(A); }));
+    }
+
+    return 0;
+}
diff --git a/experiments/pivoting_vs_no_pivoting.cpp b/experiments/pivoting_vs_no_pivoting.cpp
index ce2b9cc..0e71dc1 100644
--- a/experiments/pivoting_vs_no_pivoting.cpp
+++ b/experiments/pivoting_vs_no_pivoting.cpp
@@ -246,15 +246,5 @@ int main() {
     exp_amplified_multiplier();
     exp_permutation();
 
-    std::cout << "\n" << std::string(60, '=') << "\n";
-    std::cout << "  Conclusion\n";
-    std::cout << std::string(60, '=') << "\n";
-    std::cout <<
-        "Partial pivoting keeps multipliers bounded by 1 in magnitude.\n"
-        "Without it, a near-zero pivot inflates multipliers and destroys\n"
-        "accuracy via catastrophic cancellation.  The 'epsilon pathology'\n"
-        "case is the textbook example: a pivot of 1e-15 makes no-pivot LU\n"
-        "compute x ≈ [0, 0.5] instead of the exact [1, 1].\n\n";
-
     return 0;
 }
diff --git a/include/qr.hpp b/include/qr.hpp
new file mode 100644
index 0000000..4d6b4ec
--- /dev/null
+++ b/include/qr.hpp
@@ -0,0 +1,45 @@
+#pragma once
+
+#include "matrix.hpp"
+#include "vector.hpp"
+
+namespace linalg {
+
+// Thin QR factorization result.
+//
+// For an m x n matrix A with m >= n:
+//   Q  is m x n with orthonormal columns  (Q^T Q = I_n)
+//   R  is n x n upper triangular
+//   A  = Q * R
+struct QRResult {
+    Matrix Q;
+    Matrix R;
+};
+
+// Classical Gram-Schmidt.
+// Mathematically natural but numerically fragile: orthogonality of Q
+// degrades rapidly on ill-conditioned inputs.
+// Provided for comparison — prefer modified_gs or householder in practice.
+//
+// Throws DimensionMismatchError if rows < cols.
+// Throws SingularMatrixError   if a column is (nearly) linearly dependent.
+QRResult qr_classical_gs(const Matrix& A, double zero_tolerance = 1e-14);
+
+// Modified Gram-Schmidt.
+// Subtracts each projection immediately on the running vector rather than
+// on the original column.  Algebraically equivalent to classical GS but
+// numerically much better — round-off stays local instead of accumulating.
+//
+// Same exceptions as classical GS.
+QRResult qr_modified_gs(const Matrix& A, double zero_tolerance = 1e-14);
+
+// Householder QR.
+// Applies a sequence of orthogonal reflections to zero out below-diagonal
+// entries column by column.  Backward-stable and the standard choice for
+// dense QR.  Works correctly on rank-deficient matrices (zero pivots
+// produce zero diagonal entries in R without throwing).
+//
+// Throws DimensionMismatchError if rows < cols.
+QRResult qr_householder(const Matrix& A);
+
+}  // namespace linalg
diff --git a/src/qr.cpp b/src/qr.cpp
new file mode 100644
index 0000000..57e4c2b
--- /dev/null
+++ b/src/qr.cpp
@@ -0,0 +1,218 @@
+#include "qr.hpp"
+
+#include <cmath>
+#include <sstream>
+
+#include "linalg_error.hpp"
+
+namespace linalg {
+
+namespace {
+
+void require_tall(const Matrix& A, const char* name) {
+    if (A.rows() < A.cols()) {
+        std::ostringstream oss;
+        oss << name << " requires rows >= cols, got " << A.rows() << "x" << A.cols();
+        throw DimensionMismatchError(oss.str());
+    }
+}
+
+// ||column j of M||_2
+double col_norm(const Matrix& M, std::size_t j) {
+    double s = 0.0;
+    for (std::size_t i = 0; i < M.rows(); ++i) {
+        s += M(i, j) * M(i, j);
+    }
+    return std::sqrt(s);
+}
+
+// dot product of column j of M with column k of N (same number of rows)
+double col_dot(const Matrix& M, std::size_t j, const Matrix& N, std::size_t k) {
+    double s = 0.0;
+    for (std::size_t i = 0; i < M.rows(); ++i) {
+        s += M(i, j) * N(i, k);
+    }
+    return s;
+}
+
+}  // namespace
+
+// ---------------------------------------------------------------------------
+// Classical Gram-Schmidt
+// ---------------------------------------------------------------------------
+//
+// For column j:
+//   R[i][j] = <a_j, q_i>   for i < j   
+//   v       = a_j - sum_i R[i][j] * q_i
+//   R[j][j] = ||v||
+//   q_j     = v / R[j][j]
+//
+
+QRResult qr_classical_gs(const Matrix& A, double zero_tolerance) {
+    require_tall(A, "qr_classical_gs");
+    const std::size_t m = A.rows();
+    const std::size_t n = A.cols();
+
+    Matrix Q = Matrix::zeros(m, n);
+    Matrix R = Matrix::zeros(n, n);
+
+    for (std::size_t j = 0; j < n; ++j) {
+        // Start with column j of A.
+        for (std::size_t i = 0; i < m; ++i) Q(i, j) = A(i, j);
+
+        // Project out existing basis vectors using the *original* A column.
+        for (std::size_t k = 0; k < j; ++k) {
+            R(k, j) = col_dot(A, j, Q, k);       // <a_j, q_k>
+            for (std::size_t i = 0; i < m; ++i) {
+                Q(i, j) -= R(k, j) * Q(i, k);
+            }
+        }
+
+        const double norm = col_norm(Q, j);
+        if (norm <= zero_tolerance) {
+            std::ostringstream oss;
+            oss << "qr_classical_gs: column " << j
+                << " is (nearly) linearly dependent (norm = " << norm << ")";
+            throw SingularMatrixError(oss.str());
+        }
+        R(j, j) = norm;
+        for (std::size_t i = 0; i < m; ++i) Q(i, j) /= norm;
+    }
+
+    return QRResult{std::move(Q), std::move(R)};
+}
+
+// ---------------------------------------------------------------------------
+// Modified Gram-Schmidt
+// ---------------------------------------------------------------------------
+//
+// For column j:
+//   v = a_j
+//   For k = 0 .. j-1:
+//     R[k][j] = <v, q_k> 
+//     v       = v - R[k][j] * q_k
+//   R[j][j] = ||v||
+//   q_j     = v / R[j][j]
+//
+// Each subtraction uses the already-updated v, so round-off is re-corrected
+// at every sub-step rather than compounding into one subtraction.
+
+QRResult qr_modified_gs(const Matrix& A, double zero_tolerance) {
+    require_tall(A, "qr_modified_gs");
+    const std::size_t m = A.rows();
+    const std::size_t n = A.cols();
+
+    Matrix Q = Matrix::zeros(m, n);
+    Matrix R = Matrix::zeros(n, n);
+
+    for (std::size_t j = 0; j < n; ++j) {
+        for (std::size_t i = 0; i < m; ++i) Q(i, j) = A(i, j);
+
+        for (std::size_t k = 0; k < j; ++k) {
+            R(k, j) = col_dot(Q, j, Q, k);       // <v_running, q_k>
+            for (std::size_t i = 0; i < m; ++i) {
+                Q(i, j) -= R(k, j) * Q(i, k);
+            }
+        }
+
+        const double norm = col_norm(Q, j);
+        if (norm <= zero_tolerance) {
+            std::ostringstream oss;
+            oss << "qr_modified_gs: column " << j
+                << " is (nearly) linearly dependent (norm = " << norm << ")";
+            throw SingularMatrixError(oss.str());
+        }
+        R(j, j) = norm;
+        for (std::size_t i = 0; i < m; ++i) Q(i, j) /= norm;
+    }
+
+    return QRResult{std::move(Q), std::move(R)};
+}
+
+// ---------------------------------------------------------------------------
+// Householder QR
+// ---------------------------------------------------------------------------
+//
+// At step k, build a Householder reflector H_k that maps R[k:, k] to
+// -sign(R[k,k]) * ||R[k:,k]|| * e_1.
+//
+// H = I - (2 / (u^T u)) * u * u^T
+//   where u = x + sign(x_0) * ||x|| * e_1   (sign chosen to avoid cancellation)
+//
+// H is never formed explicitly.  It is applied via the rank-1 update:
+//   M[k:, :] -= u * (2/(u^T u) * (u^T M[k:, :]))
+//
+// After n reflections, the working copy of A has become R (upper triangular).
+// Q is accumulated by applying each H_k to an identity matrix from the left.
+// The thin Q (m x n) is the first n columns of the full m x m orthogonal Q.
+
+QRResult qr_householder(const Matrix& A) {
+    require_tall(A, "qr_householder");
+    const std::size_t m = A.rows();
+    const std::size_t n = A.cols();
+
+    // Will become R.
+    Matrix work = A;
+
+    // Q accumulated as full m x m orthogonal matrix; trim to m x n.
+    Matrix Q_full = Matrix::identity(m);
+
+    for (std::size_t k = 0; k < n; ++k) {
+        const std::size_t p = m - k;  // length of the subvector
+
+        // Build Householder vector u from the subcolumn work[k:, k].
+        std::vector<double> u(p);
+        for (std::size_t i = 0; i < p; ++i) u[i] = work(k + i, k);
+
+        const double x_norm = [&] {
+            double s = 0.0;
+            for (double v : u) s += v * v;
+            return std::sqrt(s);
+        }();
+
+        if (x_norm == 0.0) continue;
+
+        // sigma = sign(u[0]) * ||x||
+        const double sigma = (u[0] >= 0.0 ? 1.0 : -1.0) * x_norm;
+        u[0] += sigma;
+
+        const double utu = [&] {
+            double s = 0.0;
+            for (double v : u) s += v * v;
+            return s;
+        }();
+        const double tau = 2.0 / utu;
+
+        // Apply H_k to work[k:, k:n] 
+        for (std::size_t j = k; j < n; ++j) {
+            double dot = 0.0;
+            for (std::size_t i = 0; i < p; ++i) dot += u[i] * work(k + i, j);
+            const double coeff = tau * dot;
+            for (std::size_t i = 0; i < p; ++i) work(k + i, j) -= coeff * u[i];
+        }
+
+        // Apply H_k to Q_full[k:, 0:m] 
+        for (std::size_t j = 0; j < m; ++j) {
+            double dot = 0.0;
+            for (std::size_t i = 0; i < p; ++i) dot += u[i] * Q_full(k + i, j);
+            const double coeff = tau * dot;
+            for (std::size_t i = 0; i < p; ++i) Q_full(k + i, j) -= coeff * u[i];
+        }
+    }
+
+    // Thin Q: first n columns of Q_full^T
+    Matrix Q(m, n);
+    for (std::size_t i = 0; i < m; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            Q(i, j) = Q_full(j, i);
+
+    // Thin R: first n rows of work
+    Matrix R(n, n);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            R(i, j) = work(i, j);
+
+    return QRResult{std::move(Q), std::move(R)};
+}
+
+}  // namespace linalg
diff --git a/tests/test_qr.cpp b/tests/test_qr.cpp
new file mode 100644
index 0000000..a837132
--- /dev/null
+++ b/tests/test_qr.cpp
@@ -0,0 +1,280 @@
+#include "linalg_error.hpp"
+#include "matrix.hpp"
+#include "norms.hpp"
+#include "qr.hpp"
+#include "vector.hpp"
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <functional>
+#include <random>
+
+using linalg::DimensionMismatchError;
+using linalg::Matrix;
+using linalg::QRResult;
+using linalg::SingularMatrixError;
+
+// ---------------------------------------------------------------------------
+// Helpers
+// ---------------------------------------------------------------------------
+
+namespace {
+
+// ||A - QR||_F
+double reconstruction_error(const Matrix& A, const QRResult& qr) {
+    const Matrix diff = A - qr.Q * qr.R;
+    double err = 0.0;
+    for (std::size_t i = 0; i < diff.rows(); ++i)
+        for (std::size_t j = 0; j < diff.cols(); ++j)
+            err += diff(i, j) * diff(i, j);
+    return std::sqrt(err);
+}
+
+// ||Q^T Q - I||_F  (should be ~0 for orthonormal Q)
+double orthogonality_error(const QRResult& qr) {
+    const Matrix& Q = qr.Q;
+    const std::size_t n = Q.cols();
+    // Compute Q^T Q
+    Matrix QtQ(n, n);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j) {
+            double s = 0.0;
+            for (std::size_t k = 0; k < Q.rows(); ++k) s += Q(k, i) * Q(k, j);
+            QtQ(i, j) = s;
+        }
+    // ||QtQ - I||_F
+    double err = 0.0;
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j) {
+            const double d = QtQ(i, j) - (i == j ? 1.0 : 0.0);
+            err += d * d;
+        }
+    return std::sqrt(err);
+}
+
+// R must be upper triangular (strict lower triangle near zero).
+bool r_is_upper_triangular(const Matrix& R, double tol = 1e-12) {
+    for (std::size_t i = 1; i < R.rows(); ++i)
+        for (std::size_t j = 0; j < i; ++j)
+            if (std::abs(R(i, j)) > tol) return false;
+    return true;
+}
+
+Matrix random_matrix(std::size_t m, std::size_t n, unsigned seed = 42) {
+    std::mt19937 rng(seed);
+    std::uniform_real_distribution<double> dist(-5.0, 5.0);
+    Matrix M(m, n);
+    for (std::size_t i = 0; i < m; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            M(i, j) = dist(rng);
+    return M;
+}
+
+// Run all checks for a given QR function and matrix.
+using QRFn = std::function<QRResult(const Matrix&)>;
+
+void check_qr(const Matrix& A, QRFn fn,
+              double recon_tol, double ortho_tol,
+              const std::string& /*label*/) {
+    const QRResult qr = fn(A);
+    CHECK(qr.Q.rows() == A.rows());
+    CHECK(qr.Q.cols() == A.cols());
+    CHECK(qr.R.rows() == A.cols());
+    CHECK(qr.R.cols() == A.cols());
+    CHECK(reconstruction_error(A, qr) == Catch::Approx(0.0).margin(recon_tol));
+    CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(ortho_tol));
+    CHECK(r_is_upper_triangular(qr.R));
+}
+
+}  // namespace
+
+// ---------------------------------------------------------------------------
+// Macro to run the same test body for all three methods
+// ---------------------------------------------------------------------------
+
+#define FOR_ALL_METHODS(A, recon_tol, ortho_tol)                                  \
+    SECTION("classical_gs") {                                                      \
+        check_qr(A, [](const Matrix& M) { return linalg::qr_classical_gs(M); },  \
+                 recon_tol, ortho_tol, "classical_gs");                            \
+    }                                                                              \
+    SECTION("modified_gs") {                                                       \
+        check_qr(A, [](const Matrix& M) { return linalg::qr_modified_gs(M); },   \
+                 recon_tol, ortho_tol, "modified_gs");                             \
+    }                                                                              \
+    SECTION("householder") {                                                       \
+        check_qr(A, [](const Matrix& M) { return linalg::qr_householder(M); },   \
+                 recon_tol, ortho_tol, "householder");                             \
+    }
+
+// ---------------------------------------------------------------------------
+// Basic correctness: square matrices
+// ---------------------------------------------------------------------------
+
+TEST_CASE("QR: 3x3 known matrix", "[qr]") {
+    const Matrix A{
+        {1.0, 2.0, 3.0},
+        {4.0, 5.0, 6.0},
+        {7.0, 8.0, 10.0}   // not exactly singular
+    };
+    FOR_ALL_METHODS(A, 1e-12, 1e-12)
+}
+
+TEST_CASE("QR: identity matrix", "[qr]") {
+    const Matrix I = Matrix::identity(4);
+    FOR_ALL_METHODS(I, 1e-14, 1e-14)
+}
+
+TEST_CASE("QR: diagonal matrix", "[qr]") {
+    const Matrix D{
+        {3.0, 0.0, 0.0},
+        {0.0, 1.0, 0.0},
+        {0.0, 0.0, 2.0}
+    };
+    FOR_ALL_METHODS(D, 1e-14, 1e-14)
+}
+
+// ---------------------------------------------------------------------------
+// Rectangular (tall) matrices
+// ---------------------------------------------------------------------------
+
+TEST_CASE("QR: tall 5x3 random matrix", "[qr]") {
+    const Matrix A = random_matrix(5, 3, 7u);
+    FOR_ALL_METHODS(A, 1e-12, 1e-12)
+}
+
+TEST_CASE("QR: tall 10x4 random matrix", "[qr]") {
+    const Matrix A = random_matrix(10, 4, 99u);
+    FOR_ALL_METHODS(A, 1e-12, 1e-12)
+}
+
+// ---------------------------------------------------------------------------
+// Random square matrices
+// ---------------------------------------------------------------------------
+
+TEST_CASE("QR: random 6x6", "[qr]") {
+    const Matrix A = random_matrix(6, 6, 123u);
+    FOR_ALL_METHODS(A, 1e-12, 1e-12)
+}
+
+TEST_CASE("QR: random 12x12", "[qr]") {
+    const Matrix A = random_matrix(12, 12, 456u);
+    FOR_ALL_METHODS(A, 1e-11, 1e-11)
+}
+
+// ---------------------------------------------------------------------------
+// Nearly dependent columns — GS methods degrade; Householder stays clean
+// ---------------------------------------------------------------------------
+
+TEST_CASE("QR: nearly dependent columns", "[qr]") {
+    // Column 1 = column 0 + epsilon * e_1.
+    // Classical GS will lose most of Q's orthogonality here.
+    // Modified GS is better. Householder is unaffected.
+    constexpr double eps = 1e-7;
+    const Matrix A{
+        {1.0,      1.0 + eps, 0.0},
+        {1.0,      1.0,       1.0},
+        {0.0,      eps,       1.0},
+        {0.0,      0.0,       1.0}
+    };
+
+    // All three should reconstruct A accurately.
+    SECTION("classical_gs reconstruction") {
+        const QRResult qr = linalg::qr_classical_gs(A);
+        CHECK(reconstruction_error(A, qr) == Catch::Approx(0.0).margin(1e-10));
+        // Orthogonality will be poor for classical GS on this input.
+        // We only assert it's not catastrophically wrong (< 0.01).
+        CHECK(orthogonality_error(qr) < 0.01);
+    }
+    SECTION("modified_gs reconstruction") {
+        const QRResult qr = linalg::qr_modified_gs(A);
+        CHECK(reconstruction_error(A, qr) == Catch::Approx(0.0).margin(1e-10));
+        CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-8));
+    }
+    SECTION("householder reconstruction") {
+        const QRResult qr = linalg::qr_householder(A);
+        CHECK(reconstruction_error(A, qr) == Catch::Approx(0.0).margin(1e-13));
+        CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-13));
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Hilbert-like ill-conditioned matrix
+// ---------------------------------------------------------------------------
+
+TEST_CASE("QR: 4x4 Hilbert matrix", "[qr]") {
+    // H[i][j] = 1 / (i + j + 1)
+    const std::size_t n = 4;
+    Matrix H(n, n);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            H(i, j) = 1.0 / static_cast<double>(i + j + 1);
+
+    SECTION("classical_gs") {
+        const QRResult qr = linalg::qr_classical_gs(H);
+        CHECK(reconstruction_error(H, qr) == Catch::Approx(0.0).margin(1e-12));
+        // Orthogonality is imperfect on Hilbert matrices with classical GS.
+        CHECK(orthogonality_error(qr) < 1e-8);
+    }
+    SECTION("modified_gs") {
+        const QRResult qr = linalg::qr_modified_gs(H);
+        CHECK(reconstruction_error(H, qr) == Catch::Approx(0.0).margin(1e-12));
+        CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-10));
+    }
+    SECTION("householder") {
+        const QRResult qr = linalg::qr_householder(H);
+        CHECK(reconstruction_error(H, qr) == Catch::Approx(0.0).margin(1e-13));
+        CHECK(orthogonality_error(qr) == Catch::Approx(0.0).margin(1e-13));
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Structure checks
+// ---------------------------------------------------------------------------
+
+TEST_CASE("QR: R is upper triangular", "[qr]") {
+    const Matrix A = random_matrix(5, 5, 555u);
+    CHECK(r_is_upper_triangular(linalg::qr_classical_gs(A).R));
+    CHECK(r_is_upper_triangular(linalg::qr_modified_gs(A).R));
+    CHECK(r_is_upper_triangular(linalg::qr_householder(A).R));
+}
+
+TEST_CASE("QR: Q columns are unit length", "[qr]") {
+    const Matrix A = random_matrix(6, 4, 321u);
+    for (QRFn fn : {QRFn{[](const Matrix& M) { return linalg::qr_classical_gs(M); }},
+                    QRFn{[](const Matrix& M) { return linalg::qr_modified_gs(M); }},
+                    QRFn{[](const Matrix& M) { return linalg::qr_householder(M); }}}) {
+        const QRResult qr = fn(A);
+        for (std::size_t j = 0; j < qr.Q.cols(); ++j) {
+            double norm2 = 0.0;
+            for (std::size_t i = 0; i < qr.Q.rows(); ++i)
+                norm2 += qr.Q(i, j) * qr.Q(i, j);
+            CHECK(std::sqrt(norm2) == Catch::Approx(1.0).margin(1e-13));
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Failure cases
+// ---------------------------------------------------------------------------
+
+TEST_CASE("QR: fat matrix throws DimensionMismatchError", "[qr]") {
+    const Matrix A(3, 5);   // rows < cols
+    CHECK_THROWS_AS(linalg::qr_classical_gs(A), DimensionMismatchError);
+    CHECK_THROWS_AS(linalg::qr_modified_gs(A), DimensionMismatchError);
+    CHECK_THROWS_AS(linalg::qr_householder(A), DimensionMismatchError);
+}
+
+TEST_CASE("QR: linearly dependent columns throw from GS methods", "[qr]") {
+    const Matrix A{
+        {1.0, 2.0, 2.0},   // col 2 = 2 * col 0
+        {2.0, 4.0, 4.0},
+        {3.0, 6.0, 6.0}
+    };
+    CHECK_THROWS_AS(linalg::qr_classical_gs(A), SingularMatrixError);
+    CHECK_THROWS_AS(linalg::qr_modified_gs(A), SingularMatrixError);
+    // Householder handles rank-deficient input gracefully (R gets a zero diagonal entry).
+    CHECK_NOTHROW(linalg::qr_householder(A));
+}