From 750f276a1403c5defd58b06f13c703b5e3d59245 Mon Sep 17 00:00:00 2001 From: y-jan137 Date: Sat, 16 May 2026 12:10:05 +0300 Subject: Finish TODOs --- tests/test_expm.cpp | 144 ++++++++++++++++++++++++++++++ tests/test_iterative.cpp | 211 ++++++++++++++++++++++++++++++++++++++++++++ tests/test_precond.cpp | 142 +++++++++++++++++++++++++++++ tests/test_qr.cpp | 2 - tests/test_qr_iteration.cpp | 2 - tests/test_svd.cpp | 140 +++++++++++++++++++++++++++++ tests/test_sym_eigen.cpp | 168 +++++++++++++++++++++++++++++++++++ 7 files changed, 805 insertions(+), 4 deletions(-) create mode 100644 tests/test_expm.cpp create mode 100644 tests/test_iterative.cpp create mode 100644 tests/test_precond.cpp create mode 100644 tests/test_svd.cpp create mode 100644 tests/test_sym_eigen.cpp (limited to 'tests') diff --git a/tests/test_expm.cpp b/tests/test_expm.cpp new file mode 100644 index 0000000..a0d3ea0 --- /dev/null +++ b/tests/test_expm.cpp @@ -0,0 +1,144 @@ +import linalgebra; + +#include +#include + +#include +#include +#include +#include + +using linalgebra::DimensionMismatchError; +using linalgebra::Matrix; + +namespace { + +double frobenius_diff(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); +} + +// Taylor series +Matrix expm_taylor5(const Matrix& A) { + const std::size_t n = A.rows(); + Matrix result = Matrix::identity(n); + Matrix power = Matrix::identity(n); + double fact = 1.0; + for (int k = 1; k <= 5; ++k) { + power = power * A; + fact *= static_cast(k); + const double inv_fact = 1.0 / fact; + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + result(i, j) += inv_fact * power(i, j); + } + return result; +} + +} // namespace + +TEST_CASE("expm: zero matrix gives identity", "[expm]") { + Matrix Z = Matrix::zeros(3, 3); + auto E = linalgebra::expm(Z); + auto I = Matrix::identity(3); + REQUIRE(frobenius_diff(E, I) < 1e-12); +} + +TEST_CASE("expm: scalar multiple of identity", "[expm]") { + // expm(s*I) = exp(s)*I + const double s = 2.0; + Matrix A = Matrix::zeros(3, 3); + A(0, 0) = s; A(1, 1) = s; A(2, 2) = s; + auto E = linalgebra::expm(A); + + const double expected = std::exp(s); + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + REQUIRE(E(i, j) == Catch::Approx(i == j ? expected : 0.0).margin(1e-10)); +} + +TEST_CASE("expm: 2x2 nilpotent", "[expm]") { + // A = [[0,1],[0,0]], expm(A) = [[1,1],[0,1]] exactly. + Matrix A{{0.0, 1.0}, {0.0, 0.0}}; + auto E = linalgebra::expm(A); + REQUIRE(E(0, 0) == Catch::Approx(1.0).margin(1e-12)); + REQUIRE(E(0, 1) == Catch::Approx(1.0).margin(1e-12)); + REQUIRE(E(1, 0) == Catch::Approx(0.0).margin(1e-12)); + REQUIRE(E(1, 1) == Catch::Approx(1.0).margin(1e-12)); +} + +TEST_CASE("expm: 2x2 rotation generator", "[expm]") { + // A = [[0,-t],[t,0]], expm(A) = [[cos(t), -sin(t)],[sin(t), cos(t)]]. + const double t = std::numbers::pi / 4.0; + Matrix A{{0.0, -t}, {t, 0.0}}; + auto E = linalgebra::expm(A); + + REQUIRE(E(0, 0) == Catch::Approx(std::cos(t)).epsilon(1e-10)); + REQUIRE(E(0, 1) == Catch::Approx(-std::sin(t)).epsilon(1e-10)); + REQUIRE(E(1, 0) == Catch::Approx(std::sin(t)).epsilon(1e-10)); + REQUIRE(E(1, 1) == Catch::Approx(std::cos(t)).epsilon(1e-10)); +} + +TEST_CASE("expm: diagonal matrix", "[expm]") { + // A = diag(1, 2), expm(A) = diag(e, e^2). + Matrix A{{1.0, 0.0}, {0.0, 2.0}}; + auto E = linalgebra::expm(A); + REQUIRE(E(0, 0) == Catch::Approx(std::exp(1.0)).epsilon(1e-10)); + REQUIRE(E(1, 1) == Catch::Approx(std::exp(2.0)).epsilon(1e-10)); + REQUIRE(E(0, 1) == Catch::Approx(0.0).margin(1e-12)); + REQUIRE(E(1, 0) == Catch::Approx(0.0).margin(1e-12)); +} + +TEST_CASE("expm: comparison with Taylor series (small-norm A)", "[expm]") { + // For small ||A||, expm(A) ≈ Taylor series to order 5. + std::mt19937 rng(42); + std::uniform_real_distribution dist(-0.01, 0.01); + Matrix A(4, 4); + for (std::size_t i = 0; i < 4; ++i) + for (std::size_t j = 0; j < 4; ++j) + A(i, j) = dist(rng); + + auto E = linalgebra::expm(A); + auto E_taylor = expm_taylor5(A); + REQUIRE(frobenius_diff(E, E_taylor) < 1e-10); +} + +TEST_CASE("expm: large norm requires scaling (10*I)", "[expm]") { + // expm(10*I) = exp(10)*I. + Matrix A = Matrix::zeros(3, 3); + A(0, 0) = 10.0; A(1, 1) = 10.0; A(2, 2) = 10.0; + auto E = linalgebra::expm(A); + + const double e10 = std::exp(10.0); // ≈ 22026.47 + for (std::size_t i = 0; i < 3; ++i) + REQUIRE(E(i, i) == Catch::Approx(e10).epsilon(1e-8)); +} + +TEST_CASE("expm: non-square throws", "[expm]") { + Matrix A(2, 3); + REQUIRE_THROWS_AS(linalgebra::expm(A), DimensionMismatchError); +} + +TEST_CASE("expm: symmetric A gives symmetric expm(A)", "[expm]") { + std::mt19937 rng(7); + std::uniform_real_distribution dist(-1.0, 1.0); + Matrix B(4, 4); + for (std::size_t i = 0; i < 4; ++i) + for (std::size_t j = 0; j < 4; ++j) + B(i, j) = dist(rng); + // Make symmetric: A = B + B^T. + Matrix A(4, 4, 0.0); + for (std::size_t i = 0; i < 4; ++i) + for (std::size_t j = 0; j < 4; ++j) + A(i, j) = B(i, j) + B(j, i); + + auto E = linalgebra::expm(A); + for (std::size_t i = 0; i < 4; ++i) + for (std::size_t j = 0; j < 4; ++j) + REQUIRE(E(i, j) == Catch::Approx(E(j, i)).margin(1e-9)); +} diff --git a/tests/test_iterative.cpp b/tests/test_iterative.cpp new file mode 100644 index 0000000..64ee85d --- /dev/null +++ b/tests/test_iterative.cpp @@ -0,0 +1,211 @@ +import linalgebra; + +#include +#include + +#include +#include +#include + +using linalgebra::DimensionMismatchError; +using linalgebra::LinAlgError; +using linalgebra::Matrix; +using linalgebra::Vector; + +namespace { + +double frobenius_diff(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); +} + +Matrix make_dd(std::size_t n, std::mt19937& rng) { + std::uniform_real_distribution dist(0.0, 1.0); + Matrix A(n, n); + for (std::size_t i = 0; i < n; ++i) { + double row_sum = 0.0; + for (std::size_t j = 0; j < n; ++j) { + A(i, j) = dist(rng); + if (i != j) row_sum += std::abs(A(i, j)); + } + A(i, i) = row_sum + 1.0; // strictly diagonally dominant + } + return A; +} + +Matrix make_spd(std::size_t n, std::mt19937& rng) { + std::uniform_real_distribution dist(-1.0, 1.0); + Matrix B(n, n); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + B(i, j) = dist(rng); + Matrix A = linalgebra::transpose(B) * B; + for (std::size_t i = 0; i < n; ++i) A(i, i) += static_cast(n); + return A; +} + +} // namespace + +// arnoldi + +TEST_CASE("arnoldi: orthonormality and AQ=QH relation", "[iterative][arnoldi]") { + Matrix A{{2.0, 1.0, 0.0}, + {1.0, 3.0, 1.0}, + {0.0, 1.0, 4.0}}; + Vector b{1.0, 0.0, 0.0}; + const int k = 2; + auto res = linalgebra::arnoldi(A, b, k); + + const std::size_t n = 3; + const auto steps = static_cast(res.steps_taken); + + // Q columns must be orthonormal. + for (std::size_t i = 0; i <= steps; ++i) { + for (std::size_t j = 0; j <= steps; ++j) { + double dot = 0.0; + for (std::size_t r = 0; r < n; ++r) dot += res.Q(r, i) * res.Q(r, j); + const double expected = (i == j) ? 1.0 : 0.0; + REQUIRE(dot == Catch::Approx(expected).margin(1e-10)); + } + } + + for (std::size_t j = 0; j < steps; ++j) { + Vector qj(n); + for (std::size_t i = 0; i < n; ++i) qj[i] = res.Q(i, j); + const Vector Aqj = A * qj; + + for (std::size_t i = 0; i <= steps; ++i) { + double qh = 0.0; + for (std::size_t r = 0; r < n; ++r) qh += res.Q(r, i) * res.H(i, j); + } + + + // Direct check: ||A*qj - Q*H[:,j]|| + Vector Hcol(steps + 1); + for (std::size_t i = 0; i <= steps; ++i) Hcol[i] = res.H(i, j); + double resid = 0.0; + for (std::size_t row = 0; row < n; ++row) { + double qh_row = 0.0; + for (std::size_t i = 0; i <= steps; ++i) qh_row += res.Q(row, i) * Hcol[i]; + const double d = Aqj[row] - qh_row; + resid += d * d; + } + REQUIRE(std::sqrt(resid) < 1e-10); + } +} + +TEST_CASE("arnoldi: breakdown on scaled identity", "[iterative][arnoldi]") { + // A = 2*I → Krylov space is one-dimensional. + Matrix A(3, 3, 0.0); + A(0, 0) = 2.0; A(1, 1) = 2.0; A(2, 2) = 2.0; + Vector b{1.0, 0.0, 0.0}; + auto res = linalgebra::arnoldi(A, b, 3); + REQUIRE(res.breakdown == true); + REQUIRE(res.steps_taken <= 3); +} + +TEST_CASE("arnoldi: dimension mismatch throws", "[iterative][arnoldi]") { + Matrix A(3, 3, 0.0); + Vector b(4, 0.0); + REQUIRE_THROWS_AS(linalgebra::arnoldi(A, b, 2), DimensionMismatchError); +} + +// solve_cg + +TEST_CASE("solve_cg: 2x2 SPD", "[iterative][cg]") { + Matrix A{{4.0, 1.0}, {1.0, 3.0}}; + Vector b{1.0, 2.0}; + auto res = linalgebra::solve_cg(A, b); + REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-10); +} + +TEST_CASE("solve_cg: identity system converges in 1 iteration", "[iterative][cg]") { + auto I = Matrix::identity(5); + Vector b{1.0, 2.0, 3.0, 4.0, 5.0}; + auto res = linalgebra::solve_cg(I, b); + for (std::size_t i = 0; i < 5; ++i) + REQUIRE(res.x[i] == Catch::Approx(b[i]).margin(1e-10)); +} + +TEST_CASE("solve_cg: 5x5 random SPD", "[iterative][cg]") { + std::mt19937 rng(99); + Matrix A = make_spd(5, rng); + std::uniform_real_distribution dist(-2.0, 2.0); + Vector b(5); + for (std::size_t i = 0; i < 5; ++i) b[i] = dist(rng); + + auto res = linalgebra::solve_cg(A, b); + REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-8); + + // Cross-check vs LU. + auto lu = linalgebra::lu_factor(A); + auto x_lu = linalgebra::lu_solve(lu, b); + REQUIRE(linalgebra::norm2(res.x - x_lu) < 1e-8); +} + +TEST_CASE("solve_cg: non-symmetric throws", "[iterative][cg]") { + Matrix A{{1.0, 2.0}, {0.0, 1.0}}; + Vector b{1.0, 1.0}; + REQUIRE_THROWS_AS(linalgebra::solve_cg(A, b), LinAlgError); +} + +// solve_gmres + +TEST_CASE("solve_gmres: 2x2 non-symmetric", "[iterative][gmres]") { + Matrix A{{2.0, 1.0}, {1.0, 3.0}}; + Vector b{5.0, 7.0}; + auto res = linalgebra::solve_gmres(A, b); + REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-9); +} + +TEST_CASE("solve_gmres: 5x5 diagonally dominant", "[iterative][gmres]") { + std::mt19937 rng(11); + Matrix A = make_dd(5, rng); + std::uniform_real_distribution dist(-3.0, 3.0); + Vector b(5); + for (std::size_t i = 0; i < 5; ++i) b[i] = dist(rng); + + auto res = linalgebra::solve_gmres(A, b); + REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-8); + + auto lu = linalgebra::lu_factor(A); + auto x_lu = linalgebra::lu_solve(lu, b); + REQUIRE(linalgebra::norm2(res.x - x_lu) < 1e-8); +} + +TEST_CASE("solve_gmres: zero rhs gives zero solution", "[iterative][gmres]") { + Matrix A{{2.0, 1.0}, {1.0, 3.0}}; + Vector b(2, 0.0); + auto res = linalgebra::solve_gmres(A, b); + REQUIRE(linalgebra::norm2(res.x) < 1e-12); +} + + +// solve_bicgstab + +TEST_CASE("solve_bicgstab: 2x2 non-symmetric", "[iterative][bicgstab]") { + Matrix A{{2.0, 1.0}, {1.0, 3.0}}; + Vector b{5.0, 7.0}; + auto res = linalgebra::solve_bicgstab(A, b); + REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-9); +} + +TEST_CASE("solve_bicgstab: 5x5 diagonally dominant", "[iterative][bicgstab]") { + std::mt19937 rng(22); + Matrix A = make_dd(5, rng); + std::uniform_real_distribution dist(-3.0, 3.0); + Vector b(5); + for (std::size_t i = 0; i < 5; ++i) b[i] = dist(rng); + + auto res = linalgebra::solve_bicgstab(A, b); + REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-8); + + auto lu = linalgebra::lu_factor(A); + auto x_lu = linalgebra::lu_solve(lu, b); + REQUIRE(linalgebra::norm2(res.x - x_lu) < 1e-8); +} diff --git a/tests/test_precond.cpp b/tests/test_precond.cpp new file mode 100644 index 0000000..160badf --- /dev/null +++ b/tests/test_precond.cpp @@ -0,0 +1,142 @@ +import linalgebra; + +#include +#include +#include + +using linalgebra::DimensionMismatchError; +using linalgebra::Matrix; +using linalgebra::SingularMatrixError; +using linalgebra::Vector; + +// condition_number_1norm + +TEST_CASE("condition_number_1norm: identity", "[precond][condition]") { + for (std::size_t n : {1, 2, 5}) { + auto I = Matrix::identity(n); + REQUIRE(linalgebra::condition_number_1norm(I) == Catch::Approx(1.0).epsilon(1e-10)); + } +} + +TEST_CASE("condition_number_1norm: diagonal matrix", "[precond][condition]") { + // diag(1, 10, 100): ||A||_1 = 100, ||A^{-1}||_1 = 1, so cond = 100. + Matrix A(3, 3, 0.0); + A(0, 0) = 1.0; A(1, 1) = 10.0; A(2, 2) = 100.0; + const double c = linalgebra::condition_number_1norm(A); + // Power-iteration estimator gives a lower bound; for simple diagonal it should be exact. + REQUIRE(c == Catch::Approx(100.0).epsilon(1e-8)); +} + +TEST_CASE("condition_number_1norm: ill-conditioned Hilbert 4x4", "[precond][condition]") { + // Hilbert matrix H(i,j) = 1/(i+j+1). + Matrix H(4, 4); + for (std::size_t i = 0; i < 4; ++i) + for (std::size_t j = 0; j < 4; ++j) + H(i, j) = 1.0 / static_cast(i + j + 1); + const double c = linalgebra::condition_number_1norm(H); + REQUIRE(c > 1e3); // Hilbert matrices are notoriously ill-conditioned +} + +TEST_CASE("condition_number_1norm: non-square throws", "[precond][condition]") { + Matrix A(2, 3); + REQUIRE_THROWS_AS(linalgebra::condition_number_1norm(A), DimensionMismatchError); +} + +// precond_jacobi + +TEST_CASE("precond_jacobi: diagonal matrix", "[precond][jacobi]") { + Matrix A(3, 3, 0.0); + A(0, 0) = 2.0; A(1, 1) = 4.0; A(2, 2) = 8.0; + auto P = linalgebra::precond_jacobi(A); + + REQUIRE(P.inv_diag[0] == Catch::Approx(0.5)); + REQUIRE(P.inv_diag[1] == Catch::Approx(0.25)); + REQUIRE(P.inv_diag[2] == Catch::Approx(0.125)); + + Vector x{1.0, 1.0, 1.0}; + auto y = linalgebra::apply(P, x); + REQUIRE(y[0] == Catch::Approx(0.5)); + REQUIRE(y[1] == Catch::Approx(0.25)); + REQUIRE(y[2] == Catch::Approx(0.125)); +} + +TEST_CASE("precond_jacobi: zero diagonal throws", "[precond][jacobi]") { + Matrix A{{1.0, 0.0}, {0.0, 0.0}}; + REQUIRE_THROWS_AS(linalgebra::precond_jacobi(A), SingularMatrixError); +} + +TEST_CASE("precond_jacobi: non-square throws", "[precond][jacobi]") { + Matrix A(2, 3); + REQUIRE_THROWS_AS(linalgebra::precond_jacobi(A), DimensionMismatchError); +} + +// precond_ilu0 + +TEST_CASE("precond_ilu0: identity", "[precond][ilu0]") { + auto I = Matrix::identity(3); + auto P = linalgebra::precond_ilu0(I); + Vector b{3.0, 1.0, 4.0}; + auto y = linalgebra::apply(P, b); + for (std::size_t i = 0; i < 3; ++i) + REQUIRE(y[i] == Catch::Approx(b[i]).margin(1e-12)); +} + +TEST_CASE("precond_ilu0: 3x3 SPD", "[precond][ilu0]") { + // Dense ILU0 = exact LU without pivoting, so for any nonsingular A, + // apply(P, A*x) should recover x. + Matrix A{{4.0, 2.0, 0.0}, + {2.0, 3.0, 1.0}, + {0.0, 1.0, 2.0}}; + auto P = linalgebra::precond_ilu0(A); + + Vector x_true{1.0, 2.0, 3.0}; + Vector rhs = A * x_true; + auto x_rec = linalgebra::apply(P, rhs); + + for (std::size_t i = 0; i < 3; ++i) + REQUIRE(x_rec[i] == Catch::Approx(x_true[i]).margin(1e-10)); +} + +TEST_CASE("precond_ilu0: near-singular pivot throws", "[precond][ilu0]") { + Matrix A{{0.0, 1.0}, {1.0, 1.0}}; + REQUIRE_THROWS_AS(linalgebra::precond_ilu0(A), SingularMatrixError); +} + +// lstsq + +TEST_CASE("lstsq: square full-rank system", "[precond][lstsq]") { + Matrix A{{2.0, 1.0}, {1.0, 3.0}}; + Vector b{5.0, 7.0}; + auto res = linalgebra::lstsq(A, b); + REQUIRE(res.rank == 2); + REQUIRE(linalgebra::norm2(A * res.x - b) < 1e-10); +} + +TEST_CASE("lstsq: overdetermined full-rank", "[precond][lstsq]") { + // A is 3x2, consistent overdetermined system. + Matrix A{{1.0, 1.0}, {2.0, 1.0}, {3.0, 1.0}}; + Vector b{6.0, 5.0, 7.0}; + auto res = linalgebra::lstsq(A, b); + REQUIRE(res.rank == 2); + // Residual should be the minimum achievable (verify normal equations: A^T A x = A^T b). + Vector AtAx = linalgebra::transpose(A) * (A * res.x); + Vector Atb = linalgebra::transpose(A) * b; + for (std::size_t i = 0; i < 2; ++i) + REQUIRE(AtAx[i] == Catch::Approx(Atb[i]).margin(1e-8)); +} + +TEST_CASE("lstsq: rank-deficient", "[precond][lstsq]") { + // col2 = 2*col1 → rank 1. + Matrix A{{1.0, 2.0}, {2.0, 4.0}, {3.0, 6.0}}; + Vector b{1.0, 2.0, 3.0}; + auto res = linalgebra::lstsq(A, b); + REQUIRE(res.rank < 2); + // The residual should be the minimum achievable. + REQUIRE(res.residual_norm < 1e-8); +} + +TEST_CASE("lstsq: non-tall matrix throws", "[precond][lstsq]") { + Matrix A(2, 3); + Vector b(2, 0.0); + REQUIRE_THROWS_AS(linalgebra::lstsq(A, b), DimensionMismatchError); +} diff --git a/tests/test_qr.cpp b/tests/test_qr.cpp index d78fbb1..cbde901 100644 --- a/tests/test_qr.cpp +++ b/tests/test_qr.cpp @@ -224,9 +224,7 @@ TEST_CASE("QR: linearly dependent columns throw from GS methods", "[qr]") { CHECK_NOTHROW(linalgebra::qr_householder(A)); } -// --------------------------------------------------------------------------- // qr_colpiv tests -// --------------------------------------------------------------------------- TEST_CASE("QR ColPiv: full rank reconstruction", "[qr][colpiv]") { const Matrix A = random_matrix(6, 4, 100u); diff --git a/tests/test_qr_iteration.cpp b/tests/test_qr_iteration.cpp index 64a5dd6..95bb044 100644 --- a/tests/test_qr_iteration.cpp +++ b/tests/test_qr_iteration.cpp @@ -423,9 +423,7 @@ TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n", std::cout << "=============================================\n"; } -// --------------------------------------------------------------------------- // Francis double-shift QR tests -// --------------------------------------------------------------------------- TEST_CASE("Francis QR: 2x2 real eigenvalues", "[qr_iteration][francis]") { Matrix A{{3.0, 1.0}, {0.0, 2.0}}; diff --git a/tests/test_svd.cpp b/tests/test_svd.cpp new file mode 100644 index 0000000..d2121d6 --- /dev/null +++ b/tests/test_svd.cpp @@ -0,0 +1,140 @@ +import linalgebra; + +#include +#include + +#include +#include +#include + +using linalgebra::DimensionMismatchError; +using linalgebra::Matrix; +using linalgebra::SVDResult; + +namespace { + +double frobenius_diff(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); +} + +// Reconstruct A = U * diag(sigma) * Vt and measure error +double reconstruction_error(const Matrix& A, const SVDResult& svd) { + const std::size_t m = A.rows(); + const std::size_t n = A.cols(); + const std::size_t p = svd.sigma.size(); + + // Build U_thin (m x p) — first p columns of U. + Matrix U_thin(m, p); + for (std::size_t i = 0; i < m; ++i) + for (std::size_t j = 0; j < p; ++j) + U_thin(i, j) = svd.U(i, j); + + // Build Sigma_diag (p x p). + Matrix S(p, p, 0.0); + for (std::size_t i = 0; i < p; ++i) S(i, i) = svd.sigma[i]; + + // Build Vt_thin (p x n) — first p rows of Vt. + Matrix Vt_thin(p, n); + for (std::size_t i = 0; i < p; ++i) + for (std::size_t j = 0; j < n; ++j) + Vt_thin(i, j) = svd.Vt(i, j); + + const Matrix recon = U_thin * S * Vt_thin; + return frobenius_diff(A, recon); +} + +// Check orthogonality of the first k columns of M +double col_ortho_error(const Matrix& M, std::size_t k) { + const std::size_t m = M.rows(); + double s = 0.0; + for (std::size_t i = 0; i < k; ++i) { + for (std::size_t j = 0; j < k; ++j) { + double dot = 0.0; + for (std::size_t r = 0; r < m; ++r) dot += M(r, i) * M(r, j); + const double expected = (i == j) ? 1.0 : 0.0; + const double d = dot - expected; + s += d * d; + } + } + return std::sqrt(s); +} + +} // namespace + +TEST_CASE("svd: 2x2 diagonal", "[svd]") { + Matrix A{{3.0, 0.0}, {0.0, -2.0}}; + auto res = linalgebra::svd(A); + + REQUIRE(res.sigma[0] == Catch::Approx(3.0).epsilon(1e-10)); + REQUIRE(res.sigma[1] == Catch::Approx(2.0).epsilon(1e-10)); + REQUIRE(reconstruction_error(A, res) < 1e-10); +} + +TEST_CASE("svd: 3x2 tall matrix", "[svd]") { + Matrix A{{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}}; + auto res = linalgebra::svd(A); + + REQUIRE(reconstruction_error(A, res) < 1e-10); + + REQUIRE(res.sigma[0] >= res.sigma[1]); + REQUIRE(res.sigma[1] >= 0.0); + + REQUIRE(col_ortho_error(res.U, 2) < 1e-10); + + REQUIRE(col_ortho_error(linalgebra::transpose(res.Vt), 2) < 1e-10); +} + +TEST_CASE("svd: rank-1 matrix", "[svd]") { + // A = u * v^T for u = [1,2,3], v = [2,1]. + Matrix A{{2.0, 1.0}, {4.0, 2.0}, {6.0, 3.0}}; + auto res = linalgebra::svd(A); + + REQUIRE(reconstruction_error(A, res) < 1e-10); + // Second singular value should be near zero. + REQUIRE(res.sigma[0] > 1e-10); + REQUIRE(res.sigma[1] < 1e-8); +} + +TEST_CASE("svd: identity 4x4", "[svd]") { + auto I = Matrix::identity(4); + auto res = linalgebra::svd(I); + + for (std::size_t i = 0; i < 4; ++i) + REQUIRE(res.sigma[i] == Catch::Approx(1.0).epsilon(1e-10)); + REQUIRE(reconstruction_error(I, res) < 1e-10); +} + +TEST_CASE("svd: 4x4 random", "[svd]") { + std::mt19937 rng(42); + std::uniform_real_distribution dist(-2.0, 2.0); + Matrix A(4, 4); + for (std::size_t i = 0; i < 4; ++i) + for (std::size_t j = 0; j < 4; ++j) + A(i, j) = dist(rng); + + auto res = linalgebra::svd(A); + REQUIRE(reconstruction_error(A, res) < 1e-9); + REQUIRE(col_ortho_error(res.U, 4) < 1e-10); + REQUIRE(col_ortho_error(linalgebra::transpose(res.Vt), 4) < 1e-10); +} + +TEST_CASE("svd: Hilbert 4x4", "[svd]") { + Matrix H(4, 4); + for (std::size_t i = 0; i < 4; ++i) + for (std::size_t j = 0; j < 4; ++j) + H(i, j) = 1.0 / static_cast(i + j + 1); + auto res = linalgebra::svd(H); + REQUIRE(reconstruction_error(H, res) < 1e-10); + for (std::size_t i = 0; i < 4; ++i) REQUIRE(res.sigma[i] >= 0.0); +} + +TEST_CASE("svd: wide matrix throws", "[svd]") { + Matrix A(2, 4); + REQUIRE_THROWS_AS(linalgebra::svd(A), DimensionMismatchError); +} diff --git a/tests/test_sym_eigen.cpp b/tests/test_sym_eigen.cpp new file mode 100644 index 0000000..3906785 --- /dev/null +++ b/tests/test_sym_eigen.cpp @@ -0,0 +1,168 @@ +import linalgebra; + +#include +#include + +#include +#include +#include + +using linalgebra::DimensionMismatchError; +using linalgebra::LinAlgError; +using linalgebra::Matrix; +using linalgebra::Vector; + +namespace { + +double frobenius_diff(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); +} + +// Check orthogonality ||Q^T Q - I||_F < tol. +double ortho_error(const Matrix& Q) { + const std::size_t n = Q.rows(); + const std::size_t m = Q.cols(); + const Matrix Qt = linalgebra::transpose(Q); + const Matrix QtQ = Qt * Q; + const Matrix I = Matrix::identity(m); + return frobenius_diff(QtQ, I); +} + +// Build a symmetric matrix from B^T B + shift * I. +Matrix make_sym(std::size_t n, double shift, std::mt19937& rng) { + std::uniform_real_distribution dist(-1.0, 1.0); + Matrix B(n, n); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + B(i, j) = dist(rng); + Matrix A = linalgebra::transpose(B) * B; + for (std::size_t i = 0; i < n; ++i) A(i, i) += shift; + return A; +} + +} // namespace + +// tridiagonalize + +TEST_CASE("tridiagonalize: 2x2 symmetric", "[sym_eigen][tridiagonalize]") { + Matrix A{{4.0, 2.0}, {2.0, 3.0}}; + auto res = linalgebra::tridiagonalize(A); + + // Reconstruction: Q T Q^T == A. + const Matrix QTQT = res.Q * res.T * linalgebra::transpose(res.Q); + REQUIRE(frobenius_diff(A, QTQT) < 1e-12); + + // Q must be orthogonal. + REQUIRE(ortho_error(res.Q) < 1e-12); + + // T must be tridiagonal: T(i,j) == 0 for |i-j| > 1. + const std::size_t n = res.T.rows(); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + if (i > j + 1 || j > i + 1) + REQUIRE(std::abs(res.T(i, j)) < 1e-12); +} + +TEST_CASE("tridiagonalize: 4x4 symmetric", "[sym_eigen][tridiagonalize]") { + Matrix A{{6.0, 2.0, 1.0, 0.0}, + {2.0, 5.0, 3.0, 1.0}, + {1.0, 3.0, 4.0, 2.0}, + {0.0, 1.0, 2.0, 3.0}}; + auto res = linalgebra::tridiagonalize(A); + + const Matrix QTQT = res.Q * res.T * linalgebra::transpose(res.Q); + REQUIRE(frobenius_diff(A, QTQT) < 1e-10); + REQUIRE(ortho_error(res.Q) < 1e-12); + + const std::size_t n = res.T.rows(); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + if (i > j + 1 || j > i + 1) + REQUIRE(std::abs(res.T(i, j)) < 1e-10); +} + +TEST_CASE("tridiagonalize: identity", "[sym_eigen][tridiagonalize]") { + auto I = Matrix::identity(5); + auto res = linalgebra::tridiagonalize(I); + REQUIRE(frobenius_diff(I, res.Q * res.T * linalgebra::transpose(res.Q)) < 1e-12); +} + +TEST_CASE("tridiagonalize: random symmetric", "[sym_eigen][tridiagonalize]") { + std::mt19937 rng(77); + Matrix A = make_sym(8, 5.0, rng); + auto res = linalgebra::tridiagonalize(A); + REQUIRE(frobenius_diff(A, res.Q * res.T * linalgebra::transpose(res.Q)) < 1e-9); + REQUIRE(ortho_error(res.Q) < 1e-11); +} + +TEST_CASE("tridiagonalize: non-symmetric throws", "[sym_eigen][tridiagonalize]") { + Matrix A{{1.0, 2.0}, {3.0, 4.0}}; + REQUIRE_THROWS_AS(linalgebra::tridiagonalize(A), LinAlgError); +} + +TEST_CASE("tridiagonalize: non-square throws", "[sym_eigen][tridiagonalize]") { + Matrix A(2, 3); + REQUIRE_THROWS_AS(linalgebra::tridiagonalize(A), DimensionMismatchError); +} + +// eigenvectors_inverse_iteration + +TEST_CASE("eigenvectors_inverse_iteration: 2x2 diagonal", "[sym_eigen][inverse_iter]") { + Matrix A(2, 2, 0.0); + A(0, 0) = 3.0; A(1, 1) = 7.0; + Vector lambdas{3.0, 7.0}; + + auto res = linalgebra::eigenvectors_inverse_iteration(A, lambdas); + + // Each column should satisfy A*v ≈ lambda*v. + for (std::size_t col = 0; col < 2; ++col) { + Vector v(2); + v[0] = res.eigenvectors(0, col); + v[1] = res.eigenvectors(1, col); + const Vector Av = A * v; + const double lam = lambdas[col]; + double resid = 0.0; + for (std::size_t i = 0; i < 2; ++i) { + const double d = Av[i] - lam * v[i]; + resid += d * d; + } + REQUIRE(std::sqrt(resid) < 1e-8); + } +} + +TEST_CASE("eigenvectors_inverse_iteration: 3x3 symmetric", "[sym_eigen][inverse_iter]") { + // Known symmetric matrix: compute eigenvalues with francis, then eigenvectors. + Matrix A{{6.0, 2.0, 1.0}, + {2.0, 3.0, 1.0}, + {1.0, 1.0, 1.0}}; + + auto eig = linalgebra::eigenvalues_francis(A); + // Use real eigenvalues only. + auto res = linalgebra::eigenvectors_inverse_iteration(A, eig.eigenvalues_real); + + for (std::size_t col = 0; col < 3; ++col) { + Vector v(3); + for (std::size_t i = 0; i < 3; ++i) v[i] = res.eigenvectors(i, col); + const Vector Av = A * v; + const double lam = eig.eigenvalues_real[col]; + double resid = 0.0; + for (std::size_t i = 0; i < 3; ++i) { + const double d = Av[i] - lam * v[i]; + resid += d * d; + } + REQUIRE(std::sqrt(resid) < 1e-6); + } +} + +TEST_CASE("eigenvectors_inverse_iteration: non-square throws", "[sym_eigen][inverse_iter]") { + Matrix A(2, 3); + Vector lambdas{1.0}; + REQUIRE_THROWS_AS(linalgebra::eigenvectors_inverse_iteration(A, lambdas), + DimensionMismatchError); +} -- cgit v1.2.3