diff options
| -rw-r--r-- | CMakeLists.txt | 21 | ||||
| -rw-r--r-- | README.md | 45 | ||||
| -rw-r--r-- | experiments/hilbert_qr.cpp | 119 | ||||
| -rw-r--r-- | experiments/matmul.cpp | 130 | ||||
| -rw-r--r-- | experiments/pivoting_vs_no_pivoting.cpp | 196 | ||||
| -rw-r--r-- | src/expm.cpp | 171 | ||||
| -rw-r--r-- | src/iterative.cpp | 404 | ||||
| -rw-r--r-- | src/linalgebra.cpp | 5 | ||||
| -rw-r--r-- | src/precond.cpp | 279 | ||||
| -rw-r--r-- | src/qr_iteration.cpp | 4 | ||||
| -rw-r--r-- | src/svd.cpp | 323 | ||||
| -rw-r--r-- | src/sym_eigen.cpp | 208 | ||||
| -rw-r--r-- | tests/test_expm.cpp | 144 | ||||
| -rw-r--r-- | tests/test_iterative.cpp | 211 | ||||
| -rw-r--r-- | tests/test_precond.cpp | 142 | ||||
| -rw-r--r-- | tests/test_qr.cpp | 2 | ||||
| -rw-r--r-- | tests/test_qr_iteration.cpp | 2 | ||||
| -rw-r--r-- | tests/test_svd.cpp | 140 | ||||
| -rw-r--r-- | tests/test_sym_eigen.cpp | 168 |
19 files changed, 2222 insertions, 492 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 70f03a1..d4e2faa 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -37,6 +37,11 @@ target_sources(linear_algebra src/qr.cpp src/qr_iteration.cpp src/cholesky.cpp + src/sym_eigen.cpp + src/svd.cpp + src/iterative.cpp + src/precond.cpp + src/expm.cpp ) target_compile_features(linear_algebra PUBLIC cxx_std_23) @@ -59,17 +64,6 @@ elseif(MSVC) endif() 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) - target_link_libraries(pivoting_vs_no_pivoting PRIVATE linear_algebra::core) - - 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) include(FetchContent) @@ -92,6 +86,11 @@ if(LINEAR_ALGEBRA_BUILD_TESTS) tests/test_qr.cpp tests/test_qr_iteration.cpp tests/test_cholesky.cpp + tests/test_sym_eigen.cpp + tests/test_svd.cpp + tests/test_iterative.cpp + tests/test_precond.cpp + tests/test_expm.cpp ) target_link_libraries(linear_algebra_tests @@ -4,7 +4,7 @@ I mostly follow Trefethen & Bau, "Numerical Linear Algebra" and Golub & Van Loan Computations." The implementation uses `NEON` SIMD on ARM64 systems when available. -## Build +### Build The library is packaged as a C++20 named module (`linalgebra`): @@ -14,8 +14,7 @@ The library is packaged as a C++20 named module (`linalgebra`): doesn't yet support module dependency scanning) ```bash -cmake -S . -B build -G Ninja \ - -DCMAKE_CXX_COMPILER=/opt/homebrew/opt/llvm/bin/clang++ +cmake -S . -B build -G Ninja -DCMAKE_CXX_COMPILER=/opt/homebrew/opt/llvm/bin/clang++ cmake --build build ``` @@ -29,7 +28,7 @@ cmake -S . -B build -G Ninja \ Valid values are `AUTO` (uses available SIMD) and `NONE` (forces scalar fallback). -## Usage +### Usage Import the module: @@ -44,13 +43,13 @@ int main() { } ``` -## Run tests +### Run tests ```bash ctest --test-dir build --output-on-failure ``` -## What's implemented +### What's implemented - Matrix / Vector core with SIMD matmul - Triangular solvers (forward / backward substitution) @@ -67,25 +66,15 @@ ctest --test-dir build --output-on-failure - Francis double-shift QR (`eigenvalues_francis`) — implicit bulge chasing on Hessenberg form; handles complex conjugate eigenvalue pairs without complex arithmetic; robust subdiagonal + 2×2 block deflation with exceptional shifts (GVL §7.5) -- Cholesky factorization (`cholesky_factor`, `cholesky_solve`) — for symmetric positive definite systems - -## TODO: -- [ ] Symmetric tridiagonalization — Householder reduction before symmetric QR (tridiagonalize) -- [ ] Eigenvectors via inverse iteration (eigenvectors_inverse_iteration) -- [ ] SVD — Golub-Kahan bidiagonalization + QR (svd) -- [ ] Conjugate Gradient (solve_cg) — for symmetric positive definite systems -- [ ] GMRES (solve_gmres) — for general non-symmetric systems -- [ ] BiCGSTAB (solve_bicgstab) — lighter alternative to GMRES -- [ ] Condition number estimation — norm-based LINPACK estimator -- [ ] Preconditioners (precond_jacobi, precond_ilu0) — diagonal and ILU(0) -- [ ] Least squares solver (lstsq) — via QR or SVD with rank-deficient handling -- [ ] Arnoldi iteration (arnoldi) — falls out naturally from GMRES -- [ ] Matrix exponential (expm) — via Padé approximation - -## Run experiments - -```bash -./build/matmul -./build/pivoting_vs_no_pivoting -./build/hilbert_qr -``` +- Cholesky factorization (`cholesky_factor`, `cholesky_solve`) +- Symmetric tridiagonalization +- Eigenvectors via inverse iteration (eigenvectors_inverse_iteration) +- SVD — Golub-Kahan bidiagonalization + QR (svd) +- Conjugate Gradient (solve_cg) +- GMRES (solve_gmres) +- BiCGSTAB (solve_bicgstab) +- Condition number estimation +- Preconditioners (precond_jacobi, precond_ilu0) +- Least squares solver (lstsq) +- Arnoldi iteration (arnoldi) +- Matrix exponential (expm) diff --git a/experiments/hilbert_qr.cpp b/experiments/hilbert_qr.cpp deleted file mode 100644 index 3a52194..0000000 --- a/experiments/hilbert_qr.cpp +++ /dev/null @@ -1,119 +0,0 @@ -import linalgebra; -import std; - -using linalgebra::Matrix; -using linalgebra::QRResult; - -Matrix hilbert(std::size_t n) { - Matrix H(n, n); - for (std::size_t i = 0; i < n; ++i) - for (std::size_t j = 0; j < n; ++j) - H(i, j) = 1.0 / static_cast<double>(i + j + 1); - return H; -} - -double reconstruction_error(const Matrix& A, const QRResult& qr) { - const std::size_t m = A.rows(); - const std::size_t n = A.cols(); - const Matrix QR = qr.Q * qr.R; - double err = 0.0; - for (std::size_t i = 0; i < m; ++i) - for (std::size_t j = 0; j < n; ++j) { - const double d = A(i, j) - QR(i, j); - err += d * d; - } - return std::sqrt(err); -} - -double orthogonality_error(const QRResult& qr) { - const Matrix& Q = qr.Q; - const std::size_t n = Q.cols(); - 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 < Q.rows(); ++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); -} - -using Clock = std::chrono::high_resolution_clock; -using Seconds = std::chrono::duration<double>; - -template<typename Fn> -double min_time(Fn fn, int trials = 5) { - double best = 1e18; - for (int t = 0; t < trials; ++t) { - const auto t0 = Clock::now(); - fn(); - const auto t1 = Clock::now(); - best = std::min(best, Seconds(t1 - t0).count()); - } - return best; -} - -using QRFn = std::function<QRResult(const Matrix&)>; - -struct Result { double recon, ortho, time_s; }; - -std::optional<Result> measure(const Matrix& A, QRFn fn) { - try { - const QRResult qr = fn(A); - const double re = reconstruction_error(A, qr); - const double oe = orthogonality_error(qr); - const double t = min_time([&] { fn(A); }); - return Result{re, oe, t}; - } catch (const std::exception&) { - return std::nullopt; - } -} - -void print_row(const std::string& method, std::optional<Result> r) { - std::cout << std::left << std::setw(16) << method; - if (!r) { - std::cout << " FAILED (linearly dependent columns)\n"; - return; - } - std::cout << std::scientific << std::setprecision(2) - << std::setw(14) << r->recon - << std::setw(14) << r->ortho - << std::fixed << std::setprecision(3) - << std::setw(10) << r->time_s * 1e6 << " µs\n"; -} - -int main() { - std::cout << std::string(70, '*') << "\n"; - std::cout << " Hilbert QR Experiment: comparing GS variants and Householder\n"; - std::cout << std::string(70, '*') << "\n\n"; - std::cout << - "H[i][j] = 1/(i+j+1). Condition number grows ~exponentially with n.\n" - "Orthogonality loss in classical GS tracks condition number directly.\n" - "Modified GS recovers ~half the lost digits. Householder is unaffected.\n\n"; - - const std::size_t sizes[] = {2, 3, 4, 5, 6, 7, 8, 10, 12}; - - for (std::size_t n : sizes) { - const Matrix H = hilbert(n); - - std::cout << std::string(70, '-') << "\n"; - std::cout << " n = " << n << "\n"; - std::cout << std::string(70, '-') << "\n"; - std::cout << std::left - << std::setw(16) << "Method" - << std::setw(14) << "||A-QR||_F" - << std::setw(14) << "||QtQ-I||_F" - << std::setw(10) << "Time\n"; - std::cout << std::string(70, ' ') << "\n"; - - print_row("classical_gs", - measure(H, [](const Matrix& A) { return linalgebra::qr_classical_gs(A); })); - print_row("modified_gs", - measure(H, [](const Matrix& A) { return linalgebra::qr_modified_gs(A); })); - print_row("householder", - measure(H, [](const Matrix& A) { return linalgebra::qr_householder(A); })); - } - - return 0; -} diff --git a/experiments/matmul.cpp b/experiments/matmul.cpp deleted file mode 100644 index a33b501..0000000 --- a/experiments/matmul.cpp +++ /dev/null @@ -1,130 +0,0 @@ -import linalgebra; -import std; - -#if !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 linalgebra::Matrix; -using Clock = std::chrono::high_resolution_clock; -using Seconds = std::chrono::duration<double>; - - -Matrix naive_matmul(const Matrix& lhs, const Matrix& rhs) { - if (lhs.cols() != rhs.rows()) { - throw linalgebra::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; -} - -Matrix make_matrix(std::size_t n) { - Matrix M(n, n); - const double inv = 1.0 / static_cast<double>(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<double>(i + j + 1) * inv; - return M; -} - -template <typename Fn> -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<double>(n); - return 2.0 * nd * nd * nd; -} - - -int main() { - const std::vector<std::size_t> 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"; - - 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 deleted file mode 100644 index 2f043cb..0000000 --- a/experiments/pivoting_vs_no_pivoting.cpp +++ /dev/null @@ -1,196 +0,0 @@ -import linalgebra; -import std; - -using linalgebra::Matrix; -using linalgebra::Vector; - -struct NoPivotLU { - Matrix L; - Matrix U; - bool failed = false; - std::size_t fail_step = 0; -}; - -NoPivotLU lu_no_pivot(const Matrix& A, double tol = 1e-14) { - const std::size_t n = A.rows(); - Matrix work = A; - Matrix L = Matrix::zeros(n, n); - for (std::size_t i = 0; i < n; ++i) L(i, i) = 1.0; - Matrix U = Matrix::zeros(n, n); - - for (std::size_t k = 0; k < n; ++k) { - if (std::abs(work(k, k)) <= tol) { - return NoPivotLU{std::move(L), std::move(U), true, k}; - } - for (std::size_t j = k; j < n; ++j) U(k, j) = work(k, j); - for (std::size_t i = k + 1; i < n; ++i) { - L(i, k) = work(i, k) / work(k, k); - for (std::size_t j = k + 1; j < n; ++j) { - work(i, j) -= L(i, k) * work(k, j); - } - } - } - return NoPivotLU{std::move(L), std::move(U), false, 0}; -} - -std::optional<Vector> solve_no_pivot(const NoPivotLU& f, const Vector& b) { - if (f.failed) return std::nullopt; - try { - const Vector y = linalgebra::forward_substitution(f.L, b, 1e-14, /*unit_diagonal=*/true); - return linalgebra::backward_substitution(f.U, y); - } catch (...) { - return std::nullopt; - } -} - -double solve_residual(const Matrix& A, const Vector& x, const Vector& b) { - return linalgebra::norm2(A * x - b); -} - -double reconstruction_error(const Matrix& A, const linalgebra::LUResult& lu) { - const std::size_t n = A.rows(); - Matrix PA(n, n); - for (std::size_t i = 0; i < n; ++i) - for (std::size_t j = 0; j < n; ++j) - PA(i, j) = A(lu.perm[i], j); - const Matrix LU_prod = lu.L * lu.U; - double err = 0.0; - for (std::size_t i = 0; i < n; ++i) - for (std::size_t j = 0; j < n; ++j) { - const double d = PA(i, j) - LU_prod(i, j); - err += d * d; - } - return std::sqrt(err); -} - -void print_header(const std::string& title) { - std::cout << "\n" << std::string(60, '=') << "\n"; - std::cout << " " << title << "\n"; - std::cout << std::string(60, '=') << "\n"; - std::cout << std::left - << std::setw(22) << "Method" - << std::setw(20) << "||Ax - b||" - << std::setw(20) << "||PA - LU||" - << "\n"; - std::cout << std::string(60, '-') << "\n"; -} - -void report_pivoted(const Matrix& A, const Vector& b) { - try { - const linalgebra::LUResult lu = linalgebra::lu_factor(A); - const Vector x = linalgebra::lu_solve(lu, b); - std::cout << std::left << std::setw(22) << "Pivoted LU" - << std::setw(20) << std::scientific << std::setprecision(3) - << solve_residual(A, x, b) - << std::setw(20) << reconstruction_error(A, lu) - << "\n"; - } catch (const std::exception& e) { - std::cout << std::left << std::setw(22) << "Pivoted LU" - << "FAILED: " << e.what() << "\n"; - } -} - -void report_no_pivot(const Matrix& A, const Vector& b) { - const NoPivotLU f = lu_no_pivot(A); - if (f.failed) { - std::cout << std::left << std::setw(22) << "No-pivot LU" - << "FAILED at step " << f.fail_step << " (zero pivot)\n"; - return; - } - const auto x_opt = solve_no_pivot(f, b); - if (!x_opt) { - std::cout << std::left << std::setw(22) << "No-pivot LU" - << "FAILED during solve (singular U)\n"; - return; - } - const Matrix LU_prod = f.L * f.U; - double rec_err = 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) - LU_prod(i, j); - rec_err += d * d; - } - rec_err = std::sqrt(rec_err); - - std::cout << std::left << std::setw(22) << "No-pivot LU" - << std::setw(20) << std::scientific << std::setprecision(3) - << solve_residual(A, *x_opt, b) - << std::setw(20) << rec_err - << "\n"; -} - -void run_case(const std::string& label, const Matrix& A, const Vector& b) { - print_header(label); - report_pivoted(A, b); - report_no_pivot(A, b); -} - -void exp_random(std::size_t n = 8) { - std::mt19937 rng(42); - std::uniform_real_distribution<double> dist(-5.0, 5.0); - Matrix A(n, n); - for (std::size_t i = 0; i < n; ++i) - for (std::size_t j = 0; j < n; ++j) - A(i, j) = dist(rng); - - Vector b(n); - for (std::size_t i = 0; i < n; ++i) b[i] = dist(rng); - - run_case("Random 8x8 (well-conditioned)", A, b); -} - -void exp_badly_scaled() { - const Matrix A{ - {1e-14, 1.0, 2.0 }, - {1.0, 3.0, 4.0 }, - {2.0, 5.0, 7.0 } - }; - const Vector b{1e-14 + 3.0, 8.0, 14.0}; - run_case("Badly scaled (row norms differ by 10^14)", A, b); -} - -void exp_epsilon_pathology() { - constexpr double eps = 1e-15; - const Matrix A{{eps, 1.0}, {1.0, 2.0}}; - 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"; -} - -void exp_amplified_multiplier() { - const Matrix A{ - {0.001, 1.0, 0.0, 0.0 }, - {1.0, 2.0, 1.0, 0.0 }, - {0.0, 1.0, 3.0, 1.0 }, - {0.0, 0.0, 1.0, 4.0 } - }; - const Vector b = A * Vector{1.0, 2.0, 3.0, 4.0}; - run_case("Amplified multiplier (small (1,1) pivot, 4x4)", A, b); - std::cout << " Note: exact solution is x = [1, 2, 3, 4]\n"; -} - -void exp_permutation() { - const Matrix A{ - {0.0, 0.0, 3.0}, - {0.0, 2.0, 1.0}, - {5.0, 1.0, 0.0} - }; - const Vector b = A * Vector{1.0, -1.0, 2.0}; - run_case("Multiple row swaps required (zeros in pivot positions)", A, b); -} - -int main() { - std::cout << std::string(60, '*') << "\n"; - std::cout << " Pivoting vs No-Pivoting LU Experiment\n"; - std::cout << std::string(60, '*') << "\n"; - std::cout << "Residual = ||Ax - b||_2 (solve accuracy)\n"; - std::cout << "Recon err = ||PA - LU||_F (factorization accuracy)\n"; - - exp_random(); - exp_badly_scaled(); - exp_epsilon_pathology(); - exp_amplified_multiplier(); - exp_permutation(); - - return 0; -} diff --git a/src/expm.cpp b/src/expm.cpp new file mode 100644 index 0000000..71bdbd9 --- /dev/null +++ b/src/expm.cpp @@ -0,0 +1,171 @@ +export module linalgebra:expm; +import std; +import :error; +import :vector; +import :matrix; +import :lu; + +// References used throughout this file: +// Higham — "Functions of Matrices: Theory and Computation" (SIAM 2008) +// Al-Mohy & Higham (2009), SIAM J. Matrix Anal. Appl. 31(3) + +export namespace linalgebra { + +// Matrix exponential via Padé [13/13] approximant + scaling and squaring. +[[nodiscard]] Matrix expm(const Matrix& A); + +} // namespace linalgebra + +namespace { + +// Padé [13/13] coefficients, Higham (2008) Table 10.4. +constexpr double pade_b[14] = { + 64764752532480000.0, + 32382376266240000.0, + 7771770303897600.0, + 1187353796428800.0, + 129060195264000.0, + 10559470521600.0, + 670442572800.0, + 33522128640.0, + 1323241920.0, + 40840800.0, + 960960.0, + 16380.0, + 182.0, + 1.0, +}; + +// theta_13: ||A||_1 threshold below which no scaling is needed. +constexpr double theta_13 = 5.371920351148152; + +// 1-norm of a matrix. +double one_norm(const linalgebra::Matrix& M) { + const std::size_t n = M.cols(); + const std::size_t m = M.rows(); + double result = 0.0; + for (std::size_t j = 0; j < n; ++j) { + double col_sum = 0.0; + for (std::size_t i = 0; i < m; ++i) col_sum += std::abs(M(i, j)); + result = std::max(result, col_sum); + } + return result; +} + +// Evaluate the Padé [13/13] numerator U and denominator V for matrix B. +// Uses the factored evaluation from Higham Algorithm 10.20. +// +// W1 = b[13]*A6 + b[11]*A4 + b[9]*A2 +// W2 = b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*I +// Z1 = b[12]*A6 + b[10]*A4 + b[8]*A2 +// Z2 = b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*I +// W = A6*W1 + W2 +// U = B * W +// V = A6*Z1 + Z2 +// +// expm(B) ≈ (V - U)^{-1} * (V + U) +std::pair<linalgebra::Matrix, linalgebra::Matrix> +pade13(const linalgebra::Matrix& B) { + const std::size_t n = B.rows(); + const linalgebra::Matrix I = linalgebra::Matrix::identity(n); + + const linalgebra::Matrix A2 = B * B; + const linalgebra::Matrix A4 = A2 * A2; + const linalgebra::Matrix A6 = A2 * A4; + + // W1 = b[13]*A6 + b[11]*A4 + b[9]*A2 + linalgebra::Matrix W1(n, n, 0.0); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + W1(i, j) = pade_b[13] * A6(i, j) + pade_b[11] * A4(i, j) + pade_b[9] * A2(i, j); + + // W2 = b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*I + linalgebra::Matrix W2(n, n, 0.0); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + W2(i, j) = pade_b[7] * A6(i, j) + pade_b[5] * A4(i, j) + + pade_b[3] * A2(i, j) + pade_b[1] * I(i, j); + + // Z1 = b[12]*A6 + b[10]*A4 + b[8]*A2 + linalgebra::Matrix Z1(n, n, 0.0); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + Z1(i, j) = pade_b[12] * A6(i, j) + pade_b[10] * A4(i, j) + pade_b[8] * A2(i, j); + + // Z2 = b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*I + linalgebra::Matrix Z2(n, n, 0.0); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + Z2(i, j) = pade_b[6] * A6(i, j) + pade_b[4] * A4(i, j) + + pade_b[2] * A2(i, j) + pade_b[0] * I(i, j); + + // W = A6*W1 + W2 + const linalgebra::Matrix W = A6 * W1 + W2; + + // U = B * W (numerator) + const linalgebra::Matrix U = B * W; + + // V = A6*Z1 + Z2 (denominator) + const linalgebra::Matrix V = A6 * Z1 + Z2; + + return {U, V}; +} + +} // namespace + +namespace linalgebra { + +Matrix expm(const Matrix& A) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "expm requires a square matrix, got " << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + + const std::size_t n = A.rows(); + if (n == 0) return Matrix::identity(0); + + // Determine scaling factor s such that ||A / 2^s||_1 <= theta_13. + const double norm_A = one_norm(A); + int s = 0; + if (norm_A > theta_13) { + s = static_cast<int>(std::ceil(std::log2(norm_A / theta_13))); + if (s < 0) s = 0; + } + + // Scale B = A / 2^s. + const double scale = 1.0 / std::ldexp(1.0, s); + Matrix B(n, n, 0.0); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + B(i, j) = A(i, j) * scale; + + // Compute Padé [13/13] approximant: R = (V - U)^{-1} * (V + U). + auto [U, V] = pade13(B); + + // Numerator = V + U, Denominator = V - U. + Matrix Numerator(n, n, 0.0); + Matrix Denominator(n, n, 0.0); + for (std::size_t i = 0; i < n; ++i) { + for (std::size_t j = 0; j < n; ++j) { + Numerator(i, j) = V(i, j) + U(i, j); + Denominator(i, j) = V(i, j) - U(i, j); + } + } + + const LUResult lu_denom = lu_factor(Denominator); + + Matrix R(n, n, 0.0); + for (std::size_t j = 0; j < n; ++j) { + Vector col(n); + for (std::size_t i = 0; i < n; ++i) col[i] = Numerator(i, j); + const Vector sol = lu_solve(lu_denom, col); + for (std::size_t i = 0; i < n; ++i) R(i, j) = sol[i]; + } + + for (int i = 0; i < s; ++i) R = R * R; + + return R; +} + +} // namespace linalgebra diff --git a/src/iterative.cpp b/src/iterative.cpp new file mode 100644 index 0000000..a0ff41d --- /dev/null +++ b/src/iterative.cpp @@ -0,0 +1,404 @@ +export module linalgebra:iterative; +import std; +import :error; +import :vector; +import :matrix; +import :norms; +import :triangular_solve; + +// References used throughout this file: +// T&B — Trefethen & Bau, "Numerical Linear Algebra" +// GVL — Golub & Van Loan, "Matrix Computations" 4th ed. + +export namespace linalgebra { + +// Arnoldi iteration builds an orthonormal Krylov basis and the +// corresponding upper Hessenberg matrix. +// Reference: T&B Algorithm 33.1; GVL §6.3 + +struct ArnoldiResult { + Matrix Q; // n × (steps_taken + 1) orthonormal columns + Matrix H; // (steps_taken + 1) × steps_taken upper Hessenberg + int steps_taken; + bool breakdown; // true if invariant subspace found early +}; + +struct ArnoldiOptions { + double breakdown_tolerance = 1e-14; +}; + +[[nodiscard]] ArnoldiResult arnoldi(const Matrix& A, const Vector& b, int k, + ArnoldiOptions opts = {}); + +// Conjugate Gradient — for symmetric positive definite systems +// Reference: T&B Algorithm 38.1; GVL §11.3 + +struct CGOptions { + double tolerance = 1e-10; + int max_iterations = 1000; +}; + +struct CGResult { + Vector x; + int iterations; + double final_residual; +}; + +[[nodiscard]] CGResult solve_cg(const Matrix& A, const Vector& b, CGOptions opts = {}); + +// Restarted GMRES — for general square systems +// Reference: T&B Algorithm 35.1; GVL §11.4.2 + +struct GMRESOptions { + double tolerance = 1e-10; + int max_iterations = 200; + int restart = 50; +}; + +struct GMRESResult { + Vector x; + int iterations; + double final_residual; +}; + +[[nodiscard]] GMRESResult solve_gmres(const Matrix& A, const Vector& b, + GMRESOptions opts = {}); + +// BiCGSTAB — for general square systems (van der Vorst 1992) +// Reference: GVL §11.5.3 + +struct BiCGSTABOptions { + double tolerance = 1e-10; + int max_iterations = 1000; +}; + +struct BiCGSTABResult { + Vector x; + int iterations; + double final_residual; +}; + +[[nodiscard]] BiCGSTABResult solve_bicgstab(const Matrix& A, const Vector& b, + BiCGSTABOptions opts = {}); + +} // namespace linalgebra + +namespace linalgebra { + +ArnoldiResult arnoldi(const Matrix& A, const Vector& b, int k, ArnoldiOptions opts) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "arnoldi requires a square matrix, got " << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + const std::size_t n = A.rows(); + if (b.size() != n) { + std::ostringstream oss; + oss << "arnoldi: b size " << b.size() << " must match matrix dimension " << n; + throw DimensionMismatchError(oss.str()); + } + if (k <= 0) throw std::invalid_argument("arnoldi: k must be >= 1"); + + const auto kk = static_cast<std::size_t>(k); + + Matrix Q(n, kk + 1, 0.0); + Matrix H(kk + 1, kk, 0.0); + + const double b_norm = norm2(b); + if (b_norm == 0.0) { + return ArnoldiResult{std::move(Q), std::move(H), 0, true}; + } + + // q_0 = b / ||b|| + for (std::size_t i = 0; i < n; ++i) Q(i, 0) = b[i] / b_norm; + + int steps = 0; + for (std::size_t j = 0; j < kk; ++j) { + // z = A * Q[:, j] + Vector qj(n); + for (std::size_t i = 0; i < n; ++i) qj[i] = Q(i, j); + Vector z = A * qj; + + // Modified Gram-Schmidt orthogonalization. + for (std::size_t i = 0; i <= j; ++i) { + double h = 0.0; + for (std::size_t row = 0; row < n; ++row) h += Q(row, i) * z[row]; + H(i, j) = h; + for (std::size_t row = 0; row < n; ++row) z[row] -= h * Q(row, i); + } + + const double z_norm = norm2(z); + H(j + 1, j) = z_norm; + ++steps; + + if (z_norm < opts.breakdown_tolerance) { + // Lucky breakdown: invariant subspace found. + return ArnoldiResult{std::move(Q), std::move(H), steps, true}; + } + + for (std::size_t i = 0; i < n; ++i) Q(i, j + 1) = z[i] / z_norm; + } + + return ArnoldiResult{std::move(Q), std::move(H), steps, false}; +} + +CGResult solve_cg(const Matrix& A, const Vector& b, CGOptions opts) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "solve_cg requires a square matrix, got " << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + const std::size_t n = A.rows(); + if (b.size() != n) { + std::ostringstream oss; + oss << "solve_cg: rhs size " << b.size() << " does not match dimension " << n; + throw DimensionMismatchError(oss.str()); + } + + // Check symmetry. + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = i + 1; j < n; ++j) + if (std::abs(A(i, j) - A(j, i)) > 1e-10 * (std::abs(A(i, j)) + 1.0)) + throw LinAlgError("solve_cg: matrix is not symmetric"); + + Vector x(n, 0.0); + Vector r = b; // r = b - A*x0, x0 = 0 + Vector p = r; + double rr = dot(r, r); + + if (std::sqrt(rr) < opts.tolerance) { + return CGResult{std::move(x), 0, std::sqrt(rr)}; + } + + for (int iter = 1; iter <= opts.max_iterations; ++iter) { + const Vector Ap = A * p; + const double pAp = dot(p, Ap); + if (std::abs(pAp) == 0.0) break; + const double alpha = rr / pAp; + + for (std::size_t i = 0; i < n; ++i) x[i] += alpha * p[i]; + for (std::size_t i = 0; i < n; ++i) r[i] -= alpha * Ap[i]; + + const double rr_new = dot(r, r); + const double res = std::sqrt(rr_new); + + if (res < opts.tolerance) { + return CGResult{std::move(x), iter, res}; + } + + const double beta = rr_new / rr; + for (std::size_t i = 0; i < n; ++i) p[i] = r[i] + beta * p[i]; + rr = rr_new; + } + + const double final_res = norm2(A * x - b); + std::ostringstream oss; + oss << "solve_cg: did not converge in " << opts.max_iterations + << " iterations (final residual = " << final_res << ")"; + throw NonConvergenceError(oss.str()); +} + +GMRESResult solve_gmres(const Matrix& A, const Vector& b, GMRESOptions opts) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "solve_gmres requires a square matrix, got " << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + const std::size_t n = A.rows(); + if (b.size() != n) { + std::ostringstream oss; + oss << "solve_gmres: rhs size " << b.size() << " does not match dimension " << n; + throw DimensionMismatchError(oss.str()); + } + + const int m = std::min(opts.restart, static_cast<int>(n)); + + Vector x(n, 0.0); + int total_iters = 0; + + while (total_iters < opts.max_iterations) { + // Compute residual r = b - A*x. + const Vector Ax = A * x; + Vector r(n); + for (std::size_t i = 0; i < n; ++i) r[i] = b[i] - Ax[i]; + const double beta = norm2(r); + + if (beta < opts.tolerance) { + return GMRESResult{std::move(x), total_iters, beta}; + } + + // Arnoldi to build Krylov basis. + const ArnoldiResult ar = arnoldi(A, r, m); + const int steps = ar.steps_taken; + + if (steps == 0) break; + + // Solve least-squares problem: min ||beta*e1 - H_hat * y|| + // where H_hat is (steps+1) x steps. + // Apply Givens rotations to reduce H_hat to upper triangular. + const auto sz = static_cast<std::size_t>(steps); + + // Work on a copy of the relevant submatrix of H and rhs g. + std::vector<std::vector<double>> Hwork(sz + 1, std::vector<double>(sz, 0.0)); + for (std::size_t i = 0; i <= sz; ++i) + for (std::size_t j = 0; j < sz; ++j) + Hwork[i][j] = ar.H(i, j); + + std::vector<double> g(sz + 1, 0.0); + g[0] = beta; + + // Accumulated Givens rotations. + std::vector<double> cs(sz), sn(sz); + + for (std::size_t j = 0; j < sz; ++j) { + // Givens to zero H[j+1, j]. + const double f = Hwork[j][j]; + const double hh = Hwork[j + 1][j]; + const double r_val = std::hypot(f, hh); + if (r_val == 0.0) { cs[j] = 1.0; sn[j] = 0.0; continue; } + cs[j] = f / r_val; + sn[j] = hh / r_val; + + // Apply to column j of H (only the two relevant rows). + Hwork[j][j] = cs[j] * f + sn[j] * hh; + Hwork[j + 1][j] = 0.0; + + // Apply to remaining columns. + for (std::size_t l = j + 1; l < sz; ++l) { + const double t0 = Hwork[j][l]; + const double t1 = Hwork[j + 1][l]; + Hwork[j][l] = cs[j] * t0 + sn[j] * t1; + Hwork[j + 1][l] = -sn[j] * t0 + cs[j] * t1; + } + + // Apply to g. + const double g0 = g[j]; + const double g1 = g[j + 1]; + g[j] = cs[j] * g0 + sn[j] * g1; + g[j + 1] = -sn[j] * g0 + cs[j] * g1; + } + + // Backward substitution: solve the sz×sz upper triangular system. + std::vector<double> y(sz, 0.0); + for (int i = static_cast<int>(sz) - 1; i >= 0; --i) { + double sum = g[static_cast<std::size_t>(i)]; + for (std::size_t j = static_cast<std::size_t>(i) + 1; j < sz; ++j) + sum -= Hwork[static_cast<std::size_t>(i)][j] * y[j]; + if (std::abs(Hwork[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)]) > 0.0) + y[static_cast<std::size_t>(i)] = + sum / Hwork[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)]; + } + + // Update x = x + Q[:, 0:steps] * y. + for (std::size_t j = 0; j < sz; ++j) { + for (std::size_t i = 0; i < n; ++i) { + x[i] += y[j] * ar.Q(i, j); + } + } + + ++total_iters; + + // Check convergence. + const double final_res = std::abs(g[sz]); + if (final_res < opts.tolerance || ar.breakdown) { + return GMRESResult{std::move(x), total_iters, final_res}; + } + } + + const double final_res = norm2(A * x - b); + std::ostringstream oss; + oss << "solve_gmres: did not converge in " << opts.max_iterations + << " restarts (final residual = " << final_res << ")"; + throw NonConvergenceError(oss.str()); +} + +BiCGSTABResult solve_bicgstab(const Matrix& A, const Vector& b, BiCGSTABOptions opts) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "solve_bicgstab requires a square matrix, got " + << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + const std::size_t n = A.rows(); + if (b.size() != n) { + std::ostringstream oss; + oss << "solve_bicgstab: rhs size " << b.size() << " does not match dimension " << n; + throw DimensionMismatchError(oss.str()); + } + + Vector x(n, 0.0); + Vector r = b; // r = b - A*x0, x0 = 0 + Vector r_hat = r; // shadow residual, fixed throughout + + double rho_old = 1.0, alpha = 1.0, omega = 1.0; + Vector v(n, 0.0), p(n, 0.0); + + for (int iter = 1; iter <= opts.max_iterations; ++iter) { + const double rho_new = dot(r_hat, r); + + if (std::abs(rho_new) < 1e-300) { + std::ostringstream oss; + oss << "solve_bicgstab: breakdown (rho near zero) at iteration " << iter; + throw NonConvergenceError(oss.str()); + } + + const double beta = (rho_new / rho_old) * (alpha / omega); + + for (std::size_t i = 0; i < n; ++i) + p[i] = r[i] + beta * (p[i] - omega * v[i]); + + v = A * p; + + const double denom_alpha = dot(r_hat, v); + if (std::abs(denom_alpha) < 1e-300) { + std::ostringstream oss; + oss << "solve_bicgstab: breakdown (r_hat·v near zero) at iteration " << iter; + throw NonConvergenceError(oss.str()); + } + alpha = rho_new / denom_alpha; + + Vector s(n); + for (std::size_t i = 0; i < n; ++i) s[i] = r[i] - alpha * v[i]; + + const double s_norm = norm2(s); + if (s_norm < opts.tolerance) { + for (std::size_t i = 0; i < n; ++i) x[i] += alpha * p[i]; + return BiCGSTABResult{std::move(x), iter, s_norm}; + } + + const Vector t = A * s; + const double tt = dot(t, t); + + if (std::abs(tt) < 1e-300) { + std::ostringstream oss; + oss << "solve_bicgstab: breakdown (t·t near zero) at iteration " << iter; + throw NonConvergenceError(oss.str()); + } + + omega = dot(t, s) / tt; + + if (std::abs(omega) < 1e-300) { + std::ostringstream oss; + oss << "solve_bicgstab: breakdown (omega near zero) at iteration " << iter; + throw NonConvergenceError(oss.str()); + } + + for (std::size_t i = 0; i < n; ++i) x[i] += alpha * p[i] + omega * s[i]; + for (std::size_t i = 0; i < n; ++i) r[i] = s[i] - omega * t[i]; + + const double r_norm = norm2(r); + if (r_norm < opts.tolerance) { + return BiCGSTABResult{std::move(x), iter, r_norm}; + } + + rho_old = rho_new; + } + + const double final_res = norm2(A * x - b); + std::ostringstream oss; + oss << "solve_bicgstab: did not converge in " << opts.max_iterations + << " iterations (final residual = " << final_res << ")"; + throw NonConvergenceError(oss.str()); +} + +} // namespace linalgebra diff --git a/src/linalgebra.cpp b/src/linalgebra.cpp index 18e9cf3..d6d5652 100644 --- a/src/linalgebra.cpp +++ b/src/linalgebra.cpp @@ -9,3 +9,8 @@ export import :lu; export import :qr; export import :qr_iteration; export import :cholesky; +export import :sym_eigen; +export import :svd; +export import :iterative; +export import :precond; +export import :expm; diff --git a/src/precond.cpp b/src/precond.cpp new file mode 100644 index 0000000..a4ed503 --- /dev/null +++ b/src/precond.cpp @@ -0,0 +1,279 @@ +export module linalgebra:precond; +import std; +import :error; +import :vector; +import :matrix; +import :norms; +import :triangular_solve; +import :lu; +import :qr; + +// References used throughout this file: +// T&B — Trefethen & Bau, "Numerical Linear Algebra" +// GVL — Golub & Van Loan, "Matrix Computations" 4th ed. +// Hager — Hager, W.W. (1984), SIAM J. Sci. Stat. Comput. 5(2):311-316 + +export namespace linalgebra { + +// Condition number estimation (1-norm, Hager/LINPACK power iteration) +// Reference: Hager (1984); GVL §2.3.3 + +[[nodiscard]] double condition_number_1norm(const Matrix& A, + double singular_tolerance = 1e-12); + +// Jacobi (diagonal) preconditioner + +struct JacobiPrecond { + Vector inv_diag; +}; + +[[nodiscard]] JacobiPrecond precond_jacobi(const Matrix& A, + double zero_tolerance = 1e-14); + +[[nodiscard]] Vector apply(const JacobiPrecond& P, const Vector& x); + +// ILU(0) preconditioner (for dense matrices = LU without pivoting) +// Reference: Saad, "Iterative Methods for Sparse Linear Systems" §10.3 + +struct ILU0Precond { + Matrix LU; // combined: strict lower = L multipliers, upper = U +}; + +[[nodiscard]] ILU0Precond precond_ilu0(const Matrix& A, + double zero_tolerance = 1e-14); + +[[nodiscard]] Vector apply(const ILU0Precond& P, const Vector& b); + +// Least-squares solver via column-pivoting QR +// Reference: T&B Lecture 11; GVL §5.5 + +struct LstsqResult { + Vector x; + std::size_t rank; + double residual_norm; +}; + +struct LstsqOptions { + double rank_tolerance = 1e-12; +}; + +[[nodiscard]] LstsqResult lstsq(const Matrix& A, const Vector& b, + LstsqOptions opts = {}); + +} // namespace linalgebra + +namespace { + +// Solve A^T z = rhs given a pre-computed LU factorization of A. +// PA = LU => A^T = U^T L^T P +// Steps: (1) solve U^T q = rhs, (2) solve L^T w = q, (3) z[perm[i]] = w[i] +linalgebra::Vector solve_transpose(const linalgebra::LUResult& lu, + const linalgebra::Vector& rhs) { + const std::size_t n = lu.L.rows(); + + // Solve U^T q = rhs (U^T is lower triangular) + const linalgebra::Matrix Ut = linalgebra::transpose(lu.U); + const linalgebra::Vector q = linalgebra::forward_substitution(Ut, rhs); + + // Solve L^T w = q (L^T is upper triangular, unit diagonal) + const linalgebra::Matrix Lt = linalgebra::transpose(lu.L); + const linalgebra::Vector w = linalgebra::backward_substitution(Lt, q, 1e-14, true); + + // Apply inverse permutation: z[perm[i]] = w[i] + linalgebra::Vector z(n); + for (std::size_t i = 0; i < n; ++i) z[lu.perm[i]] = w[i]; + return z; +} + +} // namespace + +namespace linalgebra { + +double condition_number_1norm(const Matrix& A, double singular_tolerance) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "condition_number_1norm requires a square matrix, got " + << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + const std::size_t n = A.rows(); + + // Exact 1-norm of A: max column sum of absolute values. + double norm_A = 0.0; + for (std::size_t j = 0; j < n; ++j) { + double col_sum = 0.0; + for (std::size_t i = 0; i < n; ++i) col_sum += std::abs(A(i, j)); + norm_A = std::max(norm_A, col_sum); + } + + const LUResult lu = lu_factor(A, singular_tolerance); + + // Estimate ||A^{-1}||_1 via the power-iteration method (Hager 1984). + // Start with x = [1/n, ..., 1/n]. + Vector x(n, 1.0 / static_cast<double>(n)); + double est = 0.0; + + for (int iter = 0; iter < 5; ++iter) { + const Vector y = lu_solve(lu, x); // y = A^{-1} x + + // 1-norm of y. + double y1 = 0.0; + for (std::size_t i = 0; i < n; ++i) y1 += std::abs(y[i]); + + if (y1 <= est) break; + est = y1; + + Vector xi(n); + for (std::size_t i = 0; i < n; ++i) xi[i] = (y[i] >= 0.0) ? 1.0 : -1.0; + + // z = A^{-T} xi + const Vector z = solve_transpose(lu, xi); + + // Find the index maximizing |z[j]|. + std::size_t j_max = 0; + double max_z = std::abs(z[0]); + for (std::size_t i = 1; i < n; ++i) { + if (std::abs(z[i]) > max_z) { + max_z = std::abs(z[i]); + j_max = i; + } + } + + // Convergence check. + double xz = 0.0; + for (std::size_t i = 0; i < n; ++i) xz += std::abs(z[i]) / static_cast<double>(n); + if (max_z <= xz) break; + + // New starting vector: e_{j_max}. + x.fill(0.0); + x[j_max] = 1.0; + } + + return norm_A * est; +} + +JacobiPrecond precond_jacobi(const Matrix& A, double zero_tolerance) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "precond_jacobi requires a square matrix, got " + << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + const std::size_t n = A.rows(); + Vector inv_diag(n); + for (std::size_t i = 0; i < n; ++i) { + if (std::abs(A(i, i)) <= zero_tolerance) { + std::ostringstream oss; + oss << "precond_jacobi: zero diagonal entry at index " << i; + throw SingularMatrixError(oss.str()); + } + inv_diag[i] = 1.0 / A(i, i); + } + return JacobiPrecond{std::move(inv_diag)}; +} + +Vector apply(const JacobiPrecond& P, const Vector& x) { + const std::size_t n = x.size(); + if (n != P.inv_diag.size()) { + throw DimensionMismatchError("apply(JacobiPrecond): size mismatch"); + } + Vector result(n); + for (std::size_t i = 0; i < n; ++i) result[i] = P.inv_diag[i] * x[i]; + return result; +} + +ILU0Precond precond_ilu0(const Matrix& A, double zero_tolerance) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "precond_ilu0 requires a square matrix, got " + << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + const std::size_t n = A.rows(); + Matrix LU = A; + + for (std::size_t k = 0; k < n; ++k) { + if (std::abs(LU(k, k)) <= zero_tolerance) { + std::ostringstream oss; + oss << "precond_ilu0: near-zero pivot at step " << k; + throw SingularMatrixError(oss.str()); + } + for (std::size_t i = k + 1; i < n; ++i) { + LU(i, k) /= LU(k, k); + for (std::size_t j = k + 1; j < n; ++j) { + LU(i, j) -= LU(i, k) * LU(k, j); + } + } + } + return ILU0Precond{std::move(LU)}; +} + +Vector apply(const ILU0Precond& P, const Vector& b) { + const std::size_t n = P.LU.rows(); + if (b.size() != n) { + throw DimensionMismatchError("apply(ILU0Precond): size mismatch"); + } + + // Extract L (unit lower) and U (upper) from combined storage. + Matrix L = Matrix::zeros(n, n); + Matrix U = Matrix::zeros(n, n); + for (std::size_t i = 0; i < n; ++i) { + L(i, i) = 1.0; + for (std::size_t j = 0; j < i; ++j) L(i, j) = P.LU(i, j); + for (std::size_t j = i; j < n; ++j) U(i, j) = P.LU(i, j); + } + + const Vector y = forward_substitution(L, b, 1e-14, true); + return backward_substitution(U, y); +} + +LstsqResult lstsq(const Matrix& A, const Vector& b, LstsqOptions opts) { + const std::size_t m = A.rows(); + const std::size_t n = A.cols(); + + if (m < n) { + std::ostringstream oss; + oss << "lstsq requires rows >= cols, got " << m << "x" << n; + throw DimensionMismatchError(oss.str()); + } + if (b.size() != m) { + std::ostringstream oss; + oss << "lstsq: rhs size " << b.size() << " does not match rows " << m; + throw DimensionMismatchError(oss.str()); + } + + const QRColPivResult qr = qr_colpiv(A, opts.rank_tolerance); + const std::size_t r = qr.rank; + + // c = Q^T b (Q is m×n with orthonormal columns) + Vector c(n, 0.0); + for (std::size_t j = 0; j < n; ++j) { + double d = 0.0; + for (std::size_t i = 0; i < m; ++i) d += qr.Q(i, j) * b[i]; + c[j] = d; + } + + Vector x(n, 0.0); + + if (r > 0) { + Matrix Rr(r, r); + for (std::size_t i = 0; i < r; ++i) + for (std::size_t j = 0; j < r; ++j) + Rr(i, j) = qr.R(i, j); + + Vector cr(r); + for (std::size_t i = 0; i < r; ++i) cr[i] = c[i]; + + const Vector y = backward_substitution(Rr, cr); + + // Permuted solution: x_perm[0:r] = y, x_perm[r:n] = 0, then un-permute. + for (std::size_t j = 0; j < r; ++j) x[qr.perm[j]] = y[j]; + } + + const Vector res = A * x - b; + const double residual_norm = norm2(res); + + return LstsqResult{std::move(x), r, residual_norm}; +} + +} // namespace linalgebra diff --git a/src/qr_iteration.cpp b/src/qr_iteration.cpp index f8a3be5..7f50097 100644 --- a/src/qr_iteration.cpp +++ b/src/qr_iteration.cpp @@ -14,9 +14,7 @@ import :qr; export namespace linalgebra { -// --------------------------------------------------------------------------- // Options -// --------------------------------------------------------------------------- struct QRIterationOptions { double tolerance = 1e-10; @@ -509,9 +507,7 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A, QRIterationOptions opt return result; } -// --------------------------------------------------------------------------- // Francis double-shift QR with implicit bulge chasing -// --------------------------------------------------------------------------- QRIterationResult eigenvalues_francis(const Matrix& A, QRIterationOptions opts) { require_square(A, "eigenvalues_francis"); diff --git a/src/svd.cpp b/src/svd.cpp new file mode 100644 index 0000000..8157845 --- /dev/null +++ b/src/svd.cpp @@ -0,0 +1,323 @@ +export module linalgebra:svd; +import std; +import :error; +import :vector; +import :matrix; +import :qr_iteration; + +// References used throughout this file: +// GVL — Golub & Van Loan, "Matrix Computations" 4th ed. +// T&B — Trefethen & Bau, "Numerical Linear Algebra" + +export namespace linalgebra { + +struct SVDResult { + Matrix U; // m × m orthogonal (left singular vectors) + Vector sigma; // min(m,n) singular values, sorted descending + Matrix Vt; // n × n orthogonal (V^T, right singular vectors transposed) +}; + +struct SVDOptions { + double tolerance = 1e-12; + int max_iterations = 1000; +}; + +// Golub-Kahan bidiagonalization + Golub-Reinsch QR sweeps. +// Requires rows >= cols; throws DimensionMismatchError otherwise. +// Reference: GVL §8.6; T&B Lecture 31. +[[nodiscard]] SVDResult svd(const Matrix& A, SVDOptions opts = {}); + +} // namespace linalgebra + +namespace linalgebra { + +SVDResult svd(const Matrix& A, SVDOptions opts) { + const std::size_t m = A.rows(); + const std::size_t n = A.cols(); + + if (m < n) { + std::ostringstream oss; + oss << "svd requires rows >= cols, got " << m << "x" << n; + throw DimensionMismatchError(oss.str()); + } + if (n == 0) { + return SVDResult{Matrix::identity(m), Vector(0), Matrix::identity(0)}; + } + + // GK bidiagonalization A = U_acc * B * Vt_acc + // where B is m×n upper bidiagonal (nonzero on diagonal and superdiagonal). + + Matrix work = A; + Matrix U_acc = Matrix::identity(m); + Matrix Vt_acc = Matrix::identity(n); + + for (std::size_t k = 0; k < n; ++k) { + const std::size_t p_left = m - k; // column length below row k + + // Left Householder: zero work[k+1:, k]. + { + std::vector<double> u(p_left); + for (std::size_t i = 0; i < p_left; ++i) u[i] = work(k + i, k); + + double x_norm = 0.0; + for (double v : u) x_norm += v * v; + x_norm = std::sqrt(x_norm); + + if (x_norm > 0.0) { + const double sigma = (u[0] >= 0.0 ? 1.0 : -1.0) * x_norm; + u[0] += sigma; + double utu = 0.0; + for (double v : u) utu += v * v; + const double tau = 2.0 / utu; + + // Apply to work from left: work[k:, k:] -= tau * u * (u^T work[k:, k:]) + for (std::size_t j = k; j < n; ++j) { + double d = 0.0; + for (std::size_t i = 0; i < p_left; ++i) d += u[i] * work(k + i, j); + const double c = tau * d; + for (std::size_t i = 0; i < p_left; ++i) work(k + i, j) -= c * u[i]; + } + + // Accumulate into U_acc from right: U_acc[:, k:] -= tau * (U_acc[:, k:] u) * u^T + for (std::size_t j = 0; j < m; ++j) { + double d = 0.0; + for (std::size_t i = 0; i < p_left; ++i) d += U_acc(j, k + i) * u[i]; + const double c = tau * d; + for (std::size_t i = 0; i < p_left; ++i) U_acc(j, k + i) -= c * u[i]; + } + } + } + + // Right Householder: zero work[k, k+2:]. + if (k + 2 <= n) { + const std::size_t p_right = n - k - 1; // row length after col k+1 + + std::vector<double> v(p_right); + for (std::size_t j = 0; j < p_right; ++j) v[j] = work(k, k + 1 + j); + + double x_norm = 0.0; + for (double val : v) x_norm += val * val; + x_norm = std::sqrt(x_norm); + + if (x_norm > 0.0) { + const double sigma = (v[0] >= 0.0 ? 1.0 : -1.0) * x_norm; + v[0] += sigma; + double vtv = 0.0; + for (double val : v) vtv += val * val; + const double tau = 2.0 / vtv; + + // Apply to work from right: work[:, k+1:] -= tau * (work[:, k+1:] v) * v^T + for (std::size_t i = k; i < m; ++i) { + double d = 0.0; + for (std::size_t j = 0; j < p_right; ++j) d += work(i, k + 1 + j) * v[j]; + const double c = tau * d; + for (std::size_t j = 0; j < p_right; ++j) work(i, k + 1 + j) -= c * v[j]; + } + + // Accumulate into Vt_acc from left: Vt_acc[k+1:, :] -= tau * v * (v^T Vt_acc[k+1:, :]) + for (std::size_t j = 0; j < n; ++j) { + double d = 0.0; + for (std::size_t i = 0; i < p_right; ++i) d += v[i] * Vt_acc(k + 1 + i, j); + const double c = tau * d; + for (std::size_t i = 0; i < p_right; ++i) Vt_acc(k + 1 + i, j) -= c * v[i]; + } + } + } + } + + // Extract bidiagonal: d[i] = diag, e[i] = superdiag. + std::vector<double> d(n), e(n > 1 ? n - 1 : 0, 0.0); + for (std::size_t i = 0; i < n; ++i) d[i] = work(i, i); + for (std::size_t i = 0; i + 1 < n; ++i) e[i] = work(i, i + 1); + + // GR QR sweeps on the bidiagonal. + // Accumulate left rotations into U_f, right rotations into Vt_f. + + Matrix U_f = Matrix::identity(n); + Matrix Vt_f = Matrix::identity(n); + + std::size_t active = n; + int total_iters = 0; + + while (active > 1) { + if (total_iters >= opts.max_iterations) { + std::ostringstream oss; + oss << "svd: did not converge in " << opts.max_iterations << " iterations"; + throw NonConvergenceError(oss.str()); + } + ++total_iters; + + // Deflate small superdiagonals. + while (active > 1) { + const std::size_t i = active - 2; + if (std::abs(e[i]) <= opts.tolerance * (std::abs(d[i]) + std::abs(d[active - 1]))) { + e[i] = 0.0; + --active; + } else { + break; + } + } + if (active <= 1) break; + + // Handle zero on diagonal: if d[k] == 0 for k < active-1, chase the + // nonzero e[k] to zero using a sequence of left Givens rotations. + bool zero_diag = false; + for (std::size_t k = 0; k + 1 < active; ++k) { + if (std::abs(d[k]) <= opts.tolerance) { + zero_diag = true; + // Chase e[k] to zero using left Givens rotations in rows k and k+1..active-1. + double f = e[k]; + e[k] = 0.0; + for (std::size_t j = k + 1; j < active && f != 0.0; ++j) { + const double g = d[j]; + const double r = std::hypot(f, g); + const double c = g / r; + const double s = -f / r; + d[j] = r; + if (j + 1 < active) { + f = s * e[j]; + e[j] *= c; + } + // Accumulate into U_f (left rotation on rows k and j). + // Apply a rotation between rows k and j manually on U_f columns. + for (std::size_t col = 0; col < n; ++col) { + const double uk = U_f(k, col); + const double uj = U_f(j, col); + U_f(k, col) = c * uk - s * uj; + U_f(j, col) = s * uk + c * uj; + } + } + break; + } + } + if (zero_diag) continue; + + // Wilkinson shift from bottom 2×2 of B^T B. + // B^T B bottom-right 2×2 (indices active-2, active-1): + // [d[a-2]^2 + e[a-3]^2, d[a-2]*e[a-2]] + // [d[a-2]*e[a-2], d[a-1]^2 + e[a-2]^2] (if a>=2) + const double a = active >= 2 ? d[active - 2] : 0.0; + const double b_val = active >= 2 ? e[active - 2] : 0.0; + const double c_val = d[active - 1]; + const double t11 = a * a + (active >= 3 ? e[active - 3] * e[active - 3] : 0.0); + const double t12 = a * b_val; + const double t22 = c_val * c_val + b_val * b_val; + const double delta = 0.5 * (t11 - t22); + const double denom = std::abs(delta) + std::hypot(delta, t12); + const double mu = (denom == 0.0) ? t22 + : t22 - (delta >= 0.0 ? 1.0 : -1.0) * (t12 * t12) / denom; + + // Golub-Reinsch implicit QR step. + double f = d[0] * d[0] - mu; + double g = d[0] * e[0]; + + for (std::size_t i = 0; i + 1 < active; ++i) { + // Right Givens: eliminate g from (f, g) in columns i and i+1. + { + const double r = std::hypot(f, g); + const double cr = (r == 0.0) ? 1.0 : f / r; + const double sr = (r == 0.0) ? 0.0 : g / r; + + if (i > 0) e[i - 1] = r; + + f = cr * d[i] + sr * e[i]; + e[i] = -sr * d[i] + cr * e[i]; + g = sr * d[i + 1]; + d[i + 1] *= cr; + + // Accumulate right rotation into Vt_f (acts on cols i and i+1 of V, i.e., rows of Vt_f). + for (std::size_t row = 0; row < n; ++row) { + const double vi = Vt_f(i, row); + const double vi1 = Vt_f(i + 1, row); + Vt_f(i, row) = cr * vi + sr * vi1; + Vt_f(i + 1, row) = -sr * vi + cr * vi1; + } + } + + // Left Givens: eliminate g from (f, g) in rows i and i+1. + { + const double r = std::hypot(f, g); + const double cl = (r == 0.0) ? 1.0 : f / r; + const double sl = (r == 0.0) ? 0.0 : g / r; + + d[i] = r; + + f = cl * e[i] + sl * d[i + 1]; + d[i + 1] = -sl * e[i] + cl * d[i + 1]; + e[i] = f; + + if (i + 2 < active) { + g = sl * e[i + 1]; + e[i + 1] *= cl; + } + + // Accumulate left rotation into U_f (acts on rows i and i+1). + for (std::size_t col = 0; col < n; ++col) { + const double ui = U_f(i, col); + const double ui1 = U_f(i + 1, col); + U_f(i, col) = cl * ui + sl * ui1; + U_f(i + 1, col) = -sl * ui + cl * ui1; + } + } + } + // Set the last diagonal update. + e[active - 2] = f; + } + + // U = U_acc * U_f^T (U_f accumulates row operations on B, + // which correspond to left singular vectors relative to U_acc). + // Vt = Vt_f * Vt_acc (Vt_f accumulates row ops on Vt_acc). + // sigma = |d[i]|, flip signs into U. + + // U_f stores row-wise left rotations applied to the bidiagonal's rows. + // The actual left factor is U_acc * U_f^T (since each left Givens G was + // applied as B <- G B, meaning U_acc absorbs G^T from the right). + const Matrix Uf_t = transpose(U_f); + // U_acc is m×m, Uf_t is n×n. We need m×m U; embed Uf_t into top-left. + Matrix U_full = Matrix::zeros(m, m); + // Copy U_acc * Uf_t into first n columns; remaining m-n columns of U_acc unchanged. + for (std::size_t i = 0; i < m; ++i) { + for (std::size_t j = 0; j < n; ++j) { + double s = 0.0; + for (std::size_t l = 0; l < n; ++l) s += U_acc(i, l) * Uf_t(l, j); + U_full(i, j) = s; + } + for (std::size_t j = n; j < m; ++j) U_full(i, j) = U_acc(i, j); + } + + Matrix Vt_full = Vt_f * Vt_acc; + + Vector sigma(n); + for (std::size_t i = 0; i < n; ++i) { + sigma[i] = std::abs(d[i]); + if (d[i] < 0.0) { + // Negate corresponding column of U (row of U^T) to keep sigma positive. + for (std::size_t j = 0; j < m; ++j) U_full(j, i) = -U_full(j, i); + } + } + + // Sort singular values descending, applying same permutation to U cols and Vt rows. + std::vector<std::size_t> idx(n); + std::iota(idx.begin(), idx.end(), std::size_t{0}); + std::sort(idx.begin(), idx.end(), + [&](std::size_t a, std::size_t b) { return sigma[a] > sigma[b]; }); + + Vector sigma_sorted(n); + Matrix U_sorted(m, m); + Matrix Vt_sorted(n, n); + + for (std::size_t i = 0; i < m; ++i) + for (std::size_t j = n; j < m; ++j) + U_sorted(i, j) = U_full(i, j); + + for (std::size_t rank = 0; rank < n; ++rank) { + const std::size_t src = idx[rank]; + sigma_sorted[rank] = sigma[src]; + for (std::size_t i = 0; i < m; ++i) U_sorted(i, rank) = U_full(i, src); + for (std::size_t j = 0; j < n; ++j) Vt_sorted(rank, j) = Vt_full(src, j); + } + + return SVDResult{std::move(U_sorted), std::move(sigma_sorted), std::move(Vt_sorted)}; +} + +} // namespace linalgebra diff --git a/src/sym_eigen.cpp b/src/sym_eigen.cpp new file mode 100644 index 0000000..a38b272 --- /dev/null +++ b/src/sym_eigen.cpp @@ -0,0 +1,208 @@ +export module linalgebra:sym_eigen; +import std; +import :error; +import :vector; +import :matrix; +import :norms; +import :lu; + +// References used throughout this file: +// T&B — Trefethen & Bau, "Numerical Linear Algebra" +// GVL — Golub & Van Loan, "Matrix Computations" 4th ed. + +export namespace linalgebra { + +// Symmetric tridiagonalization via Householder reflections +// Reference: GVL §8.3.1; T&B Lecture 26 +// +// For symmetric A, computes Q and T such that A = Q T Q^T, +// where T is symmetric tridiagonal and Q is orthogonal. + +struct TridiagonalizeResult { + Matrix T; // symmetric tridiagonal (full n×n) + Matrix Q; // orthogonal: A = Q T Q^T +}; + +[[nodiscard]] TridiagonalizeResult tridiagonalize(const Matrix& A, + double symmetry_tolerance = 1e-12); + +// Eigenvectors via inverse iteration +// Reference: GVL §7.6.1; T&B Lecture 27 +// +// Given A and a vector of approximate eigenvalues, returns a matrix whose +// columns are the corresponding (approximate) eigenvectors. + +struct InverseIterationOptions { + double tolerance = 1e-10; + int max_iterations = 100; +}; + +struct InverseIterationResult { + Matrix eigenvectors; // n × k, column j is eigenvector for eigenvalues[j] + Vector residuals; // ||A*v_j - lambda_j*v_j||_2 for each j +}; + +[[nodiscard]] InverseIterationResult eigenvectors_inverse_iteration( + const Matrix& A, + const Vector& eigenvalues, + InverseIterationOptions opts = {}); + +} // namespace linalgebra + +namespace linalgebra { + +TridiagonalizeResult tridiagonalize(const Matrix& A, double symmetry_tolerance) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "tridiagonalize requires a square matrix, got " + << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + + const std::size_t n = A.rows(); + + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = i + 1; j < n; ++j) + if (std::abs(A(i, j) - A(j, i)) > symmetry_tolerance) + throw LinAlgError("tridiagonalize requires a symmetric matrix"); + + Matrix T = A; + Matrix Q = Matrix::identity(n); + + for (std::size_t k = 0; k + 2 <= n; ++k) { + const std::size_t p = n - k - 1; // length of sub-column + + // Form Householder vector u from T[k+1:, k]. + std::vector<double> u(p); + for (std::size_t i = 0; i < p; ++i) u[i] = T(k + 1 + i, k); + + double x_norm = 0.0; + for (double v : u) x_norm += v * v; + x_norm = std::sqrt(x_norm); + + if (x_norm == 0.0) continue; + + const double sigma = (u[0] >= 0.0 ? 1.0 : -1.0) * x_norm; + u[0] += sigma; + + double utu = 0.0; + for (double v : u) utu += v * v; + const double tau = 2.0 / utu; + + // Apply H = I - tau*u*u^T symmetrically: T <- H T H + // Left: T[k+1:, :] -= tau * u * (u^T T[k+1:, :]) + for (std::size_t j = 0; j < n; ++j) { + double d = 0.0; + for (std::size_t i = 0; i < p; ++i) d += u[i] * T(k + 1 + i, j); + const double coeff = tau * d; + for (std::size_t i = 0; i < p; ++i) T(k + 1 + i, j) -= coeff * u[i]; + } + + // Right: T[:, k+1:] -= tau * (T[:, k+1:] u) * u^T + for (std::size_t j = 0; j < n; ++j) { + double d = 0.0; + for (std::size_t i = 0; i < p; ++i) d += T(j, k + 1 + i) * u[i]; + const double coeff = tau * d; + for (std::size_t i = 0; i < p; ++i) T(j, k + 1 + i) -= coeff * u[i]; + } + + // Explicitly zero sub-subdiagonal entries for numerical cleanliness. + for (std::size_t i = 1; i < p; ++i) { + T(k + 1 + i, k) = 0.0; + T(k, k + 1 + i) = 0.0; + } + + // Accumulate Q: Q[:, k+1:] -= tau * (Q[:, k+1:] u) * u^T + for (std::size_t j = 0; j < n; ++j) { + double d = 0.0; + for (std::size_t i = 0; i < p; ++i) d += Q(j, k + 1 + i) * u[i]; + const double coeff = tau * d; + for (std::size_t i = 0; i < p; ++i) Q(j, k + 1 + i) -= coeff * u[i]; + } + } + + return TridiagonalizeResult{std::move(T), std::move(Q)}; +} + +InverseIterationResult eigenvectors_inverse_iteration(const Matrix& A, + const Vector& eigenvalues, + InverseIterationOptions opts) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "eigenvectors_inverse_iteration requires a square matrix, got " + << A.rows() << "x" << A.cols(); + throw DimensionMismatchError(oss.str()); + } + + const std::size_t n = A.rows(); + const std::size_t k = eigenvalues.size(); + + Matrix evecs(n, k); + Vector residuals(k, 0.0); + + // Starting vector: uniform unit vector. + Vector v0(n, 1.0 / std::sqrt(static_cast<double>(n))); + + for (std::size_t col = 0; col < k; ++col) { + double lambda = eigenvalues[col]; + + // Build shifted matrix B = A - lambda*I. + Matrix B = A; + for (std::size_t i = 0; i < n; ++i) B(i, i) -= lambda; + + // Try to factor; if near-singular (either at factorization or solve time), perturb the shift. + auto make_lu = [&]() -> LUResult { + try { + LUResult result = lu_factor(B, 1e-14); + lu_solve(result, v0); // probe: backward_substitution validates U diagonal + return result; + } catch (const SingularMatrixError&) { + lambda += 1e-7; + for (std::size_t i = 0; i < n; ++i) B(i, i) += 1e-7; + return lu_factor(B, 1e-14); + } + }; + LUResult lu_b = make_lu(); + + Vector v = v0; + + bool converged = false; + for (int iter = 0; iter < opts.max_iterations; ++iter) { + Vector w = lu_solve(lu_b, v); + + // Normalize. + double w_norm = norm2(w); + if (w_norm == 0.0) break; + v = w / w_norm; + + // Residual: ||A v - lambda v||. + const Vector Av = A * v; + double res = 0.0; + for (std::size_t i = 0; i < n; ++i) { + const double d = Av[i] - eigenvalues[col] * v[i]; + res += d * d; + } + res = std::sqrt(res); + + if (res < opts.tolerance) { + converged = true; + residuals[col] = res; + break; + } + } + + if (!converged) { + std::ostringstream oss; + oss << "eigenvectors_inverse_iteration: did not converge for eigenvalue " + << col << " (lambda = " << eigenvalues[col] << ") in " + << opts.max_iterations << " iterations"; + throw NonConvergenceError(oss.str()); + } + + for (std::size_t i = 0; i < n; ++i) evecs(i, col) = v[i]; + } + + return InverseIterationResult{std::move(evecs), std::move(residuals)}; +} + +} // namespace linalgebra 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 <catch2/catch_approx.hpp> +#include <catch2/catch_test_macros.hpp> + +#include <cmath> +#include <cstddef> +#include <numbers> +#include <random> + +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<double>(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<double> 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<double> 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 <catch2/catch_approx.hpp> +#include <catch2/catch_test_macros.hpp> + +#include <cmath> +#include <cstddef> +#include <random> + +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<double> 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<double> 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<double>(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<std::size_t>(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<double> 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<double> 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<double> 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 <catch2/catch_approx.hpp> +#include <catch2/catch_test_macros.hpp> +#include <cstddef> + +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<double>(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 <catch2/catch_approx.hpp> +#include <catch2/catch_test_macros.hpp> + +#include <cmath> +#include <cstddef> +#include <random> + +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<double> 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<double>(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 <catch2/catch_approx.hpp> +#include <catch2/catch_test_macros.hpp> + +#include <cmath> +#include <cstddef> +#include <random> + +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<double> 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); +} |