feat: add svd algo

README.md

@@ -25,7 +25,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, MatrixExponential                             |
-| [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 |
+| [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, Singular Value Decomposition (Golub-Kahan) |
 | [Nonlinear Control](doc/nonlinear_control/README.md)               | Feedback Linearization, Backstepping Control |
 | [Robust Control](doc/robust_control/README.md)                     | Active Disturbance Rejection Control (ADRC + ESO), Sliding Mode Control (SMC), Disturbance Observer (DOB) |
 | [Performance Optimization](doc/performance-optimization/README.md) | Compiler optimizations, SIMD                                         |

ROADMAP.md

@@ -26,7 +26,6 @@ Difficulty legend:
 
 | #  | Component                                            | Target module             | Difficulty |
 |----|------------------------------------------------------|---------------------------|------------|
-| 43 | Singular Value Decomposition (Golub-Kahan)           | `solvers`                 | ★★★★★      |
 | 44 | Total Least Squares                                  | `estimators/offline`      | ★★★★★      |
 | 45 | IIR filter design (Butterworth/Chebyshev + bilinear) | `filters/passive`         | ★★★★★      |
 | 46 | H∞ state-feedback control                            | `robust_control` (new)    | ★★★★★      |

doc/solvers/README.md

@@ -17,4 +17,5 @@ Numerical solvers for linear systems, polynomial roots, and matrix equations.
 | [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       |
 | [Lyapunov / Sylvester Solvers](LyapunovSylvester.md)                       | Sylvester AX+XB=C and continuous/discrete Lyapunov solvers via Kronecker vectorisation |
-| [Jacobi Eigenvalue Solver](JacobiEigenSolver.md)                           | Cyclic Jacobi rotations for the full symmetric eigenvalue/eigenvector problem          |
+| [Singular Value Decomposition](SingularValueDecomposition.md)              | Golub-Kahan bidiagonalization and implicit-shift QR for SVD, pseudo-inverse, rank, and condition number |
+| [Jacobi Eigenvalue Solver](JacobiEigenSolver.md)                           | Cyclic Jacobi rotations for the full symmetric eigenvalue/eigenvector problem                           |

doc/solvers/SingularValueDecomposition.md

@@ -0,0 +1,133 @@
+# Singular Value Decomposition
+
+## Overview & Motivation
+
+The singular value decomposition (SVD) factors any $m \times n$ matrix $A$ ($m \geq n$) into
+two orthonormal rotation factors and a non-negative diagonal, revealing the fundamental geometric
+action of $A$ as a stretch along $n$ orthogonal directions. It is the most numerically stable
+matrix factorization available and exposes rank, condition number, and the pseudo-inverse directly
+from its output. In embedded systems, SVD powers robust least-squares solvers for sensor calibration,
+redundant-actuator control allocation, and low-rank model reduction — all applications where the
+matrix $A$ may be ill-conditioned or rank-deficient and normal-equation approaches fail.
+
+## Mathematical Theory
+
+### Factorization
+
+For $A \in \mathbb{R}^{m \times n}$ with $m \geq n$, the thin SVD is
+
+$$A = U \Sigma V^\top,$$
+
+where $U \in \mathbb{R}^{m \times n}$ has orthonormal columns ($U^\top U = I_n$),
+$\Sigma = \mathrm{diag}(\sigma_1, \ldots, \sigma_n)$ with $\sigma_1 \geq \cdots \geq \sigma_n \geq 0$,
+and $V \in \mathbb{R}^{n \times n}$ is orthogonal ($V^\top V = V V^\top = I_n$).
+The scalars $\sigma_i$ are the **singular values** of $A$; the columns of $U$ and $V$ are
+the left and right **singular vectors**.
+
+### Connection to the Symmetric Eigenproblem
+
+The singular values satisfy $\sigma_i^2 = \lambda_i(A^\top A)$, where $\lambda_i$ denotes the
+$i$-th eigenvalue of $A^\top A$. This ties the SVD to the symmetric eigenproblem and provides
+a cross-check: computing the eigenvalues of $A^\top A$ via the Jacobi or tridiagonal-QR solver
+should yield the squared singular values.
+
+### Pseudo-Inverse
+
+The Moore-Penrose pseudo-inverse is
+
+$$A^+ = V \Sigma^+ U^\top, \quad \Sigma^+_{ii} = \begin{cases} 1/\sigma_i & \sigma_i > \tau \\ 0 & \text{otherwise} \end{cases}$$
+
+for a threshold $\tau > 0$. Truncating small singular values regularises the solution against
+near-null-space components, avoiding division by near-zero values that arise from rank-deficiency
+or numerical noise.
+
+### Rank and Condition Number
+
+$$\mathrm{rank}(A) = \#\{i : \sigma_i > \tau\}, \quad \kappa(A) = \sigma_1 / \sigma_n.$$
+
+The condition number $\kappa$ measures sensitivity: a system with $\kappa \gg 1$ amplifies
+input perturbations by a factor of $\kappa$ in the least-squares solution.
+
+## Complexity Analysis
+
+| Phase             | Time                                            | Space    | Notes                                        |
+|-------------------|-------------------------------------------------|----------|----------------------------------------------|
+| Bidiagonalization | $O(mn^2 - n^3/3)$                               | $O(mn)$  | Two-sided Householder reflectors             |
+| QR sweeps         | $O(n^2)$ per iteration, $O(n)$ total iterations | $O(n^2)$ | Golub-Reinsch implicit-shift; converges fast |
+| Pseudo-inverse    | $O(n^2 m)$                                      | $O(nm)$  | Matrix triple product $V \Sigma^+ U^\top$    |
+| Total             | $O(mn^2 + n^3)$                                 | $O(mn)$  | Dominated by bidiagonalization               |
+
+All storage is stack-allocated; no heap is used. The bidiagonal form is stored implicitly in
+the working copy of $A$ alongside the Householder reflector scalars.
+
+## Step-by-Step Walkthrough
+
+**Input:** $A = \begin{pmatrix} 3 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 2 \end{pmatrix}$ (diagonal matrix).
+
+**Phase 1 — Bidiagonalization:** $A$ is already diagonal (and hence upper-bidiagonal). All
+Householder reflectors are trivial ($\beta = 0$). The bidiagonal form retains $B = A$;
+$U = I_3$, $V = I_3$.
+
+**Phase 2 — QR sweeps:** The superdiagonal is already zero, so the sweep loop exits immediately.
+
+**Sign correction and sort:** $\sigma = (3, 1, 2)$ after reading the diagonal; the entry
+$-1$ has a negative sign, so it is flipped to $+1$ and the corresponding column of $U$ is
+negated. After descending sort: $\sigma = (3, 2, 1)$.
+
+## Pitfalls & Edge Cases
+
+- **Rank deficiency:** a singular value of exactly zero is harmless; the pseudo-inverse
+  threshold $\tau$ must be chosen relative to $\sigma_1$ and machine epsilon to avoid
+  false zero-detection.
+- **Ill-conditioning in the shift:** the implicit shift $\mu$ is derived from the trailing
+  $2 \times 2$ submatrix; a zero superdiagonal entry deflates the problem before the shift
+  is computed, preventing division by zero.
+- **Sign ambiguity:** singular vectors are defined up to sign; the implementation folds negative
+  diagonal entries into the sign of the corresponding $U$ column so that all $\sigma_i \geq 0$.
+- **Repeated singular values:** the corresponding singular-vector subspaces are correct but the
+  individual vectors are not unique; tests should check subspace properties, not individual vectors.
+- **Condition number with zero $\sigma_{\min}$:** a rank-deficient matrix has $\sigma_n = 0$;
+  `ConditionNumber()` returns zero as a sentinel in this case rather than dividing by zero.
+
+## Variants & Generalizations
+
+- **Full SVD:** $U \in \mathbb{R}^{m \times m}$ (all $m$ left singular vectors); the extra
+  $m - n$ columns span the left null space of $A$. The thin SVD suffices for least squares.
+- **Truncated SVD:** retain only the $k$ largest singular values for low-rank approximation
+  $A \approx U_k \Sigma_k V_k^\top$; the approximation error in Frobenius norm equals
+  $\sqrt{\sigma_{k+1}^2 + \cdots + \sigma_n^2}$ (Eckart-Young theorem).
+- **Divide-and-conquer SVD:** splits the bidiagonal matrix recursively; $O(n^2)$ per level
+  and $O(n \log n)$ overall, faster than Golub-Reinsch for large $n$ but incompatible with
+  a no-recursion embedded requirement.
+- **Jacobi SVD:** applies one-sided plane rotations directly to $A$; simpler but $O(n^3)$
+  per sweep with slower convergence than implicit-shift QR.
+
+## Applications
+
+- **Least-squares regression:** $x = A^+ b$ minimises $\|Ax - b\|_2$ and is unique among
+  minimum-norm solutions when $A$ is rank-deficient.
+- **Control allocation:** maps actuator space to output space; the pseudo-inverse distributes
+  effort across redundant actuators while respecting the null space.
+- **System identification:** SVD of a Hankel matrix (Ho-Kalman) recovers system order (rank)
+  and state-space realization from impulse-response data.
+- **Principal-component analysis (offline):** left singular vectors of the centred data matrix
+  are the principal components; singular values are the standard deviations along each axis.
+- **Condition monitoring:** a rising condition number $\kappa$ before solving warns that the
+  system is near-singular and the solution may be unreliable.
+
+## Connections to Other Algorithms
+
+- `QrDecomposition` provides the Householder reflectors reused in the bidiagonalization phase.
+- `JacobiEigenSolver` (symmetric eigenproblem) is an alternative route to the same result via
+  $\sigma_i^2 = \lambda_i(A^\top A)$; the two cross-check each other.
+- `LuDecomposition` and `GaussianElimination` solve square systems but cannot handle
+  rank-deficient or overdetermined cases; SVD supersedes them for those inputs.
+- `TotalLeastSquares` (item 44) builds directly on SVD: the solution lies in the right singular
+  vector corresponding to the smallest singular value of the augmented matrix $[A\ b]$.
+
+## References & Further Reading
+
+- G. Golub and W. Kahan, "Calculating the Singular Values and Pseudo-Inverse of a Matrix," *SIAM J. Numer. Anal.*, 2(2), 1965.
+- G. Golub and C. Reinsch, "Singular Value Decomposition and Least Squares Solutions," *Numer. Math.*, 14(5), 1970.
+- G. H. Golub and C. F. Van Loan, *Matrix Computations*, 4th ed., Johns Hopkins UP, 2013, Ch. 6 & 8.
+- W. H. Press et al., *Numerical Recipes in C*, 3rd ed., Cambridge UP, 2007, §2.6.

numerical/math/HouseholderTransform.hpp

@@ -44,4 +44,40 @@ namespace math
         for (std::size_t i = start + 1; i < N; ++i)
             v.at(i, 0) = x.at(i, 0) * invV0;
     }
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    OPTIMIZE_FOR_SPEED void ApplyReflectorLeft(Matrix<T, Rows, Cols>& a, const Vector<T, Rows>& v, T beta,
+        std::size_t rowStart, std::size_t colStart)
+    {
+        static_assert(std::is_floating_point_v<T>, "ApplyReflectorLeft supports floating-point types only");
+
+        for (std::size_t j = colStart; j < Cols; ++j)
+        {
+            T dot{ T{} };
+            for (std::size_t i = rowStart; i < Rows; ++i)
+                dot += v.at(i, 0) * a.at(i, j);
+
+            T scale = beta * dot;
+            for (std::size_t i = rowStart; i < Rows; ++i)
+                a.at(i, j) -= scale * v.at(i, 0);
+        }
+    }
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    OPTIMIZE_FOR_SPEED void ApplyReflectorRight(Matrix<T, Rows, Cols>& a, const Vector<T, Cols>& v, T beta,
+        std::size_t rowStart, std::size_t colStart)
+    {
+        static_assert(std::is_floating_point_v<T>, "ApplyReflectorRight supports floating-point types only");
+
+        for (std::size_t i = rowStart; i < Rows; ++i)
+        {
+            T dot{ T{} };
+            for (std::size_t j = colStart; j < Cols; ++j)
+                dot += a.at(i, j) * v.at(j, 0);
+
+            T scale = beta * dot;
+            for (std::size_t j = colStart; j < Cols; ++j)
+                a.at(i, j) -= scale * v.at(j, 0);
+        }
+    }
 }

numerical/solvers/CMakeLists.txt

@@ -24,6 +24,7 @@ target_sources(numerical.solver PRIVATE
     OdeSystem.hpp
     QrDecomposition.hpp
     RungeKuttaIntegrators.hpp
+    SingularValueDecomposition.hpp
     Solver.hpp
     SpectralRadius.hpp
 )
@@ -37,6 +38,7 @@ numerical_add_coverage_sources(numerical.solver
     LyapunovSylvester.cpp
     QrDecomposition.cpp
     RungeKuttaIntegrators.cpp
+    SingularValueDecomposition.cpp
     SpectralRadius.cpp
 )
 

numerical/solvers/QrDecomposition.hpp

@@ -32,29 +32,11 @@ namespace solvers
         void GivensUpdateRow(const math::Matrix<T, 1, Cols>& newRow);
 
     private:
-        static void ApplyReflector(math::Matrix<T, Rows, Cols>& block, const math::Vector<T, Rows>& v, T beta, std::size_t start);
-
         math::Matrix<T, Rows, Cols> qr{};
         math::Vector<T, Cols> betas{};
         bool factored{ false };
     };
 
-    template<typename T, std::size_t Rows, std::size_t Cols>
-    OPTIMIZE_FOR_SPEED void QrDecomposition<T, Rows, Cols>::ApplyReflector(
-        math::Matrix<T, Rows, Cols>& block, const math::Vector<T, Rows>& v, T beta, std::size_t start)
-    {
-        for (std::size_t j = start; j < Cols; ++j)
-        {
-            T dot{ T{} };
-            for (std::size_t i = start; i < Rows; ++i)
-                dot += v.at(i, 0) * block.at(i, j);
-
-            T scale = beta * dot;
-            for (std::size_t i = start; i < Rows; ++i)
-                block.at(i, j) -= scale * v.at(i, 0);
-        }
-    }
-
     template<typename T, std::size_t Rows, std::size_t Cols>
     OPTIMIZE_FOR_SPEED bool QrDecomposition<T, Rows, Cols>::Decompose(const math::Matrix<T, Rows, Cols>& a)
     {
@@ -78,7 +60,7 @@ namespace solvers
 
             if (beta != T{})
             {
-                ApplyReflector(qr, v, beta, k);
+                math::ApplyReflectorLeft(qr, v, beta, k, k);
 
                 for (std::size_t i = k + 1; i < Rows; ++i)
                     qr.at(i, k) = v.at(i, 0);
@@ -116,7 +98,7 @@ namespace solvers
             for (std::size_t i = col + 1; i < Rows; ++i)
                 v.at(i, 0) = qr.at(i, col);
 
-            ApplyReflector(result, v, beta, col);
+            math::ApplyReflectorLeft(result, v, beta, col, col);
         }
 
         return result;

numerical/solvers/SingularValueDecomposition.cpp

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