From b11123e1cd710c135d9924f263ab1808cd96c8c2 Mon Sep 17 00:00:00 2001 From: y-jan137 Date: Sun, 15 Mar 2026 09:11:40 +0300 Subject: Add QR decomposition --- CMakeLists.txt | 5 + README.md | 16 +- experiments/hilbert_qr.cpp | 167 +++++++++++++++++++ experiments/pivoting_vs_no_pivoting.cpp | 10 -- include/qr.hpp | 45 +++++ src/qr.cpp | 218 +++++++++++++++++++++++++ tests/test_qr.cpp | 280 ++++++++++++++++++++++++++++++++ 7 files changed, 730 insertions(+), 11 deletions(-) create mode 100644 experiments/hilbert_qr.cpp create mode 100644 include/qr.hpp create mode 100644 src/qr.cpp create mode 100644 tests/test_qr.cpp 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 +#include +#include +#include +#include +#include +#include +#include +#include + +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(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; + +// Run fn() `trials` times, return minimum elapsed seconds. +template +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; + +struct Result { double recon, ortho, time_s; }; + +std::optional 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 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 +#include + +#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] = 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); // + 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 = 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); // + 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 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 +#include + +#include +#include +#include +#include + +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 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; + +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(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)); +} -- cgit v1.2.3