From 7022ec3215e37fc859e240d69a46c6b060f727ca Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Wed, 29 Jul 2026 15:51:21 +0000 Subject: [PATCH 1/2] add spectral radius --- README.md | 2 +- ROADMAP.md | 1 - doc/solvers/README.md | 1 + doc/solvers/SpectralRadius.md | 91 +++++++++++++++++ numerical/solvers/CMakeLists.txt | 2 + numerical/solvers/SpectralRadius.cpp | 6 ++ numerical/solvers/SpectralRadius.hpp | 98 +++++++++++++++++++ numerical/solvers/test/CMakeLists.txt | 1 + numerical/solvers/test/TestSpectralRadius.cpp | 79 +++++++++++++++ 9 files changed, 279 insertions(+), 2 deletions(-) create mode 100644 doc/solvers/SpectralRadius.md create mode 100644 numerical/solvers/SpectralRadius.cpp create mode 100644 numerical/solvers/SpectralRadius.hpp create mode 100644 numerical/solvers/test/TestSpectralRadius.cpp diff --git a/README.md b/README.md index 28eb51c2..341141ee 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 | -| [Solvers](doc/solvers/README.md) | Gaussian Elimination, Levinson-Durbin, Durand-Kerner, Cholesky, DARE, Runge-Kutta ODE Integrators (RK4 + Dormand-Prince) | +| [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 | | [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 554ed89b..c30ea371 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -51,7 +51,6 @@ Difficulty legend: | 48 | Decibel & magnitude-response helpers | `analysis` | ★☆☆☆☆ | | 49 | Step / transient-response metrics | `math` | ★★☆☆☆ | | 50 | Matrix norms & condition number | `math` | ★★★☆☆ | -| 51 | Spectral radius / discrete stability margin | `math` | ★★★☆☆ | | 52 | Estimator consistency metrics (NEES / NIS) | `estimators` | ★★★☆☆ | Items 48–52 are the **evaluation & metrics primitives** — reusable quantities the per-family diff --git a/doc/solvers/README.md b/doc/solvers/README.md index 8dd78088..a703566c 100644 --- a/doc/solvers/README.md +++ b/doc/solvers/README.md @@ -12,3 +12,4 @@ Numerical solvers for linear systems, polynomial roots, and matrix equations. | [Durand-Kerner](DurandKerner.md) | Simultaneous iterative root-finder for polynomials | | [Discrete Algebraic Riccati Equation](DiscreteAlgebraicRiccatiEquation.md) | Iterative solver for the DARE arising in LQR and Kalman filter design | | [Runge-Kutta ODE Integrators](RungeKuttaIntegrators.md) | Fixed-step RK4 and adaptive Dormand-Prince RK45 for ODE integration | +| [Spectral Radius](SpectralRadius.md) | Dominant eigenvalue magnitude and discrete-time Schur stability margin | diff --git a/doc/solvers/SpectralRadius.md b/doc/solvers/SpectralRadius.md new file mode 100644 index 00000000..2ac6ecc3 --- /dev/null +++ b/doc/solvers/SpectralRadius.md @@ -0,0 +1,91 @@ +# Spectral Radius & Discrete Stability Margin + +## Overview & Motivation + +The spectral radius of a square matrix is the largest magnitude among all its eigenvalues. In discrete-time dynamical systems, the spectral radius of the state-transition matrix is the single most important scalar quantity: if it is strictly less than one, every state trajectory converges to zero; if it equals or exceeds one, at least one mode grows without bound or fails to decay. Rather than reporting only a binary stable/unstable flag, exposing the spectral radius and the complementary stability margin (one minus the spectral radius) gives the designer a quantitative handle on how much gain or perturbation the system can tolerate before losing stability. + +## Mathematical Theory + +### Characteristic Polynomial via Faddeev–LeVerrier + +Given an $n \times n$ matrix $A$, the characteristic polynomial is + +$$p(\lambda) = \det(\lambda I - A) = \lambda^n + c_{n-1}\lambda^{n-1} + \cdots + c_1\lambda + c_0.$$ + +The Faddeev–LeVerrier recurrence computes the coefficients $c_k$ without determinant expansion. Define auxiliary matrices $M_k$ by + +$$M_1 = A, \quad c_{n-1} = -\operatorname{tr}(M_1),$$ +$$M_k = A\bigl(M_{k-1} + c_{n-k+1}I\bigr), \quad c_{n-k} = -\frac{1}{k}\operatorname{tr}(M_k), \quad k = 2, \ldots, n.$$ + +This recurrence requires only matrix-matrix multiplication, scalar multiplication of a matrix by the identity, and the trace, making it suitable for static (stack-allocated) computation without a heap. + +### Spectral Radius + +Given the full set of eigenvalues $\{\lambda_i\}_{i=1}^{n}$ (roots of the characteristic polynomial), + +$$\rho(A) = \max_{i} |\lambda_i|.$$ + +The roots are found numerically with the Durand–Kerner simultaneous root-finding iteration, which converges for arbitrary polynomials and requires only a fixed-size array of complex iterates. + +### Schur Stability + +A matrix $A$ is Schur-stable (the discrete-time analogue of Hurwitz stability) if and only if $\rho(A) < 1$. Every eigenvalue then lies strictly inside the complex unit disk, guaranteeing exponential decay of all modes. + +### Stability Margin + +$$\delta(A) = 1 - \rho(A).$$ + +A positive margin indicates Schur stability; $\delta = 0$ corresponds to marginal stability (an eigenvalue on the unit circle); $\delta < 0$ signals instability. The margin quantifies how far $\rho(A)$ is from the stability boundary. + +## Complexity Analysis + +| Case | Time | Space | Notes | +|---------|------------------|-----------------|-----------------------------------------------| +| Best | $O(n^3)$ | $O(n^2)$ | Faddeev–LeVerrier dominates; two $n\times n$ temporaries on the stack | +| Average | $O(n^3 + n \cdot I)$ | $O(n^2 + n)$ | $I$ Durand–Kerner iterations, typically $\ll 200$ | +| Worst | $O(n^3 + n \cdot I_{\max})$ | $O(n^2 + n)$ | Near-repeated roots slow DK convergence | + +The Faddeev–LeVerrier step executes $n$ matrix multiplications each costing $O(n^3)$, giving $O(n^4)$ total; for the small dimensions typical of embedded control matrices ($n \le 6$) this is negligible. + +## Step-by-Step Walkthrough + +Consider $A = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$ (90° rotation, $\rho = 1$). + +1. **k = 1**: $M_1 = A \cdot 0 + I \cdot c_0$. With $c_0 = 1$: $M_1 = I$. Then $c_1 = -\tfrac{1}{1}\operatorname{tr}(A \cdot M_1) = -\operatorname{tr}(A) = 0$. +2. **k = 2**: $M_2 = A(M_1 + c_1 I) = A \cdot I = A$. Then $c_2 = -\tfrac{1}{2}\operatorname{tr}(A \cdot A) = -\tfrac{1}{2}\operatorname{tr}(-I) = 1$. + +Characteristic polynomial: $\lambda^2 + 0\lambda + 1 = \lambda^2 + 1$. Roots: $\pm i$, magnitudes both 1. $\rho = 1$. + +## Pitfalls & Edge Cases + +- **Near-repeated roots**: Durand–Kerner converges slowly when two or more eigenvalues are nearly identical. The stability conclusion (inside/outside unit disk) is unaffected, but the raw $\rho$ value may have larger numerical error. +- **Identity approximation in Matrix::Identity()**: This implementation builds its own exact identity to avoid reliance on library internals. +- **Ill-conditioned characteristic polynomials**: For large $n$ the polynomial coefficients can have widely different magnitudes, degrading root-finder precision. The algorithm is intended for small matrices ($n \le 6$) typical of discrete-time state-space models. +- **Size = 1 case**: Handled directly as $\rho = |a_{11}|$, bypassing polynomial computation. +- **Zero matrix**: $\rho = 0$, $\delta = 1$ (maximally stable). Handled correctly by the recurrence. + +## Variants & Generalizations + +- **Power iteration** computes only the dominant eigenvalue magnitude without the full spectrum; cheaper for large $n$ but provides no stability margin without additional eigenpair deflation. +- **QR algorithm** (Jacobi, Hessenberg QR) finds all eigenvalues more robustly for large or symmetric matrices; too expensive and heap-dependent for small embedded targets. +- **Continuous-time analogue**: Replace the unit-disk criterion with the left half-plane criterion; Routh–Hurwitz or Lyapunov-based approaches apply. + +## Applications + +- Verifying discrete-time state-feedback or observer designs remain within the unit disk after coefficient quantisation. +- Checking closed-loop stability after LQR or pole-placement synthesis before deployment. +- Computing gain/phase margins indirectly by perturbing the loop gain and re-evaluating $\rho$. +- Monitoring on-line whether a time-varying system's frozen eigenvalues have drifted outside the unit disk. + +## Connections to Other Algorithms + +- Uses **Durand–Kerner** (polynomial root-finder) for the eigenvalue computation stage. +- The characteristic polynomial recurrence is mathematically related to **Leverrier's method** for matrix inverses and adjoints. +- The Schur-stability criterion is the discrete-time counterpart of the Hurwitz criterion used in continuous-time control analysis. +- The stability margin scalar is a prerequisite for **TESTING.md** evaluation metrics that require a quantitative stability distance. + +## References & Further Reading + +- Faddeev, D. K. & LeVerrier, U. J. F. — classical matrix identity; see G. H. Golub & C. F. Van Loan, *Matrix Computations*, 4th ed., Section 7.3. +- Durand, E. (1960) and Kerner, I. O. (1966) — simultaneous root-finding; see *Numerical Recipes in C*, 3rd ed., Section 9.5. +- Schur stability: K. J. Åström & B. Wittenmark, *Computer-Controlled Systems*, 3rd ed., Chapter 3. diff --git a/numerical/solvers/CMakeLists.txt b/numerical/solvers/CMakeLists.txt index 5418a4f7..c092e6f8 100644 --- a/numerical/solvers/CMakeLists.txt +++ b/numerical/solvers/CMakeLists.txt @@ -20,6 +20,7 @@ target_sources(numerical.solver PRIVATE OdeSystem.hpp RungeKuttaIntegrators.hpp Solver.hpp + SpectralRadius.hpp ) numerical_add_coverage_sources(numerical.solver @@ -27,6 +28,7 @@ numerical_add_coverage_sources(numerical.solver DurandKerner.cpp GaussianElimination.cpp RungeKuttaIntegrators.cpp + SpectralRadius.cpp ) add_subdirectory(test) diff --git a/numerical/solvers/SpectralRadius.cpp b/numerical/solvers/SpectralRadius.cpp new file mode 100644 index 00000000..8027108d --- /dev/null +++ b/numerical/solvers/SpectralRadius.cpp @@ -0,0 +1,6 @@ +#include "numerical/solvers/SpectralRadius.hpp" + +namespace solvers +{ + template class SpectralRadius; +} diff --git a/numerical/solvers/SpectralRadius.hpp b/numerical/solvers/SpectralRadius.hpp new file mode 100644 index 00000000..a8e459cd --- /dev/null +++ b/numerical/solvers/SpectralRadius.hpp @@ -0,0 +1,98 @@ +#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/solvers/DurandKerner.hpp" +#include +#include +#include +#include +#include + +namespace solvers +{ + template + class SpectralRadius + { + static_assert(std::is_floating_point_v, + "SpectralRadius supports floating-point types only"); + + public: + OPTIMIZE_FOR_SPEED T Compute(const math::SquareMatrix& matrix) const; + OPTIMIZE_FOR_SPEED bool IsSchurStable(const math::SquareMatrix& matrix) const; + OPTIMIZE_FOR_SPEED T StabilityMargin(const math::SquareMatrix& matrix) const; + + private: + static OPTIMIZE_FOR_SPEED math::SquareMatrix MakeIdentity(); + static OPTIMIZE_FOR_SPEED std::array CharacteristicPolynomial(const math::SquareMatrix& matrix); + }; + + template + math::SquareMatrix SpectralRadius::MakeIdentity() + { + math::SquareMatrix result{}; + for (std::size_t i = 0; i < Size; ++i) + result.at(i, i) = T(1); + return result; + } + + template + std::array SpectralRadius::CharacteristicPolynomial(const math::SquareMatrix& matrix) + { + std::array coefficients{}; + coefficients[0] = T(1); + + math::SquareMatrix M{}; + const math::SquareMatrix identity = MakeIdentity(); + + for (std::size_t k = 1; k <= Size; ++k) + { + M = matrix * M + identity * coefficients[k - 1]; + coefficients[k] = -(matrix * M).Trace() / static_cast(k); + } + + return coefficients; + } + + template + T SpectralRadius::Compute(const math::SquareMatrix& matrix) const + { + if constexpr (Size == 1) + return std::abs(matrix.at(0, 0)); + + const std::array poly = CharacteristicPolynomial(matrix); + + DurandKerner solver{}; + auto roots = solver.Solve(std::span{ poly.data(), poly.size() }); + + T rho{}; + for (std::size_t i = 0; i < roots.size(); ++i) + { + T magnitude = std::abs(roots[i]); + if (magnitude > rho) + rho = magnitude; + } + + return rho; + } + + template + bool SpectralRadius::IsSchurStable(const math::SquareMatrix& matrix) const + { + return Compute(matrix) < T(1); + } + + template + T SpectralRadius::StabilityMargin(const math::SquareMatrix& matrix) const + { + return T(1) - Compute(matrix); + } + +#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD + extern template class SpectralRadius; +#endif +} diff --git a/numerical/solvers/test/CMakeLists.txt b/numerical/solvers/test/CMakeLists.txt index 61fd5772..f97b8d42 100644 --- a/numerical/solvers/test/CMakeLists.txt +++ b/numerical/solvers/test/CMakeLists.txt @@ -13,4 +13,5 @@ target_sources(numerical.solvers_test PRIVATE TestGaussianElimination.cpp TestLevinsonDurbin.cpp TestRungeKuttaIntegrators.cpp + TestSpectralRadius.cpp ) diff --git a/numerical/solvers/test/TestSpectralRadius.cpp b/numerical/solvers/test/TestSpectralRadius.cpp new file mode 100644 index 00000000..904905af --- /dev/null +++ b/numerical/solvers/test/TestSpectralRadius.cpp @@ -0,0 +1,79 @@ +#include "numerical/math/Tolerance.hpp" +#include "numerical/solvers/SpectralRadius.hpp" +#include +#include + +namespace +{ + class TestSpectralRadius : public ::testing::Test + { + protected: + solvers::SpectralRadius solver3; + solvers::SpectralRadius solver2; + solvers::SpectralRadius solver1; + }; +} + +TEST_F(TestSpectralRadius, diagonal_matrix_rho_equals_largest_magnitude_diagonal) +{ + math::SquareMatrix A{ { 1.0f, 0.0f, 0.0f }, { 0.0f, 2.0f, 0.0f }, { 0.0f, 0.0f, 3.0f } }; + + float rho = solver3.Compute(A); + + EXPECT_NEAR(rho, 3.0f, 1e-2f); +} + +TEST_F(TestSpectralRadius, known_matrix_rho_matches_analytic_dominant_eigenvalue) +{ + math::SquareMatrix A{ { 3.0f, 1.0f, 0.0f }, { 0.0f, 2.0f, 0.0f }, { 0.0f, 0.0f, 1.0f } }; + + float rho = solver3.Compute(A); + + EXPECT_NEAR(rho, 3.0f, 1e-2f); +} + +TEST_F(TestSpectralRadius, complex_conjugate_pair_rho_equals_modulus) +{ + float r = 0.5f; + float theta = std::numbers::pi_v / 4.0f; + math::SquareMatrix A{ + { r * std::cos(theta), -r * std::sin(theta) }, + { r * std::sin(theta), r * std::cos(theta) } + }; + + float rho = solver2.Compute(A); + + EXPECT_NEAR(rho, r, 1e-2f); +} + +TEST_F(TestSpectralRadius, schur_stable_true_when_all_eigenvalues_inside_unit_disk) +{ + math::SquareMatrix A{ { 0.5f, 0.0f, 0.0f }, { 0.0f, 0.3f, 0.0f }, { 0.0f, 0.0f, 0.1f } }; + + EXPECT_TRUE(solver3.IsSchurStable(A)); +} + +TEST_F(TestSpectralRadius, schur_stable_false_when_eigenvalue_outside_unit_disk) +{ + math::SquareMatrix A{ { 0.5f, 0.0f, 0.0f }, { 0.0f, 1.1f, 0.0f }, { 0.0f, 0.0f, 0.1f } }; + + EXPECT_FALSE(solver3.IsSchurStable(A)); +} + +TEST_F(TestSpectralRadius, stability_margin_equals_one_minus_rho) +{ + math::SquareMatrix A{ { 0.5f, 0.0f, 0.0f }, { 0.0f, 0.3f, 0.0f }, { 0.0f, 0.0f, 0.1f } }; + + float margin = solver3.StabilityMargin(A); + + EXPECT_NEAR(margin, 0.5f, 1e-2f); +} + +TEST_F(TestSpectralRadius, scalar_size_one_rho_equals_absolute_value_of_entry) +{ + math::SquareMatrix A{ { -3.0f } }; + + float rho = solver1.Compute(A); + + EXPECT_NEAR(rho, 3.0f, math::Tolerance()); +} From c8d07e9eab86dce4e2083e7f1bac8fae9524c97d Mon Sep 17 00:00:00 2001 From: gfs Date: Wed, 29 Jul 2026 19:40:24 +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/SpectralRadius.md | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/solvers/README.md b/doc/solvers/README.md index a703566c..18d8b2d6 100644 --- a/doc/solvers/README.md +++ b/doc/solvers/README.md @@ -12,4 +12,4 @@ Numerical solvers for linear systems, polynomial roots, and matrix equations. | [Durand-Kerner](DurandKerner.md) | Simultaneous iterative root-finder for polynomials | | [Discrete Algebraic Riccati Equation](DiscreteAlgebraicRiccatiEquation.md) | Iterative solver for the DARE arising in LQR and Kalman filter design | | [Runge-Kutta ODE Integrators](RungeKuttaIntegrators.md) | Fixed-step RK4 and adaptive Dormand-Prince RK45 for ODE integration | -| [Spectral Radius](SpectralRadius.md) | Dominant eigenvalue magnitude and discrete-time Schur stability margin | +| [Spectral Radius](SpectralRadius.md) | Dominant eigenvalue magnitude and discrete-time Schur stability margin | diff --git a/doc/solvers/SpectralRadius.md b/doc/solvers/SpectralRadius.md index 2ac6ecc3..5420a9f0 100644 --- a/doc/solvers/SpectralRadius.md +++ b/doc/solvers/SpectralRadius.md @@ -39,11 +39,11 @@ A positive margin indicates Schur stability; $\delta = 0$ corresponds to margina ## Complexity Analysis -| Case | Time | Space | Notes | -|---------|------------------|-----------------|-----------------------------------------------| -| Best | $O(n^3)$ | $O(n^2)$ | Faddeev–LeVerrier dominates; two $n\times n$ temporaries on the stack | -| Average | $O(n^3 + n \cdot I)$ | $O(n^2 + n)$ | $I$ Durand–Kerner iterations, typically $\ll 200$ | -| Worst | $O(n^3 + n \cdot I_{\max})$ | $O(n^2 + n)$ | Near-repeated roots slow DK convergence | +| Case | Time | Space | Notes | +|---------|-----------------------------|--------------|-----------------------------------------------------------------------| +| Best | $O(n^3)$ | $O(n^2)$ | Faddeev–LeVerrier dominates; two $n\times n$ temporaries on the stack | +| Average | $O(n^3 + n \cdot I)$ | $O(n^2 + n)$ | $I$ Durand–Kerner iterations, typically $\ll 200$ | +| Worst | $O(n^3 + n \cdot I_{\max})$ | $O(n^2 + n)$ | Near-repeated roots slow DK convergence | The Faddeev–LeVerrier step executes $n$ matrix multiplications each costing $O(n^3)$, giving $O(n^4)$ total; for the small dimensions typical of embedded control matrices ($n \le 6$) this is negligible.