Add QR iteration and cleanup

1 files changed, 293 insertions, 0 deletions

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 <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <iostream>
+#include <utility>
+#include <vector>
+
+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<std::pair<double,double>> — a plain container of
+// pairs, not a math vector.
+using EigPairs = std::vector<std::pair<double, double>>;
+
+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<double, double>& a,
+                 const std::pair<double, double>& 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<double, double>& a,
+                 const std::pair<double, double>& 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);
+}