From 69a18d17b17d1e1d2a17dfd087d2f20dcba05715 Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Thu, 30 Jul 2026 10:49:39 +0000 Subject: [PATCH 1/3] add lu decomposition --- README.md | 2 +- ROADMAP.md | 1 - doc/solvers/LuDecomposition.md | 151 ++++++++++++++++ doc/solvers/README.md | 1 + numerical/math/TriangularSolve.hpp | 17 ++ numerical/math/test/TestTriangularSolve.cpp | 23 +++ numerical/solvers/CMakeLists.txt | 2 + numerical/solvers/LuDecomposition.cpp | 7 + numerical/solvers/LuDecomposition.hpp | 169 ++++++++++++++++++ numerical/solvers/test/CMakeLists.txt | 1 + .../solvers/test/TestLuDecomposition.cpp | 159 ++++++++++++++++ .../solvers/LuDecomposition/explanation.md | 27 --- .../solvers/LuDecomposition/implementation.md | 78 -------- roadmap/solvers/LuDecomposition/tests.md | 58 ------ 14 files changed, 531 insertions(+), 165 deletions(-) create mode 100644 doc/solvers/LuDecomposition.md create mode 100644 numerical/solvers/LuDecomposition.cpp create mode 100644 numerical/solvers/LuDecomposition.hpp create mode 100644 numerical/solvers/test/TestLuDecomposition.cpp delete mode 100644 roadmap/solvers/LuDecomposition/explanation.md delete mode 100644 roadmap/solvers/LuDecomposition/implementation.md delete mode 100644 roadmap/solvers/LuDecomposition/tests.md diff --git a/README.md b/README.md index 44258e5e..2ef4c48a 100644 --- a/README.md +++ b/README.md @@ -27,7 +27,7 @@ Refer to the documentation to quickly integrate and utilize the library's signal | [Optimization](doc/optimization/README.md) | Gradient Descent | | [Regularization](doc/regularization/README.md) | L1 (Lasso), L2 (Ridge) | | [Math](doc/math/README.md) | CORDIC, Quaternion, MatrixNorms, Step Response Metrics | -| [Solvers](doc/solvers/README.md) | Gaussian Elimination, Levinson-Durbin, Durand-Kerner, Cholesky, DARE, Runge-Kutta ODE Integrators (RK4 + Dormand-Prince), Spectral Radius & Discrete Stability Margin, QR Decomposition (Householder / Givens) | +| [Solvers](doc/solvers/README.md) | Gaussian Elimination, Levinson-Durbin, Durand-Kerner, Cholesky, DARE, Runge-Kutta ODE Integrators (RK4 + Dormand-Prince), Spectral Radius & Discrete Stability Margin, QR Decomposition (Householder / Givens), LU Decomposition with Partial Pivoting | | [Performance Optimization](doc/performance-optimization/README.md) | Compiler optimizations, SIMD | Each category page lists its algorithms with a brief description and links to the detailed documentation. diff --git a/ROADMAP.md b/ROADMAP.md index c9473cb9..8e1f585b 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -27,7 +27,6 @@ Difficulty legend: | # | Component | Target module | Difficulty | |----|------------------------------------------------------|---------------------------|------------| -| 28 | LU decomposition with partial pivoting | `solvers` | ★★★★☆ | | 29 | Matrix exponential (scaling & squaring + Padé) | `math` | ★★★★☆ | | 30 | Continuous → discrete conversion (`c2d`) | `math` | ★★★★☆ | | 31 | Lyapunov / Sylvester equation solvers | `solvers` | ★★★★☆ | diff --git a/doc/solvers/LuDecomposition.md b/doc/solvers/LuDecomposition.md new file mode 100644 index 00000000..7365136a --- /dev/null +++ b/doc/solvers/LuDecomposition.md @@ -0,0 +1,151 @@ +# LU Decomposition with Partial Pivoting + +## Overview & Motivation + +Many embedded control and estimation algorithms reduce to solving a dense linear system +$A x = b$ for a general (non-symmetric) matrix $A$. LU decomposition with partial pivoting +factors $A$ as $P A = L U$, where $P$ is a permutation, $L$ is unit lower-triangular, and +$U$ is upper-triangular. Factoring $A$ once and reusing $L$ and $U$ turns every subsequent +solve into two cheap triangular sweeps — a critical advantage in loops where the same plant +or observer matrix must be inverted against many right-hand sides at different time steps. + +## Mathematical Theory + +### Definitions + +For an $n \times n$ matrix $A$, partial-pivoting LU factorization produces + +$$P A = L U,$$ + +where: +- $P \in \{0,1\}^{n \times n}$ is a permutation matrix encoding row swaps, +- $L \in \mathbb{R}^{n \times n}$ is unit lower-triangular ($L_{ii} = 1$, $L_{ij} = 0$ for $j > i$), +- $U \in \mathbb{R}^{n \times n}$ is upper-triangular ($U_{ij} = 0$ for $i > j$). + +### Gaussian Elimination with Partial Pivoting (GEPP) + +At step $k$ of the outer loop ($k = 0, \ldots, n-1$): + +1. **Pivot selection**: find $p = \arg\max_{i \geq k} |A_{ik}|$ (the largest-magnitude entry + in the current column from the diagonal down). +2. **Row swap**: interchange rows $k$ and $p$ in the working matrix and record the swap in $P$. +3. **Elimination**: for each $i > k$ compute the multiplier $\ell_{ik} = A_{ik} / A_{kk}$, + store it below the diagonal, and update the trailing submatrix: + +$$A_{ij} \leftarrow A_{ij} - \ell_{ik} A_{kj}, \quad j > k.$$ + +After $n - 1$ steps the working matrix holds $U$ on and above the diagonal and the multipliers +$\ell_{ik}$ below it, packed together in a single $n \times n$ array. + +### Triangular Solves + +Given $P A = L U$ and a right-hand side $b$, the system $A x = b$ becomes $L U x = P b$. +Solving in two passes: + +$$L y = P b \quad \text{(forward substitution)},$$ +$$U x = y \quad \text{(back substitution)}.$$ + +Both passes cost $O(n^2)$ operations, and both reuse the shared +`math::SolveUnitLowerTriangular` / `math::SolveUpperTriangular` kernels +(from [`math/TriangularSolve.hpp`](../../numerical/math/TriangularSolve.hpp)) +applied directly to the packed factor — the same triangular-solve building blocks +that back `QrDecomposition`. + +The factored form is exposed for reuse by other components: `L()`, `U()`, and `P()` +return the individual factors (so callers can verify $P A = L U$ or build derived +quantities), and `IsSingular()` reports whether the last `Decompose` aborted. + +### Determinant + +The determinant of $A$ equals the product of $U$'s diagonal entries multiplied by the +sign of $P$: + +$$\det(A) = \operatorname{sign}(P) \prod_{k=0}^{n-1} U_{kk}.$$ + +The sign flips each time rows are swapped; it starts at $+1$ and is multiplied by $-1$ per swap. + +### Inverse + +The inverse is formed by solving $A X = I$ column by column — $n$ back-and-forward +substitution pairs — at a total cost of $O(n^3)$. + +## Complexity Analysis + +| Operation | Time | Space | Notes | +|--------------|-----------------------------|----------|-----------------------------------------------| +| `Decompose` | $O\!\left(\tfrac{2}{3}n^3\right)$ | $O(n^2)$ | GEPP in-place; multipliers packed with $U$ | +| `Solve` | $O(n^2)$ | $O(n)$ | Two triangular sweeps; reuses stored factors | +| `Determinant`| $O(n)$ | $O(1)$ | One pass over $U$'s diagonal | +| `Inverse` | $O(n^3)$ | $O(n^2)$ | $n$ column solves; avoid when `Solve` suffices | + +All storage is stack-allocated; the packed LU matrix and pivot vector are bounded by $N$. + +## Step-by-Step Walkthrough + +Factor $A = \begin{pmatrix} 2 & 1 & -1 \\ -3 & -1 & 2 \\ -2 & 1 & 2 \end{pmatrix}$. + +**Step $k = 0$:** Largest magnitude in column 0 is $|-3| = 3$ at row 1. Swap rows 0 and 1. +Working matrix: $\begin{pmatrix} -3 & -1 & 2 \\ 2 & 1 & -1 \\ -2 & 1 & 2 \end{pmatrix}$. +Multipliers: $\ell_{10} = 2/(-3) = -2/3$, $\ell_{20} = -2/(-3) = 2/3$. +After elimination: $U_{0\cdot} = (-3,-1,2)$; rows 1 and 2 updated. + +**Step $k = 1$:** Select pivot in column 1 from rows 1 onward. Possibly swap. Compute +multiplier $\ell_{21}$ and update $U_{2\cdot}$. + +**Step $k = 2$:** Single diagonal entry remains; $U_{22}$ is read directly. + +The resulting $L$ has the stored multipliers below the diagonal (with unit diagonal), and $U$ +sits on and above it. + +## Pitfalls & Edge Cases + +- **Zero pivot after partial pivoting**: if the largest candidate in a column is below the + singularity threshold $\varepsilon$, the matrix is (numerically) rank-deficient. The + factorization aborts and returns `false`; the caller must not invoke `Solve` or `Inverse`. +- **Skipping pivoting**: without row swaps, a matrix with $A_{00} = 0$ would cause immediate + division by zero; partial pivoting is mandatory for general matrices. +- **Growth factor**: partial pivoting bounds the element growth factor at $2^{n-1}$; in + practice growth is far smaller. For adversarial matrices, complete pivoting or iterative + refinement provides stronger guarantees. +- **Reusing the factored form**: `Inverse` calls `Solve` $n$ times and is $O(n^3)$; prefer + `Solve` for a single right-hand side to avoid the extra $O(n^3)$ work and the conditioning + penalty of computing the inverse explicitly. +- **Symmetric positive-definite systems**: use `CholeskyDecomposition` instead — it is + roughly twice as fast and requires no pivoting. + +## Variants & Generalizations + +- **Complete pivoting (GECP)**: swaps both rows and columns to achieve a smaller growth + factor bound; rarely needed in practice but theoretically more stable. +- **Block LU**: partitions $A$ into blocks for cache-friendly BLAS-3 operations on large + matrices; impractical for the small $n$ typical in embedded applications. +- **Sparse LU (SuperLU / KLU)**: reorders columns to minimise fill-in; relevant for + finite-element or circuit matrices, not for dense embedded systems. +- **Recursive LU**: divides the matrix recursively for parallelism; incompatible with the + no-recursion constraint of this library. + +## Applications + +- General dense linear system solvers for observer and controller gains. +- One-shot or repeated-solve scenarios in model-predictive control warm-starting. +- Determinant computation for stability margins and matrix rank tests. +- Matrix inversion for covariance update steps where the explicit inverse is required. + +## Connections to Other Algorithms + +- `GaussianElimination` performs the same GEPP steps but does not retain the factored form, + making it a single-shot solver; `LuDecomposition` amortises the $O(n^3)$ cost over + multiple `Solve` calls. +- `CholeskyDecomposition` is a specialised $O\!\left(\tfrac{1}{3}n^3\right)$ factorization + for symmetric positive-definite matrices — roughly $2\times$ faster and without the + overhead of pivoting. +- `QrDecomposition` is preferred for least-squares (overdetermined) problems; for square + full-rank systems, LU is faster. +- `ConditionNumber` can be estimated cheaply after factorization by examining the ratio + $\max_k |U_{kk}| / \min_k |U_{kk}|$. + +## References & Further Reading + +- G. H. Golub & C. F. Van Loan, *Matrix Computations*, 4th ed., Johns Hopkins UP, 2013, Ch. 3. +- W. H. Press et al., *Numerical Recipes in C*, 3rd ed., Cambridge UP, 2007, §2.3. +- L. N. Trefethen & D. Bau, *Numerical Linear Algebra*, SIAM, 1997, Lecture 20–22. diff --git a/doc/solvers/README.md b/doc/solvers/README.md index 4a3a5598..1eb988ee 100644 --- a/doc/solvers/README.md +++ b/doc/solvers/README.md @@ -15,3 +15,4 @@ Numerical solvers for linear systems, polynomial roots, and matrix equations. | [Spectral Radius](SpectralRadius.md) | Dominant eigenvalue magnitude and discrete-time Schur stability margin | | [Condition Number](ConditionNumber.md) | 1-norm condition number estimate via column-wise inverse solve | | [QR Decomposition](QrDecomposition.md) | Householder factorization and Givens streaming row update for least-squares solves | +| [LU Decomposition](LuDecomposition.md) | PA = LU factorization with partial pivoting for general dense linear systems | diff --git a/numerical/math/TriangularSolve.hpp b/numerical/math/TriangularSolve.hpp index e2d4d39f..a456c179 100644 --- a/numerical/math/TriangularSolve.hpp +++ b/numerical/math/TriangularSolve.hpp @@ -13,6 +13,23 @@ namespace math { + template + [[nodiscard]] OPTIMIZE_FOR_SPEED Vector SolveUnitLowerTriangular(const Matrix& l, const Vector& c) + { + Vector y{}; + + for (std::size_t i = 0; i < N; ++i) + { + T sum = c.at(i, 0); + for (std::size_t j = 0; j < i; ++j) + sum = sum - l.at(i, j) * y.at(j, 0); + + y.at(i, 0) = sum; + } + + return y; + } + template [[nodiscard]] OPTIMIZE_FOR_SPEED Vector SolveUpperTriangular(const Matrix& r, const Vector& c) { diff --git a/numerical/math/test/TestTriangularSolve.cpp b/numerical/math/test/TestTriangularSolve.cpp index 7e8288b7..3af8803d 100644 --- a/numerical/math/test/TestTriangularSolve.cpp +++ b/numerical/math/test/TestTriangularSolve.cpp @@ -33,6 +33,29 @@ TEST_F(TriangularSolveTest, SolvesUpperTriangularSystem) EXPECT_NEAR(solution.at(i, 0), x.at(i, 0), math::Tolerance()); } +TEST_F(TriangularSolveTest, SolvesUnitLowerTriangularSystem) +{ + math::Matrix l{ + { 1.0f, 0.0f, 0.0f }, + { -2.0f, 1.0f, 0.0f }, + { 3.0f, 4.0f, 1.0f } + }; + + math::Vector c; + for (std::size_t i = 0; i < 3; ++i) + { + float sum = 0.0f; + for (std::size_t j = 0; j < 3; ++j) + sum += l.at(i, j) * x.at(j, 0); + c.at(i, 0) = sum; + } + + auto solution = math::SolveUnitLowerTriangular(l, c); + + for (std::size_t i = 0; i < 3; ++i) + EXPECT_NEAR(solution.at(i, 0), x.at(i, 0), math::Tolerance()); +} + TEST_F(TriangularSolveTest, IdentityReturnsRightHandSide) { auto identity = math::SquareMatrix::Identity(); diff --git a/numerical/solvers/CMakeLists.txt b/numerical/solvers/CMakeLists.txt index 697a8007..faac6736 100644 --- a/numerical/solvers/CMakeLists.txt +++ b/numerical/solvers/CMakeLists.txt @@ -18,6 +18,7 @@ target_sources(numerical.solver PRIVATE DurandKerner.hpp GaussianElimination.hpp LevinsonDurbin.hpp + LuDecomposition.hpp OdeSystem.hpp QrDecomposition.hpp RungeKuttaIntegrators.hpp @@ -29,6 +30,7 @@ numerical_add_coverage_sources(numerical.solver DiscreteAlgebraicRiccatiEquation.cpp DurandKerner.cpp GaussianElimination.cpp + LuDecomposition.cpp QrDecomposition.cpp RungeKuttaIntegrators.cpp SpectralRadius.cpp diff --git a/numerical/solvers/LuDecomposition.cpp b/numerical/solvers/LuDecomposition.cpp new file mode 100644 index 00000000..302ece81 --- /dev/null +++ b/numerical/solvers/LuDecomposition.cpp @@ -0,0 +1,7 @@ +#include "numerical/solvers/LuDecomposition.hpp" + +namespace solvers +{ + template class LuDecomposition; + template class LuDecomposition; +} diff --git a/numerical/solvers/LuDecomposition.hpp b/numerical/solvers/LuDecomposition.hpp new file mode 100644 index 00000000..fd87c0ce --- /dev/null +++ b/numerical/solvers/LuDecomposition.hpp @@ -0,0 +1,169 @@ +#pragma once + +#if defined(__GNUC__) || defined(__clang__) +#pragma GCC optimize("O3", "fast-math") +#endif + +#include "numerical/math/CompilerOptimizations.hpp" +#include "numerical/math/Matrix.hpp" +#include "numerical/math/TriangularSolve.hpp" +#include +#include +#include +#include + +namespace solvers +{ + template + class LuDecomposition + { + static_assert(std::is_floating_point_v, "LuDecomposition supports floating-point types only"); + + public: + LuDecomposition() = default; + + OPTIMIZE_FOR_SPEED bool Decompose(const math::Matrix& a); + OPTIMIZE_FOR_SPEED math::Vector Solve(const math::Vector& b) const; + math::Matrix Inverse() const; + T Determinant() const; + + bool IsSingular() const { return singular; } + math::Matrix L() const; + math::Matrix U() const; + math::Matrix P() const; + + private: + math::Matrix lu{}; + std::array piv{}; + int pivotSign{ 1 }; + bool singular{ false }; + }; + + template + OPTIMIZE_FOR_SPEED bool LuDecomposition::Decompose(const math::Matrix& a) + { + lu = a; + for (std::size_t i = 0; i < N; ++i) + piv[i] = i; + pivotSign = 1; + singular = false; + + constexpr T eps{ T{ 1e-12f } }; + + for (std::size_t k = 0; k < N; ++k) + { + std::size_t p = k; + T maxVal = std::abs(lu.at(k, k)); + for (std::size_t i = k + 1; i < N; ++i) + { + T candidate = std::abs(lu.at(i, k)); + if (candidate > maxVal) + { + maxVal = candidate; + p = i; + } + } + + if (maxVal < eps) + { + singular = true; + return false; + } + + if (p != k) + { + for (std::size_t j = 0; j < N; ++j) + { + T tmp = lu.at(p, j); + lu.at(p, j) = lu.at(k, j); + lu.at(k, j) = tmp; + } + std::size_t tmpIdx = piv[p]; + piv[p] = piv[k]; + piv[k] = tmpIdx; + pivotSign = -pivotSign; + } + + for (std::size_t i = k + 1; i < N; ++i) + { + lu.at(i, k) /= lu.at(k, k); + for (std::size_t j = k + 1; j < N; ++j) + lu.at(i, j) -= lu.at(i, k) * lu.at(k, j); + } + } + + return true; + } + + template + OPTIMIZE_FOR_SPEED math::Vector LuDecomposition::Solve(const math::Vector& b) const + { + math::Vector pb{}; + for (std::size_t i = 0; i < N; ++i) + pb.at(i, 0) = b.at(piv[i], 0); + + math::Vector y = math::SolveUnitLowerTriangular(lu, pb); + return math::SolveUpperTriangular(lu, y); + } + + template + math::Matrix LuDecomposition::Inverse() const + { + math::Matrix result{}; + for (std::size_t j = 0; j < N; ++j) + { + math::Vector e{}; + e.at(j, 0) = T{ 1 }; + auto col = Solve(e); + for (std::size_t i = 0; i < N; ++i) + result.at(i, j) = col.at(i, 0); + } + return result; + } + + template + T LuDecomposition::Determinant() const + { + T det{ static_cast(pivotSign) }; + for (std::size_t k = 0; k < N; ++k) + det *= lu.at(k, k); + return det; + } + + template + math::Matrix LuDecomposition::L() const + { + math::Matrix result{}; + for (std::size_t i = 0; i < N; ++i) + { + result.at(i, i) = T{ 1 }; + for (std::size_t j = 0; j < i; ++j) + result.at(i, j) = lu.at(i, j); + } + return result; + } + + template + math::Matrix LuDecomposition::U() const + { + math::Matrix result{}; + for (std::size_t i = 0; i < N; ++i) + for (std::size_t j = i; j < N; ++j) + result.at(i, j) = lu.at(i, j); + return result; + } + + template + math::Matrix LuDecomposition::P() const + { + math::Matrix result{}; + for (std::size_t i = 0; i < N; ++i) + result.at(i, piv[i]) = T{ 1 }; + return result; + } + +#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD + extern template class LuDecomposition; + extern template class LuDecomposition; +#endif +} diff --git a/numerical/solvers/test/CMakeLists.txt b/numerical/solvers/test/CMakeLists.txt index dc8d0182..7d8d2754 100644 --- a/numerical/solvers/test/CMakeLists.txt +++ b/numerical/solvers/test/CMakeLists.txt @@ -13,6 +13,7 @@ target_sources(numerical.solvers_test PRIVATE TestDurandKerner.cpp TestGaussianElimination.cpp TestLevinsonDurbin.cpp + TestLuDecomposition.cpp TestQrDecomposition.cpp TestRungeKuttaIntegrators.cpp TestSpectralRadius.cpp diff --git a/numerical/solvers/test/TestLuDecomposition.cpp b/numerical/solvers/test/TestLuDecomposition.cpp new file mode 100644 index 00000000..5c037592 --- /dev/null +++ b/numerical/solvers/test/TestLuDecomposition.cpp @@ -0,0 +1,159 @@ +#include "numerical/math/Tolerance.hpp" +#include "numerical/solvers/GaussianElimination.hpp" +#include "numerical/solvers/LuDecomposition.hpp" +#include + +namespace +{ + class TestLuDecomposition : public ::testing::Test + { + protected: + solvers::LuDecomposition lu{}; + solvers::LuDecomposition lu4{}; + }; +} + +TEST_F(TestLuDecomposition, reconstructs_PA_equals_LU) +{ + math::Matrix a{ + { 2.0f, 1.0f, -1.0f }, + { -3.0f, -1.0f, 2.0f }, + { -2.0f, 1.0f, 2.0f } + }; + + lu.Decompose(a); + + auto pa = lu.P() * a; + auto reconstructed = lu.L() * lu.U(); + + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + EXPECT_NEAR(pa.at(i, j), reconstructed.at(i, j), 1e-5f); +} + +TEST_F(TestLuDecomposition, solves_linear_system) +{ + math::Matrix a{ + { 2.0f, 1.0f, -1.0f }, + { -3.0f, -1.0f, 2.0f }, + { -2.0f, 1.0f, 2.0f } + }; + + math::Vector b{ { 8.0f }, { -11.0f }, { -3.0f } }; + + lu.Decompose(a); + auto x = lu.Solve(b); + + EXPECT_NEAR(x.at(0, 0), 2.0f, 1e-5f); + EXPECT_NEAR(x.at(1, 0), 3.0f, 1e-5f); + EXPECT_NEAR(x.at(2, 0), -1.0f, 1e-5f); +} + +TEST_F(TestLuDecomposition, determinant_matches_known) +{ + math::Matrix a{ + { 2.0f, 1.0f, -1.0f }, + { -3.0f, -1.0f, 2.0f }, + { -2.0f, 1.0f, 2.0f } + }; + + lu.Decompose(a); + + EXPECT_NEAR(lu.Determinant(), -1.0f, 1e-4f); +} + +TEST_F(TestLuDecomposition, inverse_times_A_is_identity) +{ + math::Matrix a{ + { 2.0f, 1.0f, -1.0f }, + { -3.0f, -1.0f, 2.0f }, + { -2.0f, 1.0f, 2.0f } + }; + + lu.Decompose(a); + auto inv = lu.Inverse(); + auto product = a * inv; + + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + EXPECT_NEAR(product.at(i, j), (i == j) ? 1.0f : 0.0f, 1e-5f); +} + +TEST_F(TestLuDecomposition, pivoting_handles_zero_leading_pivot) +{ + math::Matrix a{ + { 0.0f, 2.0f, 3.0f }, + { 4.0f, 5.0f, 6.0f }, + { 7.0f, 8.0f, 9.0f } + }; + + bool ok = lu.Decompose(a); + EXPECT_TRUE(ok); + + math::Vector b{ { 1.0f }, { 2.0f }, { 3.0f } }; + auto x = lu.Solve(b); + auto residual = a * x; + + for (std::size_t i = 0; i < 3; ++i) + EXPECT_NEAR(residual.at(i, 0), b.at(i, 0), 1e-4f); +} + +TEST_F(TestLuDecomposition, singular_matrix_returns_false) +{ + math::Matrix a{ + { 1.0f, 2.0f, 3.0f }, + { 4.0f, 5.0f, 6.0f }, + { 7.0f, 8.0f, 9.0f } + }; + + bool ok = lu.Decompose(a); + + EXPECT_FALSE(ok); +} + +TEST_F(TestLuDecomposition, multiple_rhs_reuse_factorization) +{ + math::Matrix a{ + { 2.0f, 1.0f, 0.0f, 0.0f }, + { 1.0f, 3.0f, 1.0f, 0.0f }, + { 0.0f, 1.0f, 4.0f, 1.0f }, + { 0.0f, 0.0f, 1.0f, 2.0f } + }; + + math::Vector b1{ { 1.0f }, { 2.0f }, { 3.0f }, { 4.0f } }; + math::Vector b2{ { 5.0f }, { 6.0f }, { 7.0f }, { 8.0f } }; + + lu4.Decompose(a); + auto x1 = lu4.Solve(b1); + auto x2 = lu4.Solve(b2); + + EXPECT_NEAR(x1.at(0, 0), 0.2258065f, 1e-4f); + EXPECT_NEAR(x1.at(1, 0), 0.5483871f, 1e-4f); + EXPECT_NEAR(x1.at(2, 0), 0.1290323f, 1e-4f); + EXPECT_NEAR(x1.at(3, 0), 1.9354839f, 1e-4f); + + EXPECT_NEAR(x2.at(0, 0), 1.9032258f, 1e-4f); + EXPECT_NEAR(x2.at(1, 0), 1.1935484f, 1e-4f); + EXPECT_NEAR(x2.at(2, 0), 0.5161290f, 1e-4f); + EXPECT_NEAR(x2.at(3, 0), 3.7419355f, 1e-4f); +} + +TEST_F(TestLuDecomposition, matches_gaussian_elimination) +{ + math::Matrix a{ + { 2.0f, 1.0f, -1.0f }, + { -3.0f, -1.0f, 2.0f }, + { -2.0f, 1.0f, 2.0f } + }; + + math::Vector b{ { 8.0f }, { -11.0f }, { -3.0f } }; + + lu.Decompose(a); + auto xLu = lu.Solve(b); + + solvers::GaussianElimination ge{}; + auto xGe = ge.Solve(a, b); + + for (std::size_t i = 0; i < 3; ++i) + EXPECT_NEAR(xLu.at(i, 0), xGe.at(i, 0), math::Tolerance()); +} diff --git a/roadmap/solvers/LuDecomposition/explanation.md b/roadmap/solvers/LuDecomposition/explanation.md deleted file mode 100644 index 2302fc59..00000000 --- a/roadmap/solvers/LuDecomposition/explanation.md +++ /dev/null @@ -1,27 +0,0 @@ -# LU Decomposition — Overview - -## What it is -A factorization `P·A = L·U` that splits a permuted matrix into a unit lower-triangular `L` and an -upper-triangular `U`. It is the general-purpose engine for dense linear systems, determinants, and -matrix inverses. - -## Why it matters (embedded) -Many estimators and controllers reduce to "solve `A x = b`" for a general (non-symmetric) `A`. -Factoring `A` **once** and reusing `L`, `U` turns every subsequent solve into two cheap triangular -sweeps — a large win when the same plant matrix is solved against many right-hand sides in a loop. - -## How it works (intuition) -It is Gaussian elimination with bookkeeping: each step eliminates a column below the pivot and stores -the multipliers in `L`. **Partial pivoting** swaps in the largest-magnitude pivot first, which keeps -the multipliers `≤ 1` and stops round-off from exploding. The permutation `P` records those swaps. - -## Key parameters -- **Pivoting strategy** — partial pivoting (row swaps) is the stable default. -- **Reuse** — decompose once, then `Solve` repeatedly; only form `Inverse` when it is truly needed. - -## Reference -Golub & Van Loan, *Matrix Computations*, Ch. 3 (Gaussian elimination with partial pivoting). - -## See also -`GaussianElimination` (single-shot solve), `CholeskyDecomposition` (SPD, faster), -`QrDecomposition` (least squares). diff --git a/roadmap/solvers/LuDecomposition/implementation.md b/roadmap/solvers/LuDecomposition/implementation.md deleted file mode 100644 index d7e9bca8..00000000 --- a/roadmap/solvers/LuDecomposition/implementation.md +++ /dev/null @@ -1,78 +0,0 @@ -# LU Decomposition with Partial Pivoting — Implementation Pseudocode - -> Roadmap ref: #28 (Tier 4) · Target: `numerical/solvers` · Namespace `solvers` · Type: `float` (templated on `T`, instantiated for `float` only) - -## Data structures - -``` -template # static_assert(std::is_floating_point_v); instantiated for float -class LuDecomposition: - Matrix lu # L (unit, below diag) + U (on/above diag) - std::array piv # row permutation P - int pivotSign # +1 / −1 for the determinant - bool singular -``` - -## Interface - -``` -LuDecomposition() -bool Decompose(const Matrix& A) # PA = LU; hot path; false if singular -Vector Solve(const Vector& b) const # hot path -Matrix Inverse() const # N column solves -T Determinant() const -``` - -## Algorithm (pseudocode) - -``` -function Decompose(A): # OPTIMIZE_FOR_SPEED (GEPP) - lu = A; piv = [0..N-1]; pivotSign = +1 - for k in 0 .. N-1: - p = argmax_{i>=k} |lu[i,k]| # partial pivot: largest magnitude - if |lu[p,k]| < eps: singular = true; return false - if p != k: swap rows k,p; swap piv[k],piv[p]; pivotSign = −pivotSign - for i in k+1 .. N-1: - lu[i,k] /= lu[k,k] # multiplier stored in L - lu[i, k+1..] -= lu[i,k] · lu[k, k+1..] # update trailing submatrix - return true - -function Solve(b): # OPTIMIZE_FOR_SPEED - y = ForwardSubstitute(L, P·b) # unit-lower solve - return BackSubstitute(U, y) # upper solve - -function Determinant(): return pivotSign · Π_k lu[k,k] -function Inverse(): solve A·X = I column by column # avoid unless truly needed -``` - -## Complexity & memory - -- `Decompose`: `O(N³)` (≈ `(2/3)N³`); each `Solve`: `O(N²)`; `Inverse`: `O(N³)` (N solves). -- Memory: `O(N²)` for the packed `LU` plus `O(N)` for the pivot vector — static, no heap. - -## Numerical / embedded notes - -- **Partial pivoting** (largest-magnitude pivot) bounds the growth factor and makes GEPP backward - stable; never skip it for a general matrix. -- Reuse the factored form and triangular solves shared with `GaussianElimination` — factor **once**, - then solve cheaply for many right-hand sides. -- `Determinant` from the product of `U`'s diagonal (times the pivot sign) avoids cofactor blow-up. -- Prefer `CholeskyDecomposition` when `A` is symmetric positive-definite (≈2× faster, no pivoting); - use LU for general / non-symmetric systems. -- Only materialize `Inverse()` when the inverse itself is the product — otherwise `Solve` is cheaper - and better conditioned. -- Float-only: `static_assert(std::is_floating_point_v)`; the generic `T` signature keeps a - `Q15`/`Q31` specialisation cheap to add later. - -## Deployment - -- Header: `numerical/solvers/LuDecomposition.hpp` — `#pragma once` → - `#pragma GCC optimize("O3","fast-math")`, `OPTIMIZE_FOR_SPEED` on `Decompose`/`Solve`, and - `extern template class LuDecomposition;` under `#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD`. -- Coverage: `numerical/solvers/LuDecomposition.cpp` → - `template class LuDecomposition;` and `` (the tested sizes). -- Test: `numerical/solvers/test/TestLuDecomposition.cpp` -- Doc: `doc/solvers/LuDecomposition.md` -- CMake: `.hpp` → `target_sources`; `.cpp` → `numerical_add_coverage_sources`; - `TestLuDecomposition.cpp` → the `_test` target. -- Generic pattern: see `roadmap/DEPLOYMENT.md`. diff --git a/roadmap/solvers/LuDecomposition/tests.md b/roadmap/solvers/LuDecomposition/tests.md deleted file mode 100644 index 5008aa67..00000000 --- a/roadmap/solvers/LuDecomposition/tests.md +++ /dev/null @@ -1,58 +0,0 @@ -# LU Decomposition — Unit Test Plan (Pseudocode) - -> GoogleTest · `TEST_F` (`float`) · `StrictMock` only · no heap. - -## Fixture - -``` -class TestLuDecomposition : public ::testing::Test: - LuDecomposition lu - LuDecomposition lu4 -# each case below is a TEST_F(TestLuDecomposition, ) -``` - -## Test cases (Arrange / Act / Assert) - -``` -reconstructs_PA_equals_LU: - Arrange: fixed non-symmetric 3×3 A - Act: Decompose(A) - Assert: P·A ≈ L·U - -solves_linear_system: - Arrange: A x = b with known x - Assert: Solve(b) ≈ x - -determinant_matches_known: - Arrange: 3×3 with hand-computed det - Assert: Determinant() ≈ reference (correct sign) - -inverse_times_A_is_identity: - Assert: A · Inverse() ≈ I - -pivoting_handles_zero_leading_pivot: - Arrange: A with a[0][0] == 0 (forces a row swap) - Assert: factorization still reconstructs P·A - -singular_matrix_returns_false: - Arrange: rank-deficient A - Assert: Decompose returns false, singular flag set - -multiple_rhs_reuse_factorization: - Arrange: factor once, solve for two different b - Assert: both solutions correct (no re-factor) - -matches_gaussian_elimination: - Assert: Solve(b) equals the solvers::GaussianElimination result (within tol) -``` - -## Reference vectors - -- 3×3 requiring a pivot swap ⇒ documented `P`, `L`, `U`, and `det`. -- Cross-check `Solve` against `GaussianElimination` on the same system. - -## Edge cases - -- Singular / zero-pivot column ⇒ `false`, no divide-by-zero. -- Near-singular (large growth) ⇒ solution still bounded with pivoting. -- Identity and `1×1` trivial systems. From d8b722d060ea105650119469c3e2151f6160a967 Mon Sep 17 00:00:00 2001 From: gfs Date: Thu, 30 Jul 2026 12:53:23 +0200 Subject: [PATCH 2/3] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- doc/solvers/LuDecomposition.md | 12 ++++++------ numerical/solvers/LuDecomposition.hpp | 6 +++++- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/doc/solvers/LuDecomposition.md b/doc/solvers/LuDecomposition.md index 7365136a..f130750d 100644 --- a/doc/solvers/LuDecomposition.md +++ b/doc/solvers/LuDecomposition.md @@ -71,12 +71,12 @@ substitution pairs — at a total cost of $O(n^3)$. ## Complexity Analysis -| Operation | Time | Space | Notes | -|--------------|-----------------------------|----------|-----------------------------------------------| -| `Decompose` | $O\!\left(\tfrac{2}{3}n^3\right)$ | $O(n^2)$ | GEPP in-place; multipliers packed with $U$ | -| `Solve` | $O(n^2)$ | $O(n)$ | Two triangular sweeps; reuses stored factors | -| `Determinant`| $O(n)$ | $O(1)$ | One pass over $U$'s diagonal | -| `Inverse` | $O(n^3)$ | $O(n^2)$ | $n$ column solves; avoid when `Solve` suffices | +| Operation | Time | Space | Notes | +|---------------|-----------------------------------|----------|------------------------------------------------| +| `Decompose` | $O\!\left(\tfrac{2}{3}n^3\right)$ | $O(n^2)$ | GEPP in-place; multipliers packed with $U$ | +| `Solve` | $O(n^2)$ | $O(n)$ | Two triangular sweeps; reuses stored factors | +| `Determinant` | $O(n)$ | $O(1)$ | One pass over $U$'s diagonal | +| `Inverse` | $O(n^3)$ | $O(n^2)$ | $n$ column solves; avoid when `Solve` suffices | All storage is stack-allocated; the packed LU matrix and pivot vector are bounded by $N$. diff --git a/numerical/solvers/LuDecomposition.hpp b/numerical/solvers/LuDecomposition.hpp index fd87c0ce..b348372d 100644 --- a/numerical/solvers/LuDecomposition.hpp +++ b/numerical/solvers/LuDecomposition.hpp @@ -27,7 +27,11 @@ namespace solvers math::Matrix Inverse() const; T Determinant() const; - bool IsSingular() const { return singular; } + bool IsSingular() const + { + return singular; + } + math::Matrix L() const; math::Matrix U() const; math::Matrix P() const; From 610a76a9e971972492ee61e2d6f05cb9b4aa64db Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Thu, 30 Jul 2026 10:54:07 +0000 Subject: [PATCH 3/3] fix unit test --- doc/solvers/LuDecomposition.md | 10 +++++++--- numerical/solvers/LuDecomposition.hpp | 8 ++++++-- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/doc/solvers/LuDecomposition.md b/doc/solvers/LuDecomposition.md index 7365136a..1ebaa533 100644 --- a/doc/solvers/LuDecomposition.md +++ b/doc/solvers/LuDecomposition.md @@ -99,9 +99,13 @@ sits on and above it. ## Pitfalls & Edge Cases -- **Zero pivot after partial pivoting**: if the largest candidate in a column is below the - singularity threshold $\varepsilon$, the matrix is (numerically) rank-deficient. The - factorization aborts and returns `false`; the caller must not invoke `Solve` or `Inverse`. +- **Zero pivot after partial pivoting**: if the largest candidate in a column falls at or + below a *relative* threshold ($\varepsilon_\text{rel}\,\max_k|A_{kk}|$, with + $\varepsilon_\text{rel} = 10^{-6}$ for `float`), the matrix is (numerically) rank-deficient. + The factorization aborts and returns `false`; the caller must not invoke `Solve` or `Inverse`. + A relative (rather than absolute) tolerance is essential: on a truly singular matrix the + trailing pivot is a rounding-error residual whose magnitude scales with the matrix entries + and varies across platforms/`fast-math` settings, so an absolute cutoff would miss it. - **Skipping pivoting**: without row swaps, a matrix with $A_{00} = 0$ would cause immediate division by zero; partial pivoting is mandatory for general matrices. - **Growth factor**: partial pivoting bounds the element growth factor at $2^{n-1}$; in diff --git a/numerical/solvers/LuDecomposition.hpp b/numerical/solvers/LuDecomposition.hpp index fd87c0ce..60ca1774 100644 --- a/numerical/solvers/LuDecomposition.hpp +++ b/numerical/solvers/LuDecomposition.hpp @@ -48,7 +48,8 @@ namespace solvers pivotSign = 1; singular = false; - constexpr T eps{ T{ 1e-12f } }; + constexpr T relativeEps{ T{ 1e-6f } }; + T maxPivot{ T{} }; for (std::size_t k = 0; k < N; ++k) { @@ -64,7 +65,10 @@ namespace solvers } } - if (maxVal < eps) + if (maxVal > maxPivot) + maxPivot = maxVal; + + if (maxVal <= maxPivot * relativeEps) { singular = true; return false;