feat: add spectral radius

README.md

@@ -26,8 +26,8 @@ Refer to the documentation to quickly integrate and utilize the library's signal
 | [Neural Network](doc/neural_network/README.md)                     | Layers, activations, losses, model                                   |
 | [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), Condition Number |
+| [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 |
 | [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

@@ -48,7 +48,6 @@ Difficulty legend:
 | 45 | IIR filter design (Butterworth/Chebyshev + bilinear) | `filters/passive`         | ★★★★★      |
 | 46 | H∞ state-feedback control                            | `robust_control` (new)    | ★★★★★      |
 | 47 | Model Reference Adaptive Control (MRAC)              | `nonlinear_control` (new) | ★★★★★      |
-| 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

doc/solvers/README.md

@@ -12,4 +12,5 @@ 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 |
 | [Condition Number](ConditionNumber.md)                                     | 1-norm condition number estimate via column-wise inverse solve         |

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.

numerical/solvers/CMakeLists.txt

@@ -21,13 +21,15 @@ target_sources(numerical.solver PRIVATE
     OdeSystem.hpp
     RungeKuttaIntegrators.hpp
     Solver.hpp
+    SpectralRadius.hpp
 )
 
 numerical_add_coverage_sources(numerical.solver
     DiscreteAlgebraicRiccatiEquation.cpp
     DurandKerner.cpp
     GaussianElimination.cpp
     RungeKuttaIntegrators.cpp
+    SpectralRadius.cpp
 )
 
 add_subdirectory(test)

numerical/solvers/SpectralRadius.cpp

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

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 <array>
+#include <cmath>
+#include <cstddef>
+#include <span>
+#include <type_traits>
+
+namespace solvers
+{
+    template<typename T, std::size_t Size>
+    class SpectralRadius
+    {
+        static_assert(std::is_floating_point_v<T>,
+            "SpectralRadius supports floating-point types only");
+
+    public:
+        OPTIMIZE_FOR_SPEED T Compute(const math::SquareMatrix<T, Size>& matrix) const;
+        OPTIMIZE_FOR_SPEED bool IsSchurStable(const math::SquareMatrix<T, Size>& matrix) const;
+        OPTIMIZE_FOR_SPEED T StabilityMargin(const math::SquareMatrix<T, Size>& matrix) const;
+
+    private:
+        static OPTIMIZE_FOR_SPEED math::SquareMatrix<T, Size> MakeIdentity();
+        static OPTIMIZE_FOR_SPEED std::array<T, Size + 1> CharacteristicPolynomial(const math::SquareMatrix<T, Size>& matrix);
+    };
+
+    template<typename T, std::size_t Size>
+    math::SquareMatrix<T, Size> SpectralRadius<T, Size>::MakeIdentity()
+    {
+        math::SquareMatrix<T, Size> result{};
+        for (std::size_t i = 0; i < Size; ++i)
+            result.at(i, i) = T(1);
+        return result;
+    }
+
+    template<typename T, std::size_t Size>
+    std::array<T, Size + 1> SpectralRadius<T, Size>::CharacteristicPolynomial(const math::SquareMatrix<T, Size>& matrix)
+    {
+        std::array<T, Size + 1> coefficients{};
+        coefficients[0] = T(1);
+
+        math::SquareMatrix<T, Size> M{};
+        const math::SquareMatrix<T, Size> 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<T>(k);
+        }
+
+        return coefficients;
+    }
+
+    template<typename T, std::size_t Size>
+    T SpectralRadius<T, Size>::Compute(const math::SquareMatrix<T, Size>& matrix) const
+    {
+        if constexpr (Size == 1)
+            return std::abs(matrix.at(0, 0));
+
+        const std::array<T, Size + 1> poly = CharacteristicPolynomial(matrix);
+
+        DurandKerner<T, Size> solver{};
+        auto roots = solver.Solve(std::span<const T>{ 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<typename T, std::size_t Size>
+    bool SpectralRadius<T, Size>::IsSchurStable(const math::SquareMatrix<T, Size>& matrix) const
+    {
+        return Compute(matrix) < T(1);
+    }
+
+    template<typename T, std::size_t Size>
+    T SpectralRadius<T, Size>::StabilityMargin(const math::SquareMatrix<T, Size>& matrix) const
+    {
+        return T(1) - Compute(matrix);
+    }
+
+#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
+    extern template class SpectralRadius<float, 3>;
+#endif
+}

numerical/solvers/test/CMakeLists.txt

@@ -14,4 +14,5 @@ target_sources(numerical.solvers_test PRIVATE
     TestGaussianElimination.cpp
     TestLevinsonDurbin.cpp
     TestRungeKuttaIntegrators.cpp
+    TestSpectralRadius.cpp
 )

numerical/solvers/test/TestSpectralRadius.cpp

@@ -0,0 +1,79 @@
+#include "numerical/math/Tolerance.hpp"
+#include "numerical/solvers/SpectralRadius.hpp"
+#include <gtest/gtest.h>
+#include <numbers>
+
+namespace
+{
+    class TestSpectralRadius : public ::testing::Test
+    {
+    protected:
+        solvers::SpectralRadius<float, 3> solver3;
+        solvers::SpectralRadius<float, 2> solver2;
+        solvers::SpectralRadius<float, 1> solver1;
+    };
+}
+
+TEST_F(TestSpectralRadius, diagonal_matrix_rho_equals_largest_magnitude_diagonal)
+{
+    math::SquareMatrix<float, 3> 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<float, 3> 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<float> / 4.0f;
+    math::SquareMatrix<float, 2> 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<float, 3> 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<float, 3> 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<float, 3> 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<float, 1> A{ { -3.0f } };
+
+    float rho = solver1.Compute(A);
+
+    EXPECT_NEAR(rho, 3.0f, math::Tolerance<float>());
+}