feat: add lu decomposition

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.

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`                 | ★★★★☆      |

doc/solvers/LuDecomposition.md

@@ -0,0 +1,155 @@
+# 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 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
+  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.

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       |

numerical/math/TriangularSolve.hpp

@@ -13,6 +13,23 @@
 
 namespace math
 {
+    template<typename T, std::size_t N>
+    [[nodiscard]] OPTIMIZE_FOR_SPEED Vector<T, N> SolveUnitLowerTriangular(const Matrix<T, N, N>& l, const Vector<T, N>& c)
+    {
+        Vector<T, N> 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<typename T, std::size_t N>
     [[nodiscard]] OPTIMIZE_FOR_SPEED Vector<T, N> SolveUpperTriangular(const Matrix<T, N, N>& r, const Vector<T, N>& c)
     {

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<float>());
 }
 
+TEST_F(TriangularSolveTest, SolvesUnitLowerTriangularSystem)
+{
+    math::Matrix<float, 3, 3> l{
+        { 1.0f, 0.0f, 0.0f },
+        { -2.0f, 1.0f, 0.0f },
+        { 3.0f, 4.0f, 1.0f }
+    };
+
+    math::Vector<float, 3> 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<float>());
+}
+
 TEST_F(TriangularSolveTest, IdentityReturnsRightHandSide)
 {
     auto identity = math::SquareMatrix<float, 3>::Identity();

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

numerical/solvers/LuDecomposition.cpp

@@ -0,0 +1,7 @@
+#include "numerical/solvers/LuDecomposition.hpp"
+
+namespace solvers
+{
+    template class LuDecomposition<float, 3>;
+    template class LuDecomposition<float, 4>;
+}