From 050ccc72c3c763eeffcb866948cdac15f485aa15 Mon Sep 17 00:00:00 2001 From: y-jan137 Date: Sun, 3 May 2026 17:33:09 +0300 Subject: Add cholesky solve --- CMakeLists.txt | 2 + src/cholesky.cpp | 82 +++++++++++++++++++++++++++++ src/linalgebra.cpp | 1 + tests/test_cholesky.cpp | 135 ++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 220 insertions(+) create mode 100644 src/cholesky.cpp create mode 100644 tests/test_cholesky.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 37dcaf7..70f03a1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,6 +36,7 @@ target_sources(linear_algebra src/lu.cpp src/qr.cpp src/qr_iteration.cpp + src/cholesky.cpp ) target_compile_features(linear_algebra PUBLIC cxx_std_23) @@ -90,6 +91,7 @@ if(LINEAR_ALGEBRA_BUILD_TESTS) tests/test_lu.cpp tests/test_qr.cpp tests/test_qr_iteration.cpp + tests/test_cholesky.cpp ) target_link_libraries(linear_algebra_tests diff --git a/src/cholesky.cpp b/src/cholesky.cpp new file mode 100644 index 0000000..f65d9ce --- /dev/null +++ b/src/cholesky.cpp @@ -0,0 +1,82 @@ +export module linalgebra:cholesky; +import std; +import :error; +import :vector; +import :matrix; +import :triangular_solve; + +export namespace linalgebra { + +struct CholeskyResult { + Matrix L; +}; + +CholeskyResult cholesky_factor(const Matrix& A, double tolerance = 1e-12); + +Vector cholesky_solve(const CholeskyResult& chol, const Vector& b); + +} // namespace linalgebra + +namespace linalgebra { + +CholeskyResult cholesky_factor(const Matrix& A, double tolerance) { + if (A.rows() != A.cols()) { + std::ostringstream oss; + oss << "cholesky_factor 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)) > tolerance) { + throw LinAlgError("cholesky_factor requires a symmetric matrix"); + } + } + } + + Matrix L = Matrix::zeros(n, n); + + for (std::size_t j = 0; j < n; ++j) { + double sum = A(j, j); + for (std::size_t k = 0; k < j; ++k) { + sum -= L(j, k) * L(j, k); + } + + if (sum <= tolerance) { + std::ostringstream oss; + oss << "cholesky_factor: matrix is not positive definite (diagonal became " + << sum << " at step " << j << ")"; + throw LinAlgError(oss.str()); + } + + L(j, j) = std::sqrt(sum); + + for (std::size_t i = j + 1; i < n; ++i) { + double s = A(i, j); + for (std::size_t k = 0; k < j; ++k) { + s -= L(i, k) * L(j, k); + } + L(i, j) = s / L(j, j); + } + } + + return CholeskyResult{std::move(L)}; +} + +Vector cholesky_solve(const CholeskyResult& chol, const Vector& b) { + const std::size_t n = chol.L.rows(); + + if (b.size() != n) { + std::ostringstream oss; + oss << "cholesky_solve: rhs size " << b.size() + << " does not match factorization size " << n; + throw DimensionMismatchError(oss.str()); + } + + const Vector y = forward_substitution(chol.L, b); + return backward_substitution(transpose(chol.L), y); +} + +} // namespace linalgebra diff --git a/src/linalgebra.cpp b/src/linalgebra.cpp index 302810f..18e9cf3 100644 --- a/src/linalgebra.cpp +++ b/src/linalgebra.cpp @@ -8,3 +8,4 @@ export import :triangular_solve; export import :lu; export import :qr; export import :qr_iteration; +export import :cholesky; diff --git a/tests/test_cholesky.cpp b/tests/test_cholesky.cpp new file mode 100644 index 0000000..a92563f --- /dev/null +++ b/tests/test_cholesky.cpp @@ -0,0 +1,135 @@ +import linalgebra; + +#include +#include + +#include +#include +#include + +using linalgebra::CholeskyResult; +using linalgebra::DimensionMismatchError; +using linalgebra::LinAlgError; +using linalgebra::Matrix; +using linalgebra::Vector; + +namespace { + +double reconstruction_error(const Matrix& A, const CholeskyResult& chol) { + const Matrix LLT = chol.L * linalgebra::transpose(chol.L); + const std::size_t n = A.rows(); + 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 = A(i, j) - LLT(i, j); + err += d * d; + } + } + return std::sqrt(err); +} + +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 + +TEST_CASE("Cholesky factor 2x2", "[cholesky]") { + Matrix A{{4.0, 2.0}, {2.0, 3.0}}; + auto chol = linalgebra::cholesky_factor(A); + + REQUIRE(reconstruction_error(A, chol) < 1e-12); + + REQUIRE(chol.L(0, 0) == Catch::Approx(2.0)); + REQUIRE(chol.L(1, 0) == Catch::Approx(1.0)); + REQUIRE(chol.L(1, 1) == Catch::Approx(std::sqrt(2.0))); + REQUIRE(chol.L(0, 1) == Catch::Approx(0.0).margin(1e-15)); +} + +TEST_CASE("Cholesky factor 3x3", "[cholesky]") { + Matrix A{{25.0, 15.0, -5.0}, + {15.0, 18.0, 0.0}, + {-5.0, 0.0, 11.0}}; + auto chol = linalgebra::cholesky_factor(A); + REQUIRE(reconstruction_error(A, chol) < 1e-12); +} + +TEST_CASE("Cholesky factor identity", "[cholesky]") { + auto I = Matrix::identity(5); + auto chol = linalgebra::cholesky_factor(I); + REQUIRE(reconstruction_error(I, chol) < 1e-14); + for (std::size_t i = 0; i < 5; ++i) { + REQUIRE(chol.L(i, i) == Catch::Approx(1.0)); + } +} + +TEST_CASE("Cholesky factor random SPD", "[cholesky]") { + std::mt19937 rng(42); + for (std::size_t n : {4, 8, 16, 32}) { + Matrix A = make_spd(n, rng); + auto chol = linalgebra::cholesky_factor(A); + REQUIRE(reconstruction_error(A, chol) < 1e-10); + } +} + +TEST_CASE("Cholesky solve", "[cholesky]") { + Matrix A{{4.0, 2.0}, {2.0, 3.0}}; + Vector b{8.0, 7.0}; + auto chol = linalgebra::cholesky_factor(A); + Vector x = linalgebra::cholesky_solve(chol, b); + + double residual = linalgebra::norm2(A * x - b); + REQUIRE(residual < 1e-12); +} + +TEST_CASE("Cholesky solve random SPD", "[cholesky]") { + std::mt19937 rng(123); + std::uniform_real_distribution dist(-5.0, 5.0); + + for (std::size_t n : {5, 10, 20}) { + Matrix A = make_spd(n, rng); + Vector b(n); + for (std::size_t i = 0; i < n; ++i) { + b[i] = dist(rng); + } + + auto chol = linalgebra::cholesky_factor(A); + Vector x = linalgebra::cholesky_solve(chol, b); + + double residual = linalgebra::norm2(A * x - b); + REQUIRE(residual < 1e-10); + } +} + +TEST_CASE("Cholesky rejects non-square matrix", "[cholesky]") { + Matrix A(3, 4); + REQUIRE_THROWS_AS(linalgebra::cholesky_factor(A), DimensionMismatchError); +} + +TEST_CASE("Cholesky rejects non-symmetric matrix", "[cholesky]") { + Matrix A{{1.0, 2.0}, {3.0, 4.0}}; + REQUIRE_THROWS_AS(linalgebra::cholesky_factor(A), LinAlgError); +} + +TEST_CASE("Cholesky rejects non-positive-definite matrix", "[cholesky]") { + Matrix A{{1.0, 0.0}, {0.0, -1.0}}; + REQUIRE_THROWS_AS(linalgebra::cholesky_factor(A), LinAlgError); +} + +TEST_CASE("Cholesky solve dimension mismatch", "[cholesky]") { + Matrix A{{4.0, 2.0}, {2.0, 3.0}}; + auto chol = linalgebra::cholesky_factor(A); + Vector b{1.0, 2.0, 3.0}; + REQUIRE_THROWS_AS(linalgebra::cholesky_solve(chol, b), DimensionMismatchError); +} -- cgit v1.2.3