From be01607041d2ded41e5bf39a0a8e17ffb9bcc296 Mon Sep 17 00:00:00 2001 From: y-jan137 Date: Sun, 15 Mar 2026 11:16:49 +0300 Subject: Add QR iteration and cleanup --- CMakeLists.txt | 14 +- examples/linear_system.cpp | 36 ---- examples/matmul.cpp | 32 ---- experiments/hilbert_qr.cpp | 17 -- experiments/matmul.cpp | 149 ++++++++++++++++ experiments/pivoting_vs_no_pivoting.cpp | 16 -- include/linalg_error.hpp | 6 + include/qr_iteration.hpp | 82 +++++++++ src/qr_iteration.cpp | 207 ++++++++++++++++++++++ tests/test_qr_iteration.cpp | 293 ++++++++++++++++++++++++++++++++ 10 files changed, 742 insertions(+), 110 deletions(-) delete mode 100644 examples/linear_system.cpp delete mode 100644 examples/matmul.cpp create mode 100644 experiments/matmul.cpp create mode 100644 include/qr_iteration.hpp create mode 100644 src/qr_iteration.cpp create mode 100644 tests/test_qr_iteration.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 1d2374c..1f40b0c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,7 +8,6 @@ set(CMAKE_CXX_EXTENSIONS OFF) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) option(LINEAR_ALGEBRA_BUILD_TESTS "Build unit tests" ON) -option(LINEAR_ALGEBRA_BUILD_EXAMPLES "Build example programs" ON) set(LINEAR_ALGEBRA_SIMD "AUTO" CACHE STRING "SIMD backend for matmul: AUTO, NONE, AVX, AVX2, AVX512") set_property(CACHE LINEAR_ALGEBRA_SIMD PROPERTY STRINGS AUTO NONE AVX AVX2 AVX512) @@ -19,6 +18,7 @@ add_library(linear_algebra src/triangular_solve.cpp src/lu.cpp src/qr.cpp + src/qr_iteration.cpp ) add_library(linear_algebra::core ALIAS linear_algebra) @@ -54,14 +54,6 @@ elseif(MSVC) endif() endif() -if(LINEAR_ALGEBRA_BUILD_EXAMPLES) - add_executable(linear_system examples/linear_system.cpp) - target_link_libraries(linear_system PRIVATE linear_algebra::core) - - add_executable(matmul examples/matmul.cpp) - target_link_libraries(matmul PRIVATE linear_algebra::core) -endif() - 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) @@ -69,6 +61,9 @@ if(LINEAR_ALGEBRA_BUILD_EXPERIMENTS) add_executable(hilbert_qr experiments/hilbert_qr.cpp) target_link_libraries(hilbert_qr PRIVATE linear_algebra::core) + + add_executable(matmul experiments/matmul.cpp) + target_link_libraries(matmul PRIVATE linear_algebra::core) endif() if(LINEAR_ALGEBRA_BUILD_TESTS) @@ -93,6 +88,7 @@ if(LINEAR_ALGEBRA_BUILD_TESTS) tests/test_triangular_solve.cpp tests/test_lu.cpp tests/test_qr.cpp + tests/test_qr_iteration.cpp ) target_link_libraries(linear_algebra_tests diff --git a/examples/linear_system.cpp b/examples/linear_system.cpp deleted file mode 100644 index 184f576..0000000 --- a/examples/linear_system.cpp +++ /dev/null @@ -1,36 +0,0 @@ -#include "matrix.hpp" -#include "norms.hpp" -#include "triangular_solve.hpp" -#include "vector.hpp" - -#include -#include - -int main() { - const linalg::Matrix basis = linalg::Matrix::identity(3); - const linalg::Matrix upper{ - {4.0, -2.0, 1.0}, - {0.0, 3.0, 5.0}, - {0.0, 0.0, -2.0} - }; - const linalg::Vector expected{2.0, -1.0, 3.0}; - const linalg::Vector rhs = upper * expected; - const linalg::Vector x = linalg::backward_substitution(upper, rhs); - - std::cout << "Week 3 triangular solve demo\n"; - std::cout << "Identity matrix diagonal: "; - for (std::size_t i = 0; i < basis.rows(); ++i) { - std::cout << basis(i, i) << (i + 1 == basis.rows() ? '\n' : ' '); - } - - std::cout << "Recovered solution x: "; - for (std::size_t i = 0; i < x.size(); ++i) { - std::cout << std::fixed << std::setprecision(2) << x[i] - << (i + 1 == x.size() ? '\n' : ' '); - } - - const linalg::Vector residual = (upper * x) - rhs; - std::cout << "Residual 2-norm = " << linalg::norm2(residual) << '\n'; - - return 0; -} diff --git a/examples/matmul.cpp b/examples/matmul.cpp deleted file mode 100644 index c6245f3..0000000 --- a/examples/matmul.cpp +++ /dev/null @@ -1,32 +0,0 @@ -#include "matrix.hpp" - -#include -#include - -int main() { - const linalg::Matrix a{ - {1.0, 2.0, 3.0}, - {4.0, 5.0, 6.0} - }; - const linalg::Matrix b{ - {7.0, 8.0}, - {9.0, 10.0}, - {11.0, 12.0} - }; - - const linalg::Matrix c = a * b; - - std::cout << "Matrix multiplication\n"; - std::cout << "A is " << a.rows() << " x " << a.cols() << '\n'; - std::cout << "B is " << b.rows() << " x " << b.cols() << '\n'; - std::cout << "C = A * B:\n"; - - for (std::size_t i = 0; i < c.rows(); ++i) { - for (std::size_t j = 0; j < c.cols(); ++j) { - std::cout << std::fixed << std::setprecision(2) << c(i, j) - << (j + 1 == c.cols() ? '\n' : ' '); - } - } - - return 0; -} diff --git a/experiments/hilbert_qr.cpp b/experiments/hilbert_qr.cpp index 9209a30..218b35c 100644 --- a/experiments/hilbert_qr.cpp +++ b/experiments/hilbert_qr.cpp @@ -1,9 +1,3 @@ -// 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 @@ -68,9 +62,6 @@ double orthogonality_error(const QRResult& qr) { return std::sqrt(err); } -// --------------------------------------------------------------------------- -// Timing -// --------------------------------------------------------------------------- using Clock = std::chrono::high_resolution_clock; using Seconds = std::chrono::duration; @@ -88,10 +79,6 @@ double min_time(Fn fn, int trials = 5) { 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; }; @@ -110,10 +97,6 @@ std::optional measure(const Matrix& A, QRFn fn) { } } -// --------------------------------------------------------------------------- -// Pretty printing -// --------------------------------------------------------------------------- - void print_row(const std::string& method, std::optional r) { std::cout << std::left << std::setw(16) << method; if (!r) { diff --git a/experiments/matmul.cpp b/experiments/matmul.cpp new file mode 100644 index 0000000..8ef0ac5 --- /dev/null +++ b/experiments/matmul.cpp @@ -0,0 +1,149 @@ +#include "linalg_error.hpp" +#include "matrix.hpp" + +#include +#include +#include +#include +#include + +#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX512F__) +# define MATMUL_BACKEND "AVX512" +#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX2__) +# define MATMUL_BACKEND "AVX2" +#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX__) +# define MATMUL_BACKEND "AVX" +#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && \ + defined(__ARM_NEON) && defined(__aarch64__) && \ + defined(__ARM_FEATURE_FP64_VECTOR_ARITHMETIC) +# define MATMUL_BACKEND "NEON" +#else +# define MATMUL_BACKEND "scalar" +#endif + +using linalg::Matrix; +using Clock = std::chrono::high_resolution_clock; +using Seconds = std::chrono::duration; + + +Matrix naive_matmul(const Matrix& lhs, const Matrix& rhs) { + if (lhs.cols() != rhs.rows()) { + throw linalg::DimensionMismatchError( + "naive_matmul: lhs.cols() != rhs.rows()"); + } + const std::size_t m = lhs.rows(); + const std::size_t n = rhs.cols(); + const std::size_t k = lhs.cols(); + + Matrix result(m, n, 0.0); + + const double* A = lhs.data(); + const double* B = rhs.data(); + double* C = result.data(); + + for (std::size_t i = 0; i < m; ++i) + for (std::size_t j = 0; j < n; ++j) + for (std::size_t p = 0; p < k; ++p) + C[i * n + j] += A[i * k + p] * B[p * n + j]; + + return result; +} + +// --------------------------------------------------------------------------- +// Helpers +// --------------------------------------------------------------------------- + +// Fill an n×n matrix with a deterministic pattern so the compiler cannot +// optimise multiplications away. +Matrix make_matrix(std::size_t n) { + Matrix M(n, n); + const double inv = 1.0 / static_cast(n + 1); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + M(i, j) = static_cast(i + j + 1) * inv; + return M; +} + +template +double min_time_s(Fn fn, int trials) { + double best = 1e30; + for (int t = 0; t < trials; ++t) { + const auto t0 = Clock::now(); + fn(); + const auto t1 = Clock::now(); + const double elapsed = Seconds(t1 - t0).count(); + if (elapsed < best) best = elapsed; + } + return best; +} + +double flops(std::size_t n) { + const double nd = static_cast(n); + return 2.0 * nd * nd * nd; +} + + +int main() { + const std::vector sizes = { + 8, 16, 32, 64, 128, 256, 512 + }; + + constexpr std::size_t small_threshold = 128; + constexpr int trials_small = 9; + constexpr int trials_large = 3; + + std::cout << std::string(72, '*') << "\n"; + std::cout << " Matmul benchmark: naive (ijk) vs SIMD (" MATMUL_BACKEND + ") + transpose\n"; + std::cout << " C = A * B, A and B both n×n\n"; + std::cout << std::string(72, '*') << "\n\n"; + + // Column header. + std::cout << std::left + << std::setw(6) << "n" + << std::setw(14) << "naive ms" + << std::setw(14) << "naive GFLOP/s" + << std::setw(14) << "SIMD ms" + << std::setw(14) << "SIMD GFLOP/s" + << std::setw(10) << "speedup" + << "\n"; + std::cout << std::string(72, '-') << "\n"; + + for (const std::size_t n : sizes) { + const Matrix A = make_matrix(n); + const Matrix B = make_matrix(n); + + const int trials = (n <= small_threshold) ? trials_small : trials_large; + + volatile double sink_naive = naive_matmul(A, B)(0, 0); + volatile double sink_simd = (A * B)(0, 0); + (void)sink_naive; + (void)sink_simd; + + const double t_naive = min_time_s([&] { (void)naive_matmul(A, B); }, trials); + const double t_simd = min_time_s([&] { (void)(A * B); }, trials); + + const double fp = flops(n); + const double gf_naive = fp / t_naive / 1e9; + const double gf_simd = fp / t_simd / 1e9; + const double speedup = t_naive / t_simd; + + std::cout << std::left << std::setw(6) << n + << std::fixed << std::setprecision(3) + << std::setw(14) << t_naive * 1e3 + << std::setprecision(2) + << std::setw(14) << gf_naive + << std::setprecision(3) + << std::setw(14) << t_simd * 1e3 + << std::setprecision(2) + << std::setw(14) << gf_simd + << std::setprecision(2) << std::setw(10) << speedup + << "x\n"; + } + + std::cout << "\n(each cell = minimum over " + << trials_small << " trials for n<=" << small_threshold + << ", " << trials_large << " trials for larger n)\n"; + + return 0; +} diff --git a/experiments/pivoting_vs_no_pivoting.cpp b/experiments/pivoting_vs_no_pivoting.cpp index 0e71dc1..57281d8 100644 --- a/experiments/pivoting_vs_no_pivoting.cpp +++ b/experiments/pivoting_vs_no_pivoting.cpp @@ -1,8 +1,3 @@ -// Experiment: partial pivoting vs no-pivot LU -// -// Demonstrates why partial pivoting is essential for numerical stability. -// Run the binary and inspect the residuals printed to stdout. - #include "lu.hpp" #include "matrix.hpp" #include "norms.hpp" @@ -23,7 +18,6 @@ using linalg::Vector; // --------------------------------------------------------------------------- // Local no-pivot LU for comparison only. -// This is intentionally naive — it is here to show what breaks without pivoting. // --------------------------------------------------------------------------- struct NoPivotLU { @@ -55,8 +49,6 @@ NoPivotLU lu_no_pivot(const Matrix& A, double tol = 1e-14) { return NoPivotLU{std::move(L), std::move(U), false, 0}; } -// Solve using a no-pivot LU (L unit lower triangular, U upper triangular). -// If the factorization failed or U is numerically singular, returns nullopt. std::optional solve_no_pivot(const NoPivotLU& f, const Vector& b) { if (f.failed) return std::nullopt; try { @@ -178,8 +170,6 @@ void exp_random(std::size_t n = 8) { } // 2. Badly row-scaled matrix -// Rows differ in magnitude by ~10^14. Without pivoting, tiny early pivots -// amplify round-off; with pivoting, the large-row is selected first. void exp_badly_scaled() { const Matrix A{ {1e-14, 1.0, 2.0 }, @@ -191,21 +181,15 @@ void exp_badly_scaled() { } // 3. Classic pathological example for no-pivot LU. -// With epsilon = 1e-15, no-pivot computes a huge multiplier (1/epsilon), -// which causes catastrophic cancellation in the updated rows. -// With pivoting, we swap first and the multiplier is bounded by 1. void exp_epsilon_pathology() { constexpr double eps = 1e-15; const Matrix A{{eps, 1.0}, {1.0, 2.0}}; - // True solution of [eps 1; 1 2] * x = [1+eps; 3] is x = [1; 1]. const Vector b{1.0 + eps, 3.0}; run_case("Epsilon pathology [[1e-15,1],[1,2]] (classic)", A, b); std::cout << " Note: exact solution is x = [1, 1]\n"; } // 4. Matrix where no-pivot LU diverges visibly on a 4x4 example. -// The first pivot is small (0.001) but rows below have entries ~1000. -// No pivot causes multipliers of magnitude 10^6, annihilating subdiagonal info. void exp_amplified_multiplier() { const Matrix A{ {0.001, 1.0, 0.0, 0.0 }, diff --git a/include/linalg_error.hpp b/include/linalg_error.hpp index 25c48ad..ac9b456 100644 --- a/include/linalg_error.hpp +++ b/include/linalg_error.hpp @@ -22,4 +22,10 @@ public: : LinAlgError(message) {} }; +class NonConvergenceError : public LinAlgError { +public: + explicit NonConvergenceError(const std::string& message) + : LinAlgError(message) {} +}; + } // namespace linalg diff --git a/include/qr_iteration.hpp b/include/qr_iteration.hpp new file mode 100644 index 0000000..08e17fa --- /dev/null +++ b/include/qr_iteration.hpp @@ -0,0 +1,82 @@ +#pragma once + +// QR iteration for eigenvalue computation +// +// Refs: +// Trefethen & Bau, "Numerical Linear Algebra" (T&B) +// Lecture 25 — Eigenvalue algorithms +// Lecture 26 — Schur factorisation +// Lecture 28 — The QR algorithm (unshifted) +// Lecture 29 — The QR algorithm with shifts +// Golub & Van Loan, "Matrix Computations" 4th ed. (GVL) +// §7.3 — The Unshifted QR Algorithm +// §7.4 — The Shifted QR Algorithm +// §7.4.2 — Wilkinson shift + +#include + +#include "linalg_error.hpp" +#include "matrix.hpp" +#include "vector.hpp" + +namespace linalg { + +// --------------------------------------------------------------------------- +// Options +// --------------------------------------------------------------------------- + +// All defaults are consistent with the recommendations in T&B Lecture 28. +struct QRIterationOptions { + // Convergence threshold. Iteration halts once the Frobenius norm of the + // strict lower triangle of A_k falls below this value. + // Ref: T&B §28; GVL §7.3. + double tolerance = 1e-10; + + int max_iterations = 1000; + + // When true, the Frobenius norm of the strict lower triangle is recorded + // after every QR step and returned in QRIterationResult::convergence_history. + bool track_convergence = false; +}; + + +struct QRIterationResult { + // Real and imaginary parts of the n eigenvalues. + // For symmetric inputs all imaginary parts are zero. + // Complex-conjugate pairs from 2×2 Schur blocks appear as ±imag entries. + // Both vectors always have length n (the matrix dimension). + Vector eigenvalues_real; + Vector eigenvalues_imag; + + // Total number of QR steps performed before convergence or max_iterations. + int iterations = 0; + + // Populated only when QRIterationOptions::track_convergence is true. + // Entry k is ||lower(A_k)||_F after the k-th QR step. + // std::vector is used here because linalg::Vector has no push_back; + // convergence_history is a plain time-series container, not a math object. + std::vector convergence_history; +}; + +// Algorithm (T&B Algorithm 28.1): +// +// A_0 = A +// for k = 1, 2, ...: +// factor A_{k-1} = Q_k R_k (Householder QR) +// set A_k = R_k Q_k (orthogonal similarity: preserves eigenvalues) +// +// The iterates A_k converge to the real Schur form of A: a quasi-upper- +// triangular matrix with 1×1 blocks (real eigenvalue) and 2×2 blocks +// (complex-conjugate pair) on the diagonal. +// +// Convergence rate: linear. Per-step factor ≈ |lambda_{j+1} / lambda_j| +// for the off-diagonal entries linking eigenvalue clusters j and j+1. +// (T&B Lecture 28, Theorem 28.2) +// +// Throws DimensionMismatchError if A is not square. +// Throws NonConvergenceError if convergence is not achieved within +// opts.max_iterations steps. +[[nodiscard]] QRIterationResult eigenvalues_unshifted(const Matrix& A, + QRIterationOptions opts = {}); + +} // namespace linalg diff --git a/src/qr_iteration.cpp b/src/qr_iteration.cpp new file mode 100644 index 0000000..f488a5f --- /dev/null +++ b/src/qr_iteration.cpp @@ -0,0 +1,207 @@ +#include "qr_iteration.hpp" + +#include +#include +#include + +#include "linalg_error.hpp" +#include "matrix.hpp" +#include "qr.hpp" +#include "vector.hpp" + +// References used throughout this file: +// T&B — Trefethen & Bau, "Numerical Linear Algebra" +// GVL — Golub & Van Loan, "Matrix Computations" 4th ed. + +namespace linalg { + +namespace { + +// --------------------------------------------------------------------------- +// Internal helpers +// --------------------------------------------------------------------------- + +// Frobenius norm of the strict lower triangle of an n×n matrix. +// This is the standard convergence diagnostic for QR iteration: as A_k +// approaches the real Schur form, all entries below the main diagonal +// (excluding 2×2 block sub-diagonals) tend to zero. +// +// ||lower(A)||_F = sqrt( sum_{i > j} A(i,j)^2 ) +// +// Ref: T&B §28; used as the convergence criterion in Algorithm 28.1. +double lower_triangle_norm(const Matrix& A) { + const std::size_t n = A.rows(); + double s = 0.0; + for (std::size_t i = 1; i < n; ++i) // row 1 .. n-1 + for (std::size_t j = 0; j < i; ++j) // col 0 .. i-1 (strict lower) + s += A(i, j) * A(i, j); + return std::sqrt(s); +} + +// Extract eigenvalues from a quasi-upper-triangular matrix (real Schur form). +// +// Scans the diagonal from top-left to bottom-right. At each position i: +// — |A(i+1, i)| < tol → 1×1 block: real eigenvalue A(i,i), imag = 0. +// — otherwise → 2×2 block [A(i..i+1, i..i+1)]: eigenvalues via +// quadratic formula. When the discriminant is +// negative the result is a complex-conjugate pair, +// stored as (re, +im) and (re, -im) in the real +// and imaginary part Vectors. +// +// Fills positions 0..n-1 of `real_out` and `imag_out` (pre-sized to n). +// +// Ref: T&B Lecture 28; GVL §7.4.1. +void extract_eigenvalues(const Matrix& T, double tol, + Vector& real_out, Vector& imag_out) { + const std::size_t n = T.rows(); + std::size_t out = 0; // next write position in real_out / imag_out + std::size_t i = 0; // current scan position in T + + while (i < n) { + const bool is_last = (i + 1 == n); + const bool sub_small = is_last || (std::abs(T(i + 1, i)) < tol); + + if (sub_small) { + // 1×1 block: real eigenvalue. + real_out[out] = T(i, i); + imag_out[out] = 0.0; + ++out; + ++i; + } else { + // 2×2 block: + // | a b | + // | c d | + // Characteristic polynomial: lambda^2 - (a+d)*lambda + (ad - bc) = 0. + // Discriminant: (a-d)^2 + 4*b*c. + // Ref: GVL §7.4.1. + const double a = T(i, i); + const double b = T(i, i + 1); + const double c = T(i + 1, i); + const double d = T(i + 1, i + 1); + const double tr = a + d; + const double disc = (a - d) * (a - d) + 4.0 * b * c; + + if (disc >= 0.0) { + // Real eigenvalues — unusual in converged real Schur form, but + // handled robustly in case the block didn't fully split. + const double sq = std::sqrt(disc); + real_out[out] = 0.5 * (tr + sq); + imag_out[out] = 0.0; + real_out[out + 1] = 0.5 * (tr - sq); + imag_out[out + 1] = 0.0; + } else { + // Complex-conjugate pair: real part ± imaginary part. + const double re = 0.5 * tr; + const double im = 0.5 * std::sqrt(-disc); + real_out[out] = re; + imag_out[out] = im; + real_out[out + 1] = re; + imag_out[out + 1] = -im; + } + out += 2; + i += 2; + } + } + + assert(out == n); +} + +// Verify that A is square; throw DimensionMismatchError otherwise. +void require_square(const Matrix& A, const char* fname) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << fname << ": requires a square matrix, got " + << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } +} + +} // namespace + +// --------------------------------------------------------------------------- +// Stage 1: Unshifted QR iteration +// --------------------------------------------------------------------------- +// +// Each step performs an orthogonal similarity transformation: +// A_{k-1} = Q_k R_k (Householder QR; backward-stable) +// A_k = R_k Q_k = Q_k^T A_{k-1} Q_k +// +// Similarity preserves eigenvalues (GVL §7.3.1, Theorem 7.3.1). +// The iterates converge to the real Schur form: a quasi-upper-triangular +// matrix whose 1×1 blocks give real eigenvalues and 2×2 blocks give +// complex-conjugate pairs. +// +// Convergence rate: linear, with per-step reduction factor +// |lambda_{j+1} / lambda_j| for the (j, j+1) coupling. +// (T&B Lecture 28, Theorem 28.2; GVL §7.3.2) +// +// Each iteration costs O(n^3) due to full Householder QR; Hessenberg +// reduction (Stage 3) reduces subsequent steps to O(n^2). + +QRIterationResult eigenvalues_unshifted(const Matrix& A, + QRIterationOptions opts) { + require_square(A, "eigenvalues_unshifted"); + const std::size_t n = A.rows(); + + // Threshold for classifying a sub-diagonal entry as "zero" when reading + // eigenvalues out of the converged Schur form. Using the same value as + // the convergence tolerance is appropriate; we only reach extraction once + // ||lower(A_k)||_F < opts.tolerance. Ref: GVL §7.4.1. + const double extract_tol = opts.tolerance; + + QRIterationResult result; + // Pre-size eigenvalue Vectors; they are always length n. + result.eigenvalues_real = Vector(n, 0.0); + result.eigenvalues_imag = Vector(n, 0.0); + + if (opts.track_convergence) { + result.convergence_history.reserve( + static_cast(opts.max_iterations)); + } + + // Handle the trivial 1×1 case immediately. + if (n == 1) { + result.eigenvalues_real[0] = A(0, 0); + return result; + } + + // Working copy; becomes the quasi-upper-triangular Schur form A_k. + Matrix Ak = A; + + for (int k = 0; k < opts.max_iterations; ++k) { + // --- QR step --- + // Factor A_{k-1} = Q R using backward-stable Householder reflections. + const QRResult qr = qr_householder(Ak); + + // A_k = R Q (orthogonal similarity: Q^T A_{k-1} Q) + Ak = qr.R * qr.Q; + + // --- Convergence check --- + const double lower_norm = lower_triangle_norm(Ak); + + if (opts.track_convergence) { + result.convergence_history.push_back(lower_norm); + } + ++result.iterations; + + if (lower_norm < opts.tolerance) { + extract_eigenvalues(Ak, extract_tol, + result.eigenvalues_real, + result.eigenvalues_imag); + return result; + } + } + + // Maximum iterations reached without convergence — fail loudly. + // Possible causes: eigenvalues too close in magnitude, or complex + // eigenvalue pairs that require a double shift (see Stage 2). + std::ostringstream oss; + oss << "eigenvalues_unshifted: did not converge in " + << opts.max_iterations << " iterations " + << "(final ||lower(A_k)||_F = " << lower_triangle_norm(Ak) + << ", tolerance = " << opts.tolerance << "). " + << "Try eigenvalues_shifted (Stage 2) or increase max_iterations."; + throw NonConvergenceError(oss.str()); +} + +} // namespace linalg diff --git a/tests/test_qr_iteration.cpp b/tests/test_qr_iteration.cpp new file mode 100644 index 0000000..c994c95 --- /dev/null +++ b/tests/test_qr_iteration.cpp @@ -0,0 +1,293 @@ +// Tests for qr_iteration.hpp / qr_iteration.cpp +// +// Stage 1: Unshifted QR iteration. +// +// All Stage 1 tests use symmetric matrices (only real eigenvalues) because +// the unshifted algorithm converges to upper-triangular form — not merely +// quasi-upper-triangular — only when all eigenvalues are real. A matrix +// with a complex-conjugate pair would stall: its 2×2 Schur block keeps a +// non-negligible subdiagonal entry indefinitely, so ||lower(A_k)||_F never +// falls below the tolerance. Proper handling of complex pairs requires the +// double-shift strategy introduced in Stage 2. +// +// Refs: T&B Lecture 28; GVL §7.3–7.4. + +#include "linalg_error.hpp" +#include "matrix.hpp" +#include "qr_iteration.hpp" +#include "vector.hpp" + +#include +#include + +#include +#include +#include +#include +#include +#include + +using linalg::Matrix; +using linalg::NonConvergenceError; +using linalg::QRIterationOptions; +using linalg::QRIterationResult; +using linalg::Vector; + +// --------------------------------------------------------------------------- +// Test helpers +// --------------------------------------------------------------------------- + +namespace { + +// Sort (real, imag) eigenvalue pairs by real part (ascending), then by imag. +// Returns a std::vector> — a plain container of +// pairs, not a math vector. +using EigPairs = std::vector>; + +EigPairs to_pairs(const Vector& real_v, const Vector& imag_v) { + EigPairs out; + out.reserve(real_v.size()); + for (std::size_t i = 0; i < real_v.size(); ++i) + out.emplace_back(real_v[i], imag_v[i]); + std::sort(out.begin(), out.end(), + [](const std::pair& a, + const std::pair& b) { + return a.first != b.first ? a.first < b.first + : a.second < b.second; + }); + return out; +} + +// Return true when every computed eigenvalue is within `tol` of the +// corresponding expected eigenvalue (after sorting both sets). +// `expected` is a plain std::vector of (real, imag) pairs used as test data. +bool eigs_match(const Vector& computed_real, const Vector& computed_imag, + const EigPairs& expected, double tol) { + if (computed_real.size() != expected.size()) return false; + const EigPairs computed = to_pairs(computed_real, computed_imag); + EigPairs exp_sorted = expected; + std::sort(exp_sorted.begin(), exp_sorted.end(), + [](const std::pair& a, + const std::pair& b) { + return a.first != b.first ? a.first < b.first + : a.second < b.second; + }); + for (std::size_t i = 0; i < computed.size(); ++i) { + const double dr = computed[i].first - exp_sorted[i].first; + const double di = computed[i].second - exp_sorted[i].second; + if (std::sqrt(dr * dr + di * di) > tol) return false; + } + return true; +} + +} // namespace + +// --------------------------------------------------------------------------- +// Test 1: 2×2 symmetric matrix with known eigenvalues +// --------------------------------------------------------------------------- +// +// A = | 2 1 | is symmetric positive definite. +// | 1 2 | +// +// Characteristic polynomial: (2-λ)^2 - 1 = 0 → λ = 1, 3. +// Eigenvectors: [1,-1]/√2 (λ=1) and [1,1]/√2 (λ=3). +// +// The unshifted iteration converges at rate |λ_1/λ_2| = 1/3 per step, +// so only a handful of iterations are needed. +// Ref: T&B Theorem 28.2. + +TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues", + "[qr_iteration][stage1]") { + const Matrix A{ + {2.0, 1.0}, + {1.0, 2.0} + }; + + const QRIterationResult res = linalg::eigenvalues_unshifted(A); + + REQUIRE(res.eigenvalues_real.size() == 2); + REQUIRE(res.eigenvalues_imag.size() == 2); + REQUIRE(res.iterations > 0); + + // All eigenvalues of a symmetric matrix must be real. + CHECK(std::abs(res.eigenvalues_imag[0]) < 1e-8); + CHECK(std::abs(res.eigenvalues_imag[1]) < 1e-8); + + const EigPairs expected = {{1.0, 0.0}, {3.0, 0.0}}; + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); +} + +// --------------------------------------------------------------------------- +// Test 2: 4×4 symmetric tridiagonal matrix — reference eigenvalues +// --------------------------------------------------------------------------- +// +// The n×n symmetric tridiagonal matrix with 2 on the diagonal and -1 on the +// first super- and sub-diagonals has known eigenvalues (discrete Laplacian): +// +// λ_k = 2 - 2 cos(k π / (n+1)), k = 1, …, n +// +// Ref: Golub & Van Loan §4.4.2 (discrete sine transform). +// +// For n = 4: +// λ_1 = 2 - 2 cos(π/5) ≈ 0.3820 +// λ_2 = 2 - 2 cos(2π/5) ≈ 1.3820 +// λ_3 = 2 - 2 cos(3π/5) ≈ 2.6180 +// λ_4 = 2 - 2 cos(4π/5) ≈ 3.6180 + +TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal", + "[qr_iteration][stage1]") { + const Matrix A{ + { 2.0, -1.0, 0.0, 0.0}, + {-1.0, 2.0, -1.0, 0.0}, + { 0.0, -1.0, 2.0, -1.0}, + { 0.0, 0.0, -1.0, 2.0} + }; + + const QRIterationResult res = linalg::eigenvalues_unshifted(A); + + REQUIRE(res.eigenvalues_real.size() == 4); + REQUIRE(res.eigenvalues_imag.size() == 4); + + // All eigenvalues of a symmetric matrix must be real. + for (std::size_t k = 0; k < 4; ++k) + CHECK(std::abs(res.eigenvalues_imag[k]) < 1e-8); + + // Compare against the closed-form reference. + constexpr double pi = 3.14159265358979323846; + const EigPairs expected = { + {2.0 - 2.0 * std::cos( pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(2.0 * pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0} + }; + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); +} + +// --------------------------------------------------------------------------- +// Test 3: 5×5 symmetric tridiagonal — convergence history +// --------------------------------------------------------------------------- +// +// Uses a 5×5 symmetric tridiagonal (discrete Laplacian) to guarantee all +// real eigenvalues and predictable linear convergence. The Frobenius norm +// of the strict lower triangle is printed at every step so the convergence +// rate can be observed directly. +// +// Expected behaviour: ||lower(A_k)||_F decreases geometrically each step +// (linear convergence), with ratio ≈ max_j |λ_{j+1}/λ_j|. +// Ref: T&B Theorem 28.2. +// +// 5×5 tridiagonal eigenvalues: λ_k = 2 - 2cos(kπ/6), k = 1..5. + +TEST_CASE("QR iteration (unshifted): 5x5 convergence history", + "[qr_iteration][stage1]") { + const Matrix A{ + { 2.0, -1.0, 0.0, 0.0, 0.0}, + {-1.0, 2.0, -1.0, 0.0, 0.0}, + { 0.0, -1.0, 2.0, -1.0, 0.0}, + { 0.0, 0.0, -1.0, 2.0, -1.0}, + { 0.0, 0.0, 0.0, -1.0, 2.0} + }; + + QRIterationOptions opts; + opts.track_convergence = true; + + const QRIterationResult res = linalg::eigenvalues_unshifted(A, opts); + + REQUIRE_FALSE(res.convergence_history.empty()); + REQUIRE(res.eigenvalues_real.size() == 5); + + // Print convergence history so the linear rate is visible. + std::cout << "\n=== Stage 1: Unshifted QR — 5x5 convergence history ===\n"; + std::cout << " Converged in " << res.iterations << " iteration(s)\n"; + for (std::size_t k = 0; k < res.convergence_history.size(); ++k) { + std::cout << " iter " << (k + 1) + << ": ||lower(A_k)||_F = " + << res.convergence_history[k] << "\n"; + } + std::cout << "=======================================================\n"; + + // The final recorded norm must be below the default tolerance. + CHECK(res.convergence_history.back() < opts.tolerance); +} + +// --------------------------------------------------------------------------- +// Test 4: Eigenvalue residuals below 1e-8 +// --------------------------------------------------------------------------- +// +// For several symmetric matrices with analytically known eigenvalues, verify +// that every computed eigenvalue is within 1e-8 of its expected value. +// +// Residual means the absolute error |λ_computed - λ_exact| (eigenvalue +// accuracy), not a matrix residual ||A x - λ x||, which would require +// eigenvectors unavailable in Stage 1. + +TEST_CASE("QR iteration (unshifted): residuals below 1e-8", + "[qr_iteration][stage1]") { + + SECTION("2x2: eigenvalues 1 and 3") { + const Matrix A{{2.0, 1.0}, {1.0, 2.0}}; + const QRIterationResult res = linalg::eigenvalues_unshifted(A); + const EigPairs expected = {{1.0, 0.0}, {3.0, 0.0}}; + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); + } + + SECTION("3x3 diagonal: eigenvalues 1, 4, 9") { + // Diagonal matrix — already in Schur form; converges in one step. + const Matrix D{ + {1.0, 0.0, 0.0}, + {0.0, 4.0, 0.0}, + {0.0, 0.0, 9.0} + }; + const QRIterationResult res = linalg::eigenvalues_unshifted(D); + const EigPairs expected = {{1.0, 0.0}, {4.0, 0.0}, {9.0, 0.0}}; + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); + } + + SECTION("4x4 tridiagonal: closed-form eigenvalues") { + const Matrix A{ + { 2.0, -1.0, 0.0, 0.0}, + {-1.0, 2.0, -1.0, 0.0}, + { 0.0, -1.0, 2.0, -1.0}, + { 0.0, 0.0, -1.0, 2.0} + }; + constexpr double pi = 3.14159265358979323846; + const EigPairs expected = { + {2.0 - 2.0 * std::cos( pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(2.0 * pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0} + }; + const QRIterationResult res = linalg::eigenvalues_unshifted(A); + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); + } + + SECTION("5x5 identity: all eigenvalues == 1") { + const Matrix I = Matrix::identity(5); + const QRIterationResult res = linalg::eigenvalues_unshifted(I); + REQUIRE(res.eigenvalues_real.size() == 5); + for (std::size_t k = 0; k < 5; ++k) { + CHECK(std::abs(res.eigenvalues_real[k] - 1.0) < 1e-8); + CHECK(std::abs(res.eigenvalues_imag[k]) < 1e-8); + } + } +} + +// --------------------------------------------------------------------------- +// Failure cases +// --------------------------------------------------------------------------- + +TEST_CASE("QR iteration (unshifted): non-square matrix throws", + "[qr_iteration][stage1]") { + const Matrix A(3, 4); // non-square + CHECK_THROWS_AS(linalg::eigenvalues_unshifted(A), + linalg::DimensionMismatchError); +} + +TEST_CASE("QR iteration (unshifted): max_iterations exceeded throws", + "[qr_iteration][stage1]") { + // Cap at zero iterations — any non-trivial matrix fails immediately. + const Matrix A{{2.0, 1.0}, {1.0, 2.0}}; + QRIterationOptions opts; + opts.max_iterations = 0; + CHECK_THROWS_AS(linalg::eigenvalues_unshifted(A, opts), NonConvergenceError); +} -- cgit v1.2.3