Add cholesky solve

4 files changed, 220 insertions, 0 deletions

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 <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <random>
+
+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<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
+
+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<double> 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);
+}