Add LU and experiment

5 files changed, 710 insertions, 0 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 54f9d5b..87cd6a3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -17,6 +17,7 @@ add_library(linear_algebra
     src/matrix.cpp
     src/norms.cpp
     src/triangular_solve.cpp
+    src/lu.cpp
 )
 
 add_library(linear_algebra::core ALIAS linear_algebra)
@@ -60,6 +61,12 @@ if(LINEAR_ALGEBRA_BUILD_EXAMPLES)
     target_link_libraries(matmul PRIVATE linear_algebra::core)
 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)
+endif()
+
 if(LINEAR_ALGEBRA_BUILD_TESTS)
     include(FetchContent)
 endif()
@@ -80,6 +87,7 @@ if(LINEAR_ALGEBRA_BUILD_TESTS)
         tests/test_vector.cpp
         tests/test_matrix.cpp
         tests/test_triangular_solve.cpp
+        tests/test_lu.cpp
     )
 
     target_link_libraries(linear_algebra_tests
diff --git a/experiments/pivoting_vs_no_pivoting.cpp b/experiments/pivoting_vs_no_pivoting.cpp
new file mode 100644
index 0000000..ce2b9cc
--- /dev/null
+++ b/experiments/pivoting_vs_no_pivoting.cpp
@@ -0,0 +1,260 @@
+// Experiment: partial pivoting vs no-pivot LU
+//
+// Demonstrates why partial pivoting is essential for numerical stability.
+// Run the binary and inspect the residuals printed to stdout.
+
+#include "lu.hpp"
+#include "matrix.hpp"
+#include "norms.hpp"
+#include "triangular_solve.hpp"
+#include "vector.hpp"
+
+#include <cmath>
+#include <cstddef>
+#include <iomanip>
+#include <iostream>
+#include <optional>
+#include <random>
+#include <string>
+#include <vector>
+
+using linalg::Matrix;
+using linalg::Vector;
+
+// ---------------------------------------------------------------------------
+// Local no-pivot LU for comparison only.
+// This is intentionally naive — it is here to show what breaks without pivoting.
+// ---------------------------------------------------------------------------
+
+struct NoPivotLU {
+    Matrix L;
+    Matrix U;
+    bool failed = false;   // true if a zero pivot was encountered
+    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};
+}
+
+// Solve using a no-pivot LU (L unit lower triangular, U upper triangular).
+// If the factorization failed or U is numerically singular, returns nullopt.
+std::optional<Vector> solve_no_pivot(const NoPivotLU& f, const Vector& b) {
+    if (f.failed) return std::nullopt;
+    try {
+        const Vector y = linalg::forward_substitution(f.L, b, 1e-14, /*unit_diagonal=*/true);
+        return linalg::backward_substitution(f.U, y);
+    } catch (...) {
+        return std::nullopt;
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Metrics
+// ---------------------------------------------------------------------------
+
+double solve_residual(const Matrix& A, const Vector& x, const Vector& b) {
+    return linalg::norm2(A * x - b);
+}
+
+double reconstruction_error(const Matrix& A, const linalg::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);
+}
+
+// ---------------------------------------------------------------------------
+// Reporting
+// ---------------------------------------------------------------------------
+
+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 linalg::LUResult lu = linalg::lu_factor(A);
+        const Vector x = linalg::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;
+    }
+    // Compute reconstruction error without perm (no-pivot uses A directly).
+    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);
+}
+
+// ---------------------------------------------------------------------------
+// Experiment cases
+// ---------------------------------------------------------------------------
+
+// 1. Random well-conditioned matrix
+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);
+}
+
+// 2. Badly row-scaled matrix
+//    Rows differ in magnitude by ~10^14.  Without pivoting, tiny early pivots
+//    amplify round-off; with pivoting, the large-row is selected first.
+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};  // b = A * [1, 1, 1]
+    run_case("Badly scaled (row norms differ by 10^14)", A, b);
+}
+
+// 3. Classic pathological example for no-pivot LU.
+//    With epsilon = 1e-15, no-pivot computes a huge multiplier (1/epsilon),
+//    which causes catastrophic cancellation in the updated rows.
+//    With pivoting, we swap first and the multiplier is bounded by 1.
+void exp_epsilon_pathology() {
+    constexpr double eps = 1e-15;
+    const Matrix A{{eps, 1.0}, {1.0, 2.0}};
+    //  True solution of [eps 1; 1 2] * x = [1+eps; 3] is x = [1; 1].
+    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";
+}
+
+// 4. Matrix where no-pivot LU diverges visibly on a 4x4 example.
+//    The first pivot is small (0.001) but rows below have entries ~1000.
+//    No pivot causes multipliers of magnitude 10^6, annihilating subdiagonal info.
+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";
+}
+
+// 5. Matrix requiring multiple row swaps (permutation is non-trivial).
+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);
+}
+
+// ---------------------------------------------------------------------------
+// main
+// ---------------------------------------------------------------------------
+
+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();
+
+    std::cout << "\n" << std::string(60, '=') << "\n";
+    std::cout << "  Conclusion\n";
+    std::cout << std::string(60, '=') << "\n";
+    std::cout <<
+        "Partial pivoting keeps multipliers bounded by 1 in magnitude.\n"
+        "Without it, a near-zero pivot inflates multipliers and destroys\n"
+        "accuracy via catastrophic cancellation.  The 'epsilon pathology'\n"
+        "case is the textbook example: a pivot of 1e-15 makes no-pivot LU\n"
+        "compute x ≈ [0, 0.5] instead of the exact [1, 1].\n\n";
+
+    return 0;
+}
diff --git a/include/lu.hpp b/include/lu.hpp
new file mode 100644
index 0000000..6717733
--- /dev/null
+++ b/include/lu.hpp
@@ -0,0 +1,43 @@
+#pragma once
+
+#include <cstddef>
+#include <vector>
+
+#include "matrix.hpp"
+#include "vector.hpp"
+
+namespace linalg {
+
+// Result of LU factorization with partial pivoting.
+//
+// The factorization satisfies PA = LU, where:
+//   P  is the permutation matrix encoded by `perm`
+//   L  is unit lower triangular  (L[i][i] == 1)
+//   U  is upper triangular
+//
+// `perm[i]` = index of the original row that ended up at position i.
+// Applying P to a vector b means: (Pb)[i] = b[perm[i]].
+//
+// `sign` is the sign of the permutation: +1 if an even number of row
+// swaps were made, -1 if odd.  Useful for computing det(A) = sign * prod(diag(U)).
+struct LUResult {
+    Matrix L;
+    Matrix U;
+    std::vector<std::size_t> perm;
+    int sign;
+};
+
+// Compute the LU factorization of A with partial pivoting.
+//
+// Throws DimensionMismatchError if A is not square.
+// Throws SingularMatrixError   if A is (numerically) singular, i.e. any
+//   pivot is smaller in magnitude than `singular_tolerance`.
+LUResult lu_factor(const Matrix& A, double singular_tolerance = 1e-12);
+
+// Solve Ax = b given a precomputed LU factorization.
+//
+// Applies the stored permutation, then forward / backward substitution.
+// Throws DimensionMismatchError if b.size() != lu.L.rows().
+Vector lu_solve(const LUResult& lu, const Vector& b);
+
+}  // namespace linalg
diff --git a/src/lu.cpp b/src/lu.cpp
new file mode 100644
index 0000000..9b2c950
--- /dev/null
+++ b/src/lu.cpp
@@ -0,0 +1,112 @@
+#include "lu.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+#include <sstream>
+
+#include "linalg_error.hpp"
+#include "triangular_solve.hpp"
+
+namespace linalg {
+
+LUResult lu_factor(const Matrix& A, double singular_tolerance) {
+    if (A.rows() != A.cols()) {
+        std::ostringstream oss;
+        oss << "lu_factor requires a square matrix, got " << A.rows() << "x" << A.cols();
+        throw DimensionMismatchError(oss.str());
+    }
+
+    const std::size_t n = A.rows();
+
+    // Working copy: elimination is performed in-place here.
+    Matrix work = A;
+
+    // L starts as identity; multipliers fill the strict lower triangle.
+    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);
+
+    // perm[i] = original row index now at position i.
+    std::vector<std::size_t> perm(n);
+    std::iota(perm.begin(), perm.end(), std::size_t{0});
+    int sign = 1;
+
+    for (std::size_t k = 0; k < n; ++k) {
+        // ---- Partial pivoting: find row with largest magnitude in column k ----
+        std::size_t pivot_row = k;
+        double max_val = std::abs(work(k, k));
+        for (std::size_t i = k + 1; i < n; ++i) {
+            const double val = std::abs(work(i, k));
+            if (val > max_val) {
+                max_val = val;
+                pivot_row = i;
+            }
+        }
+
+        if (pivot_row != k) {
+            // Swap rows in the working matrix.
+            for (std::size_t j = 0; j < n; ++j) {
+                std::swap(work(k, j), work(pivot_row, j));
+            }
+            // Swap already-computed multipliers in L (columns 0 .. k-1).
+            for (std::size_t j = 0; j < k; ++j) {
+                std::swap(L(k, j), L(pivot_row, j));
+            }
+            std::swap(perm[k], perm[pivot_row]);
+            sign = -sign;
+        }
+
+        // ---- Singularity check ----
+        if (std::abs(work(k, k)) <= singular_tolerance) {
+            std::ostringstream oss;
+            oss << "lu_factor: near-zero pivot " << work(k, k) << " at step " << k
+                << " (tolerance " << singular_tolerance << ")";
+            throw SingularMatrixError(oss.str());
+        }
+
+        // ---- Record U row k ----
+        for (std::size_t j = k; j < n; ++j) {
+            U(k, j) = work(k, j);
+        }
+
+        // ---- Compute multipliers and eliminate below pivot ----
+        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);
+            }
+            work(i, k) = 0.0;
+        }
+    }
+
+    return LUResult{std::move(L), std::move(U), std::move(perm), sign};
+}
+
+Vector lu_solve(const LUResult& lu, const Vector& b) {
+    const std::size_t n = lu.L.rows();
+
+    if (b.size() != n) {
+        std::ostringstream oss;
+        oss << "lu_solve: rhs size " << b.size() << " does not match factorization size " << n;
+        throw DimensionMismatchError(oss.str());
+    }
+
+    // Step 1: apply permutation P.  (Pb)[i] = b[perm[i]]
+    Vector pb(n);
+    for (std::size_t i = 0; i < n; ++i) {
+        pb[i] = b[lu.perm[i]];
+    }
+
+    // Step 2: forward substitution  Ly = Pb  (L has unit diagonal)
+    const Vector y = forward_substitution(lu.L, pb, /*singular_tolerance=*/1e-14,
+                                          /*unit_diagonal=*/true);
+
+    // Step 3: backward substitution  Ux = y
+    return backward_substitution(lu.U, y);
+}
+
+}  // namespace linalg
diff --git a/tests/test_lu.cpp b/tests/test_lu.cpp
new file mode 100644
index 0000000..e54a8a1
--- /dev/null
+++ b/tests/test_lu.cpp
@@ -0,0 +1,287 @@
+#include "linalg_error.hpp"
+#include "lu.hpp"
+#include "matrix.hpp"
+#include "norms.hpp"
+#include "vector.hpp"
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <random>
+
+using linalg::DimensionMismatchError;
+using linalg::Matrix;
+using linalg::SingularMatrixError;
+using linalg::Vector;
+using linalg::LUResult;
+
+// ---------------------------------------------------------------------------
+// Helpers
+// ---------------------------------------------------------------------------
+
+namespace {
+
+// ||PA - LU||_F  (Frobenius, computed element-wise via ||vec||_2)
+double reconstruction_error(const Matrix& A, const LUResult& lu) {
+    const std::size_t n = A.rows();
+    // Build PA by permuting rows of A.
+    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);
+        }
+    }
+    // Compute LU product.
+    const Matrix LU = lu.L * lu.U;
+    // Compute Frobenius norm of (PA - LU).
+    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(i, j);
+            err += d * d;
+        }
+    }
+    return std::sqrt(err);
+}
+
+// ||Ax - b||_2
+double solve_residual(const Matrix& A, const Vector& x, const Vector& b) {
+    return linalg::norm2(A * x - b);
+}
+
+// Generate a reproducible random nonsingular n x n matrix.
+Matrix random_matrix(std::size_t n, unsigned seed = 42) {
+    std::mt19937 rng(seed);
+    std::uniform_real_distribution<double> dist(-10.0, 10.0);
+    Matrix M(n, n);
+    for (std::size_t i = 0; i < n; ++i) {
+        for (std::size_t j = 0; j < n; ++j) {
+            M(i, j) = dist(rng);
+        }
+    }
+    return M;
+}
+
+}  // namespace
+
+// ---------------------------------------------------------------------------
+// Factorization correctness
+// ---------------------------------------------------------------------------
+
+TEST_CASE("LU factorization: 3x3 known system", "[lu]") {
+    const Matrix A{
+        {2.0, 1.0, -1.0},
+        {-3.0, -1.0, 2.0},
+        {-2.0, 1.0, 2.0}
+    };
+    const LUResult lu = linalg::lu_factor(A);
+
+    REQUIRE(lu.L.rows() == 3);
+    REQUIRE(lu.U.rows() == 3);
+    REQUIRE(lu.perm.size() == 3);
+
+    // L must have unit diagonal.
+    for (std::size_t i = 0; i < 3; ++i) {
+        CHECK(lu.L(i, i) == Catch::Approx(1.0));
+    }
+
+    // Reconstruction: ||PA - LU|| must be near zero.
+    CHECK(reconstruction_error(A, lu) == Catch::Approx(0.0).margin(1e-12));
+}
+
+TEST_CASE("LU factorization: identity matrix", "[lu]") {
+    const Matrix I = Matrix::identity(4);
+    const LUResult lu = linalg::lu_factor(I);
+    CHECK(reconstruction_error(I, lu) == Catch::Approx(0.0).margin(1e-14));
+    // U should equal I (up to row ordering already handled by PA=LU).
+    for (std::size_t i = 0; i < 4; ++i) {
+        CHECK(lu.U(i, i) == Catch::Approx(1.0));
+    }
+}
+
+TEST_CASE("LU factorization: matrix requiring row swaps", "[lu]") {
+    // First column entry is zero — no-pivot LU would immediately fail.
+    const Matrix A{
+        {0.0, 1.0, 2.0},
+        {3.0, 4.0, 5.0},
+        {6.0, 7.0, 8.0}
+    };
+    // A is singular (rows are in AP), but check that partial pivoting still
+    // proceeds and detects singularity correctly.
+    // Row 3 - row 2 = row 2 - row 1, so rank < 3.
+    CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError);
+}
+
+TEST_CASE("LU factorization: first-column zero, nonsingular", "[lu]") {
+    // [[0, 1], [1, 0]] — requires a swap at step 0.
+    const Matrix A{{0.0, 1.0}, {1.0, 0.0}};
+    const LUResult lu = linalg::lu_factor(A);
+    CHECK(reconstruction_error(A, lu) == Catch::Approx(0.0).margin(1e-14));
+    // Solving Ax = b: A swaps components.
+    const Vector b{3.0, 7.0};
+    const Vector x = linalg::lu_solve(lu, b);
+    CHECK(solve_residual(A, x, b) == Catch::Approx(0.0).margin(1e-12));
+    CHECK(x[0] == Catch::Approx(7.0));
+    CHECK(x[1] == Catch::Approx(3.0));
+}
+
+TEST_CASE("LU factorization: random nonsingular matrices", "[lu]") {
+    for (std::size_t n : {5u, 10u, 20u}) {
+        const Matrix A = random_matrix(n, 123u + static_cast<unsigned>(n));
+        const LUResult lu = linalg::lu_factor(A);
+        CHECK(reconstruction_error(A, lu) == Catch::Approx(0.0).margin(1e-10));
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Solve correctness
+// ---------------------------------------------------------------------------
+
+TEST_CASE("LU solve: known 3x3 system", "[lu]") {
+    // From Cramer / textbook: solution is x = (2, 3, -1).
+    const Matrix A{
+        {2.0, 1.0, -1.0},
+        {-3.0, -1.0, 2.0},
+        {-2.0, 1.0, 2.0}
+    };
+    const Vector expected{2.0, 3.0, -1.0};
+    const Vector b = A * expected;
+
+    const LUResult lu = linalg::lu_factor(A);
+    const Vector x = linalg::lu_solve(lu, b);
+
+    CHECK(x[0] == Catch::Approx(expected[0]).epsilon(1e-12));
+    CHECK(x[1] == Catch::Approx(expected[1]).epsilon(1e-12));
+    CHECK(x[2] == Catch::Approx(expected[2]).epsilon(1e-12));
+    CHECK(solve_residual(A, x, b) == Catch::Approx(0.0).margin(1e-12));
+}
+
+TEST_CASE("LU solve: random nonsingular systems", "[lu]") {
+    for (std::size_t n : {5u, 15u, 30u}) {
+        const Matrix A = random_matrix(n, 7u * static_cast<unsigned>(n));
+        const LUResult lu = linalg::lu_factor(A);
+
+        // Random rhs.
+        std::mt19937 rng(n);
+        std::uniform_real_distribution<double> dist(-5.0, 5.0);
+        Vector b(n);
+        for (std::size_t i = 0; i < n; ++i) b[i] = dist(rng);
+
+        const Vector x = linalg::lu_solve(lu, b);
+        CHECK(solve_residual(A, x, b) == Catch::Approx(0.0).margin(1e-9));
+    }
+}
+
+TEST_CASE("LU solve: diagonal system", "[lu]") {
+    // D = diag(2, 3, 4), b = (2, 9, 8), solution = (1, 3, 2).
+    const Matrix D{
+        {2.0, 0.0, 0.0},
+        {0.0, 3.0, 0.0},
+        {0.0, 0.0, 4.0}
+    };
+    const Vector b{2.0, 9.0, 8.0};
+    const LUResult lu = linalg::lu_factor(D);
+    const Vector x = linalg::lu_solve(lu, b);
+
+    CHECK(x[0] == Catch::Approx(1.0));
+    CHECK(x[1] == Catch::Approx(3.0));
+    CHECK(x[2] == Catch::Approx(2.0));
+}
+
+// ---------------------------------------------------------------------------
+// Failure cases
+// ---------------------------------------------------------------------------
+
+TEST_CASE("LU factorization: non-square matrix throws", "[lu]") {
+    const Matrix A(3, 4);
+    CHECK_THROWS_AS(linalg::lu_factor(A), DimensionMismatchError);
+}
+
+TEST_CASE("LU factorization: exactly singular matrix throws", "[lu]") {
+    // Zero row → singular.
+    const Matrix A{
+        {1.0, 2.0, 3.0},
+        {4.0, 5.0, 6.0},
+        {0.0, 0.0, 0.0}
+    };
+    CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError);
+}
+
+TEST_CASE("LU factorization: rank-deficient matrix throws", "[lu]") {
+    // Row 2 is a linear combination of rows 0 and 1.
+    const Matrix A{
+        {1.0, 2.0},
+        {2.0, 4.0}
+    };
+    CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError);
+}
+
+TEST_CASE("LU factorization: near-singular matrix throws at default tolerance", "[lu]") {
+    // Pivot reduced to ~1e-16, should trip the singularity check.
+    const Matrix A{
+        {1.0, 1.0},
+        {1.0, 1.0 + 1e-16}
+    };
+    CHECK_THROWS_AS(linalg::lu_factor(A), SingularMatrixError);
+}
+
+TEST_CASE("LU solve: mismatched rhs throws", "[lu]") {
+    const Matrix A = Matrix::identity(3);
+    const LUResult lu = linalg::lu_factor(A);
+    const Vector b(5, 1.0);
+    CHECK_THROWS_AS(linalg::lu_solve(lu, b), DimensionMismatchError);
+}
+
+// ---------------------------------------------------------------------------
+// L and U structure
+// ---------------------------------------------------------------------------
+
+TEST_CASE("LU factorization: L is unit lower triangular", "[lu]") {
+    const Matrix A = random_matrix(6, 999u);
+    const LUResult lu = linalg::lu_factor(A);
+    const std::size_t n = A.rows();
+
+    for (std::size_t i = 0; i < n; ++i) {
+        // Unit diagonal.
+        CHECK(lu.L(i, i) == Catch::Approx(1.0));
+        // Strict upper triangle is zero.
+        for (std::size_t j = i + 1; j < n; ++j) {
+            CHECK(lu.L(i, j) == Catch::Approx(0.0).margin(1e-15));
+        }
+    }
+}
+
+TEST_CASE("LU factorization: U is upper triangular", "[lu]") {
+    const Matrix A = random_matrix(6, 777u);
+    const LUResult lu = linalg::lu_factor(A);
+    const std::size_t n = A.rows();
+
+    for (std::size_t i = 1; i < n; ++i) {
+        for (std::size_t j = 0; j < i; ++j) {
+            CHECK(lu.U(i, j) == Catch::Approx(0.0).margin(1e-15));
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Permutation sign and determinant
+// ---------------------------------------------------------------------------
+
+TEST_CASE("LU factorization: sign of permutation is ±1", "[lu]") {
+    const Matrix A = random_matrix(5, 321u);
+    const LUResult lu = linalg::lu_factor(A);
+    CHECK((lu.sign == 1 || lu.sign == -1));
+}
+
+TEST_CASE("LU factorization: determinant via sign * prod(diag(U))", "[lu]") {
+    // det([[3,1],[2,4]]) = 12 - 2 = 10
+    const Matrix A{{3.0, 1.0}, {2.0, 4.0}};
+    const LUResult lu = linalg::lu_factor(A);
+    double det = static_cast<double>(lu.sign);
+    for (std::size_t i = 0; i < A.rows(); ++i) {
+        det *= lu.U(i, i);
+    }
+    CHECK(det == Catch::Approx(10.0).epsilon(1e-12));
+}