diff options
Diffstat (limited to 'tests/test_qr_iteration.cpp')
| -rw-r--r-- | tests/test_qr_iteration.cpp | 200 |
1 files changed, 32 insertions, 168 deletions
diff --git a/tests/test_qr_iteration.cpp b/tests/test_qr_iteration.cpp index 7126334..008f109 100644 --- a/tests/test_qr_iteration.cpp +++ b/tests/test_qr_iteration.cpp @@ -1,7 +1,4 @@ -#include "linalg_error.hpp" -#include "matrix.hpp" -#include "qr_iteration.hpp" -#include "vector.hpp" +import linalgebra; #include <catch2/catch_approx.hpp> #include <catch2/catch_test_macros.hpp> @@ -16,19 +13,14 @@ #include <utility> #include <vector> -using linalg::Matrix; -using linalg::NonConvergenceError; -using linalg::QRIterationOptions; -using linalg::QRIterationResult; -using linalg::Vector; - -// --------------------------------------------------------------------------- -// Test helpers -// --------------------------------------------------------------------------- +using linalgebra::Matrix; +using linalgebra::NonConvergenceError; +using linalgebra::QRIterationOptions; +using linalgebra::QRIterationResult; +using linalgebra::Vector; namespace { -// Sort (real, imag) eigenvalue pairs by real part (ascending), then by imag. using EigPairs = std::vector<std::pair<double, double>>; EigPairs to_pairs(const Vector& real_v, const Vector& imag_v) { @@ -66,20 +58,6 @@ bool eigs_match(const Vector& computed_real, const Vector& computed_imag, } // 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][shifted]") { const Matrix A{ @@ -87,13 +65,12 @@ TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues", {1.0, 2.0} }; - const QRIterationResult res = linalg::eigenvalues_unshifted(A); + const QRIterationResult res = linalgebra::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); @@ -101,23 +78,6 @@ TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues", 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][unshifted]") { const Matrix A{ @@ -127,16 +87,14 @@ TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal", { 0.0, 0.0, -1.0, 2.0} }; - const QRIterationResult res = linalg::eigenvalues_unshifted(A); + const QRIterationResult res = linalgebra::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}, @@ -147,21 +105,6 @@ TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal", 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][unshifted]") { const Matrix A{ @@ -175,12 +118,11 @@ TEST_CASE("QR iteration (unshifted): 5x5 convergence history", QRIterationOptions opts; opts.track_convergence = true; - const QRIterationResult res = linalg::eigenvalues_unshifted(A, opts); + const QRIterationResult res = linalgebra::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=== 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) { @@ -190,32 +132,26 @@ TEST_CASE("QR iteration (unshifted): 5x5 convergence history", } 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 -// --------------------------------------------------------------------------- - TEST_CASE("QR iteration (unshifted): residuals below 1e-8", "[qr_iteration][unshifted]") { 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 QRIterationResult res = linalgebra::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 QRIterationResult res = linalgebra::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)); } @@ -234,13 +170,13 @@ TEST_CASE("QR iteration (unshifted): residuals below 1e-8", {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); + const QRIterationResult res = linalgebra::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); + const QRIterationResult res = linalgebra::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); @@ -249,35 +185,21 @@ TEST_CASE("QR iteration (unshifted): residuals below 1e-8", } } -// --------------------------------------------------------------------------- -// Failure cases -// --------------------------------------------------------------------------- - TEST_CASE("QR iteration (unshifted): non-square matrix throws", "[qr_iteration][unshifted]") { - const Matrix A(3, 4); // non-square - CHECK_THROWS_AS(linalg::eigenvalues_unshifted(A), - linalg::DimensionMismatchError); + const Matrix A(3, 4); + CHECK_THROWS_AS(linalgebra::eigenvalues_unshifted(A), + linalgebra::DimensionMismatchError); } TEST_CASE("QR iteration (unshifted): max_iterations exceeded throws", "[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); + CHECK_THROWS_AS(linalgebra::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) { @@ -296,14 +218,6 @@ Matrix random_symmetric(std::size_t n, unsigned seed = 42) { } // 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{ @@ -318,8 +232,8 @@ TEST_CASE("QR iteration (shifted): fewer iterations than unshifted", QRIterationOptions opts; opts.track_convergence = true; - const QRIterationResult unshifted = linalg::eigenvalues_unshifted(A, opts); - const QRIterationResult shifted = linalg::eigenvalues_shifted(A, opts); + const QRIterationResult unshifted = linalgebra::eigenvalues_unshifted(A, opts); + const QRIterationResult shifted = linalgebra::eigenvalues_shifted(A, opts); std::cout << "\n=== Shifted vs Unshifted ===\n"; std::cout << " Unshifted iterations: " << unshifted.iterations << "\n"; @@ -333,25 +247,15 @@ TEST_CASE("QR iteration (shifted): fewer iterations than unshifted", 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; + opts.max_iterations = 2000; - const QRIterationResult unshifted = linalg::eigenvalues_unshifted(A, opts); - const QRIterationResult shifted = linalg::eigenvalues_shifted(A, opts); + const QRIterationResult unshifted = linalgebra::eigenvalues_unshifted(A, opts); + const QRIterationResult shifted = linalgebra::eigenvalues_shifted(A, opts); std::cout << "\n=== Hard matrix ===\n"; std::cout << " Unshifted iterations: " << unshifted.iterations << "\n"; @@ -361,10 +265,6 @@ TEST_CASE("QR iteration (shifted): converges <20 iters where unshifted needs >10 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") { @@ -381,29 +281,18 @@ TEST_CASE("QR iteration (shifted): residuals below 1e-8", {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); + const QRIterationResult res = linalgebra::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); + const QRIterationResult res = linalgebra::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) { @@ -421,7 +310,6 @@ double frobenius_norm(const Matrix& A) { 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) @@ -432,7 +320,6 @@ double diff_norm(const Matrix& A, const Matrix& B) { 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; @@ -451,29 +338,21 @@ double orthogonality_error(const Matrix& Q) { 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); + const linalgebra::HessenbergResult hr = linalgebra::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); + const Matrix QtHQ = hr.Q * hr.H * linalgebra::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); + const QRIterationResult ref = linalgebra::eigenvalues_shifted(A); + const QRIterationResult hess = linalgebra::eigenvalues_hessenberg(A); REQUIRE(hess.eigenvalues_real.size() == 8); CHECK(eigs_match(hess.eigenvalues_real, hess.eigenvalues_imag, @@ -481,10 +360,6 @@ TEST_CASE("Hessenberg QR: eigenvalues match shifted QR to 1e-6", 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{ @@ -500,19 +375,10 @@ TEST_CASE("Hessenberg QR: residuals below 1e-8 on known matrix", {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); + const QRIterationResult res = linalgebra::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; @@ -531,11 +397,11 @@ TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n", const Matrix A = random_symmetric(n, 13u); const auto t0s = Clock::now(); - { const auto tmp = linalg::eigenvalues_shifted(A); (void)tmp; } + { const auto tmp = linalgebra::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 QRIterationResult hess = linalgebra::eigenvalues_hessenberg(A); const double t_hess = Seconds(Clock::now() - t0h).count(); const double speedup = t_shifted / t_hess; @@ -547,11 +413,9 @@ TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n", << 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); + const QRIterationResult ref = linalgebra::eigenvalues_shifted(A); CHECK(eigs_match(hess.eigenvalues_real, hess.eigenvalues_imag, to_pairs(ref.eigenvalues_real, ref.eigenvalues_imag), 1e-6)); |