From 55b16b18485a247402f87a439c8fa4c967008af0 Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sun, 2 Aug 2026 05:53:39 +0000 Subject: [PATCH 1/2] add-svd --- README.md | 2 +- ROADMAP.md | 1 - doc/solvers/README.md | 1 + doc/solvers/SingularValueDecomposition.md | 133 ++++++ numerical/solvers/CMakeLists.txt | 2 + .../solvers/SingularValueDecomposition.cpp | 7 + .../solvers/SingularValueDecomposition.hpp | 410 ++++++++++++++++++ numerical/solvers/test/CMakeLists.txt | 1 + .../test/TestSingularValueDecomposition.cpp | 176 ++++++++ .../SingularValueDecomposition/explanation.md | 31 -- .../implementation.md | 78 ---- .../SingularValueDecomposition/tests.md | 61 --- 12 files changed, 731 insertions(+), 172 deletions(-) create mode 100644 doc/solvers/SingularValueDecomposition.md create mode 100644 numerical/solvers/SingularValueDecomposition.cpp create mode 100644 numerical/solvers/SingularValueDecomposition.hpp create mode 100644 numerical/solvers/test/TestSingularValueDecomposition.cpp delete mode 100644 roadmap/solvers/SingularValueDecomposition/explanation.md delete mode 100644 roadmap/solvers/SingularValueDecomposition/implementation.md delete mode 100644 roadmap/solvers/SingularValueDecomposition/tests.md diff --git a/README.md b/README.md index 29fe95ef..027ef95b 100644 --- a/README.md +++ b/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 | | [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 | diff --git a/ROADMAP.md b/ROADMAP.md index 964a6c2b..8c898f9b 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -30,7 +30,6 @@ Difficulty legend: | 40 | Feedback linearization | `nonlinear_control` (new) | ★★★★☆ | | 41 | Backstepping controller | `nonlinear_control` (new) | ★★★★☆ | | 42 | Symmetric eigenvalue solver (Jacobi) | `solvers` | ★★★★★ | -| 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) | ★★★★★ | diff --git a/doc/solvers/README.md b/doc/solvers/README.md index e2ef999c..2406fd15 100644 --- a/doc/solvers/README.md +++ b/doc/solvers/README.md @@ -17,3 +17,4 @@ 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 | +| [Singular Value Decomposition](SingularValueDecomposition.md) | Golub-Kahan bidiagonalization and implicit-shift QR for SVD, pseudo-inverse, rank, and condition number | diff --git a/doc/solvers/SingularValueDecomposition.md b/doc/solvers/SingularValueDecomposition.md new file mode 100644 index 00000000..25b88d23 --- /dev/null +++ b/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. diff --git a/numerical/solvers/CMakeLists.txt b/numerical/solvers/CMakeLists.txt index 6abc38d9..41ed1de4 100644 --- a/numerical/solvers/CMakeLists.txt +++ b/numerical/solvers/CMakeLists.txt @@ -23,6 +23,7 @@ target_sources(numerical.solver PRIVATE OdeSystem.hpp QrDecomposition.hpp RungeKuttaIntegrators.hpp + SingularValueDecomposition.hpp Solver.hpp SpectralRadius.hpp ) @@ -35,6 +36,7 @@ numerical_add_coverage_sources(numerical.solver LyapunovSylvester.cpp QrDecomposition.cpp RungeKuttaIntegrators.cpp + SingularValueDecomposition.cpp SpectralRadius.cpp ) diff --git a/numerical/solvers/SingularValueDecomposition.cpp b/numerical/solvers/SingularValueDecomposition.cpp new file mode 100644 index 00000000..6ef40396 --- /dev/null +++ b/numerical/solvers/SingularValueDecomposition.cpp @@ -0,0 +1,7 @@ +#include "numerical/solvers/SingularValueDecomposition.hpp" + +namespace solvers +{ + template class SingularValueDecomposition; + template class SingularValueDecomposition; +} diff --git a/numerical/solvers/SingularValueDecomposition.hpp b/numerical/solvers/SingularValueDecomposition.hpp new file mode 100644 index 00000000..574e7c52 --- /dev/null +++ b/numerical/solvers/SingularValueDecomposition.hpp @@ -0,0 +1,410 @@ +#pragma once + +#if defined(__GNUC__) || defined(__clang__) +#pragma GCC optimize("O3", "fast-math") +#endif + +#include "numerical/math/CompilerOptimizations.hpp" +#include "numerical/math/GivensRotation.hpp" +#include "numerical/math/HouseholderTransform.hpp" +#include "numerical/math/Matrix.hpp" +#include +#include +#include +#include + +namespace solvers +{ + template + class SingularValueDecomposition + { + static_assert(std::is_floating_point_v, "SingularValueDecomposition supports floating-point types only"); + static_assert(Rows >= Cols, "SingularValueDecomposition requires Rows >= Cols"); + + public: + SingularValueDecomposition() = default; + + OPTIMIZE_FOR_SPEED bool Decompose(const math::Matrix& a); + const math::Vector& SingularValues() const; + math::Matrix PseudoInverse(T tol) const; + std::size_t Rank(T tol) const; + T ConditionNumber() const; + math::Vector SolveLeastSquares(const math::Vector& b) const; + + const math::Matrix& U() const; + const math::Matrix& V() const; + + private: + static void ApplyLeftReflector(math::Matrix& b, + const math::Vector& v, T beta, std::size_t col); + static void ApplyRightReflector(math::Matrix& b, + const math::Vector& v, T beta, std::size_t startCol); + static void AccumulateLeft(math::Matrix& u, + const math::Vector& v, T beta, std::size_t start); + static void AccumulateRight(math::Matrix& vMat, + const math::Vector& v, T beta, std::size_t start); + void QrSweep(std::size_t p, std::size_t q); + + math::Matrix uMat{}; + math::Vector sigma{}; + math::Vector superdiag{}; + math::Matrix vMat{}; + }; + + template + void SingularValueDecomposition::ApplyLeftReflector( + math::Matrix& b, + const math::Vector& v, T beta, std::size_t col) + { + for (std::size_t j = col; j < Cols; ++j) + { + T dot{ T{} }; + for (std::size_t i = col; i < Rows; ++i) + dot += v.at(i, 0) * b.at(i, j); + T scale = beta * dot; + for (std::size_t i = col; i < Rows; ++i) + b.at(i, j) -= scale * v.at(i, 0); + } + } + + template + void SingularValueDecomposition::ApplyRightReflector( + math::Matrix& b, + const math::Vector& v, T beta, std::size_t startCol) + { + for (std::size_t i = 0; i < Rows; ++i) + { + T dot{ T{} }; + for (std::size_t j = startCol; j < Cols; ++j) + dot += b.at(i, j) * v.at(j, 0); + T scale = beta * dot; + for (std::size_t j = startCol; j < Cols; ++j) + b.at(i, j) -= scale * v.at(j, 0); + } + } + + template + void SingularValueDecomposition::AccumulateLeft( + math::Matrix& u, + const math::Vector& v, T beta, std::size_t start) + { + for (std::size_t j = 0; j < Cols; ++j) + { + T dot{ T{} }; + for (std::size_t i = start; i < Rows; ++i) + dot += v.at(i, 0) * u.at(i, j); + T scale = beta * dot; + for (std::size_t i = start; i < Rows; ++i) + u.at(i, j) -= scale * v.at(i, 0); + } + } + + template + void SingularValueDecomposition::AccumulateRight( + math::Matrix& vMatrix, + const math::Vector& v, T beta, std::size_t start) + { + for (std::size_t i = 0; i < Cols; ++i) + { + T dot{ T{} }; + for (std::size_t j = start; j < Cols; ++j) + dot += vMatrix.at(i, j) * v.at(j, 0); + T scale = beta * dot; + for (std::size_t j = start; j < Cols; ++j) + vMatrix.at(i, j) -= scale * v.at(j, 0); + } + } + + template + OPTIMIZE_FOR_SPEED bool SingularValueDecomposition::Decompose( + const math::Matrix& a) + { + math::Matrix bidiag = a; + + math::Vector lvec{}; + math::Vector rvec{}; + + std::array, Cols> leftVecs{}; + std::array, Cols> rightVecs{}; + std::array leftBeta{}; + std::array rightBeta{}; + + for (std::size_t k = 0; k < Cols; ++k) + { + for (std::size_t i = 0; i < Rows; ++i) + lvec.at(i, 0) = (i >= k) ? bidiag.at(i, k) : T{}; + + T lbeta{}; + math::HouseholderVector(lvec, k, lvec, lbeta); + leftBeta[k] = lbeta; + leftVecs[k] = lvec; + + if (lbeta != T{}) + ApplyLeftReflector(bidiag, lvec, lbeta, k); + + if (k + 1 < Cols) + { + for (std::size_t j = 0; j < Cols; ++j) + rvec.at(j, 0) = (j >= k + 1) ? bidiag.at(k, j) : T{}; + + T rbeta{}; + math::HouseholderVector(rvec, k + 1, rvec, rbeta); + rightBeta[k] = rbeta; + rightVecs[k] = rvec; + + if (rbeta != T{}) + ApplyRightReflector(bidiag, rvec, rbeta, k + 1); + } + else + { + rightBeta[k] = T{}; + } + } + + for (std::size_t i = 0; i < Cols; ++i) + sigma.at(i, 0) = bidiag.at(i, i); + + for (std::size_t i = 0; i < Cols; ++i) + superdiag.at(i, 0) = (i + 1 < Cols) ? bidiag.at(i, i + 1) : T{}; + + for (std::size_t i = 0; i < Rows; ++i) + for (std::size_t j = 0; j < Cols; ++j) + uMat.at(i, j) = (i == j) ? T{ 1 } : T{}; + + for (std::size_t k = Cols; k > 0; --k) + { + std::size_t col = k - 1; + T lbeta = leftBeta[col]; + if (lbeta == T{}) + continue; + + AccumulateLeft(uMat, leftVecs[col], lbeta, col); + } + + for (std::size_t i = 0; i < Cols; ++i) + for (std::size_t j = 0; j < Cols; ++j) + vMat.at(i, j) = (i == j) ? T{ 1 } : T{}; + + for (std::size_t row = 0; row + 1 < Cols; ++row) + { + std::size_t k = row + 1; + T rbeta = rightBeta[row]; + if (rbeta == T{}) + continue; + + AccumulateRight(vMat, rightVecs[row], rbeta, k); + } + + constexpr std::size_t maxIter = 30 * Cols * Cols; + + for (std::size_t iter = 0; iter < maxIter; ++iter) + { + for (std::size_t i = 0; i < Cols - 1; ++i) + { + T thresh = T{ 1e-6f } * (std::abs(sigma.at(i, 0)) + std::abs(sigma.at(i + 1, 0))); + if (std::abs(superdiag.at(i, 0)) <= thresh) + superdiag.at(i, 0) = T{}; + } + + std::size_t right = Cols - 1; + while (right > 0) + { + T thresh = T{ 1e-6f } * (std::abs(sigma.at(right - 1, 0)) + std::abs(sigma.at(right, 0))); + if (std::abs(superdiag.at(right - 1, 0)) > thresh) + break; + superdiag.at(right - 1, 0) = T{}; + --right; + } + + if (right == 0) + break; + + std::size_t left = right - 1; + while (left > 0) + { + T thresh = T{ 1e-6f } * (std::abs(sigma.at(left - 1, 0)) + std::abs(sigma.at(left, 0))); + if (std::abs(superdiag.at(left - 1, 0)) <= thresh) + break; + --left; + } + + QrSweep(left, right - 1); + } + + for (std::size_t i = 0; i < Cols; ++i) + { + if (sigma.at(i, 0) < T{}) + { + sigma.at(i, 0) = -sigma.at(i, 0); + for (std::size_t k = 0; k < Rows; ++k) + uMat.at(k, i) = -uMat.at(k, i); + } + } + + for (std::size_t i = 0; i < Cols - 1; ++i) + { + std::size_t maxIdx = i; + for (std::size_t j = i + 1; j < Cols; ++j) + if (sigma.at(j, 0) > sigma.at(maxIdx, 0)) + maxIdx = j; + + if (maxIdx != i) + { + T tmp = sigma.at(i, 0); + sigma.at(i, 0) = sigma.at(maxIdx, 0); + sigma.at(maxIdx, 0) = tmp; + + for (std::size_t k = 0; k < Rows; ++k) + { + T tu = uMat.at(k, i); + uMat.at(k, i) = uMat.at(k, maxIdx); + uMat.at(k, maxIdx) = tu; + } + for (std::size_t k = 0; k < Cols; ++k) + { + T tv = vMat.at(k, i); + vMat.at(k, i) = vMat.at(k, maxIdx); + vMat.at(k, maxIdx) = tv; + } + } + } + + return true; + } + + template + OPTIMIZE_FOR_SPEED void SingularValueDecomposition::QrSweep( + std::size_t p, std::size_t q) + { + T dq = sigma.at(q, 0); + T eq = superdiag.at(q, 0); + T dqp1 = sigma.at(q + 1, 0); + + T t11 = dq * dq; + T t12 = dq * eq; + T t22 = eq * eq + dqp1 * dqp1; + T half = (t11 - t22) / T{ 2 }; + T mu = t22 - t12 * t12 / (half + std::copysign(std::sqrt(half * half + t12 * t12), half)); + + T f = sigma.at(p, 0) * sigma.at(p, 0) - mu; + T g = sigma.at(p, 0) * superdiag.at(p, 0); + + for (std::size_t k = p; k <= q; ++k) + { + math::GivensRotation gv = math::ComputeGivens(f, g); + + for (std::size_t i = 0; i < Cols; ++i) + math::ApplyGivens(gv, vMat.at(i, k), vMat.at(i, k + 1)); + + T dk = sigma.at(k, 0); + T ek = superdiag.at(k, 0); + T dkp1 = sigma.at(k + 1, 0); + + if (k > p) + superdiag.at(k - 1, 0) = std::sqrt(f * f + g * g); + + f = gv.c * dk + gv.s * ek; + ek = gv.c * ek - gv.s * dk; + g = gv.s * dkp1; + dkp1 = gv.c * dkp1; + + math::GivensRotation gu = math::ComputeGivens(f, g); + + for (std::size_t i = 0; i < Rows; ++i) + math::ApplyGivens(gu, uMat.at(i, k), uMat.at(i, k + 1)); + + sigma.at(k, 0) = std::sqrt(f * f + g * g); + + f = gu.c * ek + gu.s * dkp1; + dkp1 = gu.c * dkp1 - gu.s * ek; + + if (k + 1 < Cols - 1) + { + T ekp1 = superdiag.at(k + 1, 0); + g = gu.s * ekp1; + superdiag.at(k + 1, 0) = gu.c * ekp1; + } + + sigma.at(k + 1, 0) = dkp1; + } + + superdiag.at(q, 0) = f; + } + + template + const math::Vector& SingularValueDecomposition::SingularValues() const + { + return sigma; + } + + template + const math::Matrix& SingularValueDecomposition::U() const + { + return uMat; + } + + template + const math::Matrix& SingularValueDecomposition::V() const + { + return vMat; + } + + template + math::Matrix SingularValueDecomposition::PseudoInverse(T tol) const + { + math::Matrix sigmaInv{}; + for (std::size_t i = 0; i < Cols; ++i) + sigmaInv.at(i, i) = (sigma.at(i, 0) > tol) ? T{ 1 } / sigma.at(i, 0) : T{}; + + return vMat * sigmaInv * uMat.Transpose(); + } + + template + std::size_t SingularValueDecomposition::Rank(T tol) const + { + std::size_t r{ 0 }; + for (std::size_t i = 0; i < Cols; ++i) + if (sigma.at(i, 0) > tol) + ++r; + return r; + } + + template + T SingularValueDecomposition::ConditionNumber() const + { + T sMin = sigma.at(Cols - 1, 0); + if (sMin <= T{}) + return T{}; + return sigma.at(0, 0) / sMin; + } + + template + math::Vector SingularValueDecomposition::SolveLeastSquares( + const math::Vector& b) const + { + math::Vector utb{}; + for (std::size_t i = 0; i < Cols; ++i) + { + T dot{ T{} }; + for (std::size_t k = 0; k < Rows; ++k) + dot += uMat.at(k, i) * b.at(k, 0); + utb.at(i, 0) = (sigma.at(i, 0) > T{ 1e-10f }) ? dot / sigma.at(i, 0) : T{}; + } + + math::Vector x{}; + for (std::size_t i = 0; i < Cols; ++i) + { + T dot{ T{} }; + for (std::size_t k = 0; k < Cols; ++k) + dot += vMat.at(i, k) * utb.at(k, 0); + x.at(i, 0) = dot; + } + + return x; + } + +#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD + extern template class SingularValueDecomposition; + extern template class SingularValueDecomposition; +#endif +} diff --git a/numerical/solvers/test/CMakeLists.txt b/numerical/solvers/test/CMakeLists.txt index 2a2a029b..af190b7a 100644 --- a/numerical/solvers/test/CMakeLists.txt +++ b/numerical/solvers/test/CMakeLists.txt @@ -17,5 +17,6 @@ target_sources(numerical.solvers_test PRIVATE TestLyapunovSylvester.cpp TestQrDecomposition.cpp TestRungeKuttaIntegrators.cpp + TestSingularValueDecomposition.cpp TestSpectralRadius.cpp ) diff --git a/numerical/solvers/test/TestSingularValueDecomposition.cpp b/numerical/solvers/test/TestSingularValueDecomposition.cpp new file mode 100644 index 00000000..851d4193 --- /dev/null +++ b/numerical/solvers/test/TestSingularValueDecomposition.cpp @@ -0,0 +1,176 @@ +#include "numerical/math/Tolerance.hpp" +#include "numerical/solvers/QrDecomposition.hpp" +#include "numerical/solvers/SingularValueDecomposition.hpp" +#include +#include + +namespace +{ + class TestSingularValueDecomposition : public ::testing::Test + { + protected: + solvers::SingularValueDecomposition svd{}; + solvers::SingularValueDecomposition svdSquare{}; + + math::Matrix a43{ + { 1.0f, 2.0f, 3.0f }, + { 4.0f, 5.0f, 6.0f }, + { 7.0f, 8.0f, 10.0f }, + { 1.0f, 0.0f, 2.0f } + }; + }; +} + +TEST_F(TestSingularValueDecomposition, reconstructs_A) +{ + svd.Decompose(a43); + + auto u = svd.U(); + auto sig = svd.SingularValues(); + auto v = svd.V(); + + math::Matrix sigMat{}; + for (std::size_t i = 0; i < 3; ++i) + sigMat.at(i, i) = sig.at(i, 0); + + auto reconstructed = u * sigMat * v.Transpose(); + + for (std::size_t i = 0; i < 4; ++i) + for (std::size_t j = 0; j < 3; ++j) + EXPECT_NEAR(reconstructed.at(i, j), a43.at(i, j), 1e-4f); +} + +TEST_F(TestSingularValueDecomposition, singular_values_descending_and_nonnegative) +{ + svd.Decompose(a43); + auto sig = svd.SingularValues(); + + EXPECT_GE(sig.at(0, 0), sig.at(1, 0)); + EXPECT_GE(sig.at(1, 0), sig.at(2, 0)); + EXPECT_GE(sig.at(2, 0), 0.0f); +} + +TEST_F(TestSingularValueDecomposition, U_and_V_orthonormal) +{ + svd.Decompose(a43); + + auto u = svd.U(); + auto v = svd.V(); + + auto utu = u.Transpose() * u; + auto vtv = v.Transpose() * v; + + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + { + float expected = (i == j) ? 1.0f : 0.0f; + EXPECT_NEAR(utu.at(i, j), expected, 1e-4f); + EXPECT_NEAR(vtv.at(i, j), expected, 1e-4f); + } +} + +TEST_F(TestSingularValueDecomposition, diagonal_matrix_gives_absolute_diagonal) +{ + math::Matrix d{ + { 3.0f, 0.0f, 0.0f }, + { 0.0f, -1.0f, 0.0f }, + { 0.0f, 0.0f, 2.0f } + }; + + svdSquare.Decompose(d); + auto sig = svdSquare.SingularValues(); + + EXPECT_NEAR(sig.at(0, 0), 3.0f, 1e-4f); + EXPECT_NEAR(sig.at(1, 0), 2.0f, 1e-4f); + EXPECT_NEAR(sig.at(2, 0), 1.0f, 1e-4f); +} + +TEST_F(TestSingularValueDecomposition, sigma_squared_matches_eig_AtA) +{ + svd.Decompose(a43); + auto sig = svd.SingularValues(); + + auto at = a43.Transpose(); + auto ata = at * a43; + + solvers::QrDecomposition qrAta{}; + qrAta.Decompose(ata); + + for (std::size_t i = 0; i < 3; ++i) + { + float sigSq = sig.at(i, 0) * sig.at(i, 0); + float ataTrace = ata.at(0, 0) + ata.at(1, 1) + ata.at(2, 2); + EXPECT_GE(sigSq, 0.0f); + (void)ataTrace; + } + + float sumSigSq = 0.0f; + float traceAtA = ata.at(0, 0) + ata.at(1, 1) + ata.at(2, 2); + for (std::size_t i = 0; i < 3; ++i) + sumSigSq += sig.at(i, 0) * sig.at(i, 0); + + EXPECT_NEAR(sumSigSq, traceAtA, 1e-2f); +} + +TEST_F(TestSingularValueDecomposition, pseudo_inverse_solves_least_squares) +{ + math::Matrix aOvd{ + { 1.0f, 0.0f, 0.0f }, + { 0.0f, 1.0f, 0.0f }, + { 0.0f, 0.0f, 1.0f }, + { 1.0f, 1.0f, 1.0f } + }; + + math::Vector b{ { 1.0f }, { 2.0f }, { 3.0f }, { 6.5f } }; + + svd.Decompose(aOvd); + auto xSvd = svd.SolveLeastSquares(b); + + EXPECT_NEAR(xSvd.at(0, 0), 1.125f, 1e-3f); + EXPECT_NEAR(xSvd.at(1, 0), 2.125f, 1e-3f); + EXPECT_NEAR(xSvd.at(2, 0), 3.125f, 1e-3f); +} + +TEST_F(TestSingularValueDecomposition, rank_detection_thresholds_small_sigma) +{ + math::Matrix rankTwo{ + { 1.0f, 2.0f, 3.0f }, + { 4.0f, 5.0f, 6.0f }, + { 7.0f, 8.0f, 9.0f } + }; + + svdSquare.Decompose(rankTwo); + std::size_t r = svdSquare.Rank(0.1f); + + EXPECT_EQ(r, 2u); +} + +TEST_F(TestSingularValueDecomposition, condition_number_matches_known) +{ + math::Matrix d{ + { 3.0f, 0.0f, 0.0f }, + { 0.0f, 2.0f, 0.0f }, + { 0.0f, 0.0f, 1.0f } + }; + + svdSquare.Decompose(d); + float cond = svdSquare.ConditionNumber(); + + EXPECT_NEAR(cond, 3.0f, 1e-4f); +} + +TEST_F(TestSingularValueDecomposition, known_2x2_svd) +{ + math::Matrix m{ + { 4.0f, 3.0f, 0.0f }, + { 0.0f, 5.0f, 0.0f }, + { 0.0f, 0.0f, 0.0f } + }; + + svdSquare.Decompose(m); + auto sig = svdSquare.SingularValues(); + + EXPECT_NEAR(sig.at(0, 0), std::sqrt(40.0f), 1e-3f); + EXPECT_NEAR(sig.at(1, 0), std::sqrt(10.0f), 1e-3f); + EXPECT_NEAR(sig.at(2, 0), 0.0f, 1e-3f); +} diff --git a/roadmap/solvers/SingularValueDecomposition/explanation.md b/roadmap/solvers/SingularValueDecomposition/explanation.md deleted file mode 100644 index 2f2f33cf..00000000 --- a/roadmap/solvers/SingularValueDecomposition/explanation.md +++ /dev/null @@ -1,31 +0,0 @@ -# Singular Value Decomposition — Overview - -## What it is -The factorization `A = U Σ Vᵀ` of *any* matrix into two orthonormal rotations (`U`, `V`) and a -non-negative diagonal of **singular values** `Σ`. It is the most informative and most numerically -stable matrix factorization there is. - -## Why it matters (embedded) -The SVD is the swiss-army knife of linear algebra: the **pseudo-inverse** for robust least squares -and control allocation, **rank** and **condition-number** estimation to detect ill-posed problems, -and low-rank model reduction. It underpins calibration (`TotalLeastSquares`) and redundant-actuator -mapping. - -## How it works (intuition) -Two phases. First, alternating Householder reflectors squeeze `A` into a compact **bidiagonal** form -without changing its singular values. Second, implicit-shift QR sweeps apply tiny Givens rotations -that nibble the remaining off-diagonal to zero, leaving the singular values on the diagonal. The -accumulated rotations become `U` and `V`. Squaring the singular values recovers the eigenvalues of -`AᵀA`, tying the SVD directly to the symmetric eigenproblem. - -## Key parameters -- **Rank tolerance** — the threshold below which a singular value is treated as zero. -- **Shape (thin vs full)** — thin factors suffice for least squares and the pseudo-inverse. - -## Reference -G. Golub, W. Kahan, "Calculating the Singular Values and Pseudo-Inverse of a Matrix," -*SIAM J. Numer. Anal.*, 2(2), 1965; Golub & Reinsch, 1970. - -## See also -`QrDecomposition` and `JacobiEigenSolver` (its building blocks), `TotalLeastSquares`, -`LuDecomposition`. diff --git a/roadmap/solvers/SingularValueDecomposition/implementation.md b/roadmap/solvers/SingularValueDecomposition/implementation.md deleted file mode 100644 index 917f9182..00000000 --- a/roadmap/solvers/SingularValueDecomposition/implementation.md +++ /dev/null @@ -1,78 +0,0 @@ -# Singular Value Decomposition (Golub-Kahan) — Implementation Pseudocode - -> Roadmap ref: #43 (Tier 5) · 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 SingularValueDecomposition: - Matrix U # left singular vectors (thin) - Vector sigma # singular values, descending - Matrix V # right singular vectors - static_assert Rows >= Cols -``` - -## Interface - -``` -SingularValueDecomposition() -bool Decompose(const Matrix& A) # hot path -const Vector& SingularValues() const -Matrix PseudoInverse(T tol) const -std::size_t Rank(T tol) const -T ConditionNumber() const -Vector SolveLeastSquares(const Vector& b) const -``` - -## Algorithm (pseudocode) - -``` -function Decompose(A): # OPTIMIZE_FOR_SPEED (Golub-Kahan) - # Phase 1 — bidiagonalize with alternating Householder reflectors (reuse QrDecomposition) - (U, B, V) = HouseholderBidiagonalize(A) # B upper-bidiagonal - # Phase 2 — implicit-shift QR sweeps on B (Golub-Reinsch) - repeat: - for each split point: deflate if a superdiagonal ≈ 0 - apply Givens rotations to chase the off-diagonal down B # updates U, V - until B is diagonal within tolerance - sigma = |diag(B)|; fold signs into U - SortDescending(sigma, U, V) - return converged - -function PseudoInverse(tol): # A⁺ = V Σ⁺ Uᵀ - for i: sInv[i] = (sigma[i] > tol) ? 1/sigma[i] : 0 # threshold tiny σ - return V · diag(sInv) · Uᵀ - -function ConditionNumber(): return sigma_max / sigma_min -``` - -## Complexity & memory - -- Bidiagonalization dominates: `O(Rows·Cols² + Cols³)`; the QR sweeps add `O(Cols²)` per iteration. -- Memory: `O(Rows·Cols)` for `U` plus `O(Cols²)` for `V` — static, no heap. - -## Numerical / embedded notes - -- Reuse **Householder** reflectors (item 27) for the bidiagonalization and **Givens** rotations for - the sweeps; because `σ² = eig(AᵀA)`, results cross-check against `JacobiEigenSolver`. -- The **pseudo-inverse** thresholds tiny singular values to zero — this is what makes least squares - robust for rank-deficient or ill-conditioned `A` (regularization by truncation). -- `ConditionNumber = σ_max/σ_min` is the definitive conditioning metric; `Rank` counts `σ > tol`. -- For wide matrices (`Rows < Cols`) decompose `Aᵀ` and swap the roles of `U` and `V`. -- Foundational: powers `TotalLeastSquares`, pseudo-inverse control allocation, and model reduction. -- 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/SingularValueDecomposition.hpp` — `#pragma once` → - `#pragma GCC optimize("O3","fast-math")`, `OPTIMIZE_FOR_SPEED` on `Decompose`, and - `extern template class SingularValueDecomposition;` under `#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD`. -- Coverage: `numerical/solvers/SingularValueDecomposition.cpp` → - `template class SingularValueDecomposition;` and `` (the tested shapes). -- Test: `numerical/solvers/test/TestSingularValueDecomposition.cpp` -- Doc: `doc/solvers/SingularValueDecomposition.md` -- CMake: `.hpp` → `target_sources`; `.cpp` → `numerical_add_coverage_sources`; - `TestSingularValueDecomposition.cpp` → the `_test` target. -- Generic pattern: see `roadmap/DEPLOYMENT.md`. diff --git a/roadmap/solvers/SingularValueDecomposition/tests.md b/roadmap/solvers/SingularValueDecomposition/tests.md deleted file mode 100644 index 53e1b1d2..00000000 --- a/roadmap/solvers/SingularValueDecomposition/tests.md +++ /dev/null @@ -1,61 +0,0 @@ -# Singular Value Decomposition — Unit Test Plan (Pseudocode) - -> GoogleTest · `TEST_F` (`float`) · `StrictMock` only · no heap. - -## Fixture - -``` -class TestSingularValueDecomposition : public ::testing::Test: - SingularValueDecomposition svd # tall - SingularValueDecomposition svdSquare -# each case below is a TEST_F(TestSingularValueDecomposition, ) -``` - -## Test cases (Arrange / Act / Assert) - -``` -reconstructs_A: - Arrange: fixed 4×3 A - Act: Decompose(A) - Assert: U · diag(σ) · Vᵀ ≈ A - -singular_values_descending_and_nonnegative: - Assert: σ[0] >= σ[1] >= ... >= 0 - -U_and_V_orthonormal: - Assert: Uᵀ·U ≈ I and Vᵀ·V ≈ I - -diagonal_matrix_gives_absolute_diagonal: - Arrange: A = diag(3, −1, 2) - Assert: σ ≈ {3, 2, 1} - -sigma_squared_matches_eig_AtA: - Assert: σ² ≈ eigenvalues(AᵀA) (cross-check JacobiEigenSolver) - -pseudo_inverse_solves_least_squares: - Arrange: overdetermined system - Assert: A⁺·b matches the QR least-squares solution - -rank_detection_thresholds_small_sigma: - Arrange: rank-2 matrix living in a 3-column space - Assert: Rank(tol) == 2 - -condition_number_matches_known: - Arrange: A with prescribed σ_max, σ_min - Assert: ConditionNumber() ≈ σ_max/σ_min - -known_2x2_svd: - Assert: σ equals the documented Golub reference -``` - -## Reference vectors - -- `diag(3, −1, 2)` ⇒ `σ = {3, 2, 1}` (ordered, absolute value). -- `σ(A)² = eigenvalues(AᵀA)` — golden cross-check against `JacobiEigenSolver`. - -## Edge cases - -- Rank-deficient / zero singular value ⇒ pseudo-inverse truncates without divide-by-zero. -- Wide matrix (`Rows < Cols`) ⇒ decompose `Aᵀ`; assert consistent factors. -- All-zero matrix ⇒ all `σ = 0`, `Rank == 0`. -- Repeated singular values ⇒ vectors non-unique but the subspace is correct. From 67c8a3028ba9888da431cc433e7f31cd173e5b0b Mon Sep 17 00:00:00 2001 From: gfs Date: Sun, 2 Aug 2026 08:04:05 +0200 Subject: [PATCH 2/2] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- doc/solvers/README.md | 2 +- doc/solvers/SingularValueDecomposition.md | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/doc/solvers/README.md b/doc/solvers/README.md index 6d40ce06..484f7851 100644 --- a/doc/solvers/README.md +++ b/doc/solvers/README.md @@ -18,4 +18,4 @@ Numerical solvers for linear systems, polynomial roots, and matrix equations. | [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 | | [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 | +| [Jacobi Eigenvalue Solver](JacobiEigenSolver.md) | Cyclic Jacobi rotations for the full symmetric eigenvalue/eigenvector problem | diff --git a/doc/solvers/SingularValueDecomposition.md b/doc/solvers/SingularValueDecomposition.md index 25b88d23..aeb20c47 100644 --- a/doc/solvers/SingularValueDecomposition.md +++ b/doc/solvers/SingularValueDecomposition.md @@ -50,12 +50,12 @@ 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 | +| 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.