diff --git a/README.md b/README.md index fcb9f2f7..0c9c6e5e 100644 --- a/README.md +++ b/README.md @@ -26,7 +26,7 @@ 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, MatrixExponential | +| [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 | | [Performance Optimization](doc/performance-optimization/README.md) | Compiler optimizations, SIMD | diff --git a/ROADMAP.md b/ROADMAP.md index 3a7417c5..516305ad 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -28,7 +28,6 @@ Difficulty legend: | # | Component | Target module | Difficulty | |----|------------------------------------------------------|---------------------------|------------| | 30 | Continuous → discrete conversion (`c2d`) | `math` | ★★★★☆ | -| 31 | Lyapunov / Sylvester equation solvers | `solvers` | ★★★★☆ | | 32 | Transfer-function ↔ state-space conversion | `control_analysis` | ★★★★☆ | | 33 | Madgwick / Mahony AHRS | `filters/active` | ★★★★☆ | | 34 | Sliding Mode Control (SMC) | `robust_control` (new) | ★★★★☆ | diff --git a/doc/solvers/LyapunovSylvester.md b/doc/solvers/LyapunovSylvester.md new file mode 100644 index 00000000..86f1121a --- /dev/null +++ b/doc/solvers/LyapunovSylvester.md @@ -0,0 +1,139 @@ +# Lyapunov and Sylvester Equation Solvers + +## Overview & Motivation + +The Sylvester equation $AX + XB = C$ and its symmetric specialisations — the continuous Lyapunov +$AX + XA^\top = -Q$ and the discrete Lyapunov $AXA^\top - X = -Q$ — are the central matrix +equations of stability analysis, gramian computation, and robust-control synthesis. A +positive-definite Lyapunov solution is a formal certificate that the system $\dot{x} = Ax$ is +asymptotically stable; the same equations produce the **controllability and observability gramians** +that measure how effectively a linear system can be driven and observed. + +## Mathematical Theory + +### Sylvester Equation + +Given matrices $A \in \mathbb{R}^{N \times N}$, $B \in \mathbb{R}^{M \times M}$, and +$C \in \mathbb{R}^{N \times M}$, the Sylvester equation is + +$$AX + XB = C,$$ + +with unknown $X \in \mathbb{R}^{N \times M}$. The equation has a unique solution if and only if +$A$ and $-B$ share no eigenvalue, equivalently $\lambda_i(A) + \lambda_j(B) \neq 0$ for all +$i, j$. + +### Kronecker Vectorisation + +Stacking $X$ column-by-column into $\operatorname{vec}(X) \in \mathbb{R}^{NM}$ and applying +the identities $\operatorname{vec}(AXI) = (I_M \otimes A)\operatorname{vec}(X)$ and +$\operatorname{vec}(IXB) = (B^\top \otimes I_N)\operatorname{vec}(X)$ converts the matrix +equation into the standard linear system + +$$\underbrace{(I_M \otimes A + B^\top \otimes I_N)}_{K \in \mathbb{R}^{NM \times NM}} +\operatorname{vec}(X) = \operatorname{vec}(C).$$ + +The $(iN+k, jN+l)$ entry of $K$ is $\delta_{ij} A_{kl} + B_{ji} \delta_{kl}$. The system +is solved by Gaussian elimination with partial pivoting; a near-zero pivot signals the +unsolvable (degenerate eigenvalue) case. + +### Continuous Lyapunov Equation + +Setting $B = A^\top$ and $C = -Q$ reduces the Sylvester equation to the continuous Lyapunov form + +$$AX + XA^\top = -Q.$$ + +When $A$ is Hurwitz (all eigenvalues in the open left half-plane) and $Q$ is symmetric +positive semi-definite, the unique solution $X$ is symmetric positive definite. + +### Discrete Lyapunov Equation + +The discrete counterpart + +$$AXA^\top - X = -Q$$ + +is solved by the Kronecker system $(A \otimes A - I_{N^2})\operatorname{vec}(X) = -\operatorname{vec}(Q)$, +with the $(iN+k, jN+l)$ entry of $A \otimes A - I$ equal to $A_{ij}A_{kl} - \delta_{ij}\delta_{kl}$. +When $A$ is Schur-stable ($\rho(A) < 1$) and $Q \succeq 0$, the unique solution equals the +controllability gramian $X = \sum_{k=0}^{\infty} A^k Q (A^\top)^k$. + +## Complexity Analysis + +| Case | Time | Space | Notes | +|-------------------------------|----------------|-----------------------------|---------------------------------------------------------------| +| $N{=}M{=}2$ | $O(64)$ flops | $O(NM)^2 = 16$ floats for K | 4×4 linear system | +| $N{=}M{=}3$ | $O(729)$ flops | 81 floats for K | 9×9 linear system | +| General $N{=}M$ | $O(N^6)$ flops | $O(N^4)$ words | Gaussian elimination on the $N^2 \times N^2$ Kronecker system | +| Alternative (Bartels–Stewart) | $O(N^3)$ | $O(N^2)$ | Requires Schur decomposition; preferred for $N \geq 6$ | + +For the embedded sizes ($N \leq 5$) targeted by this library the Kronecker approach is negligibly +more expensive than Bartels–Stewart and avoids the full Schur machinery. + +## Step-by-Step Walkthrough + +Solve $AX + XB = C$ for the $2 \times 2$ case +$A = \bigl[\begin{smallmatrix}-3&1\\0&-2\end{smallmatrix}\bigr]$, +$B = \bigl[\begin{smallmatrix}1&2\\0&1\end{smallmatrix}\bigr]$, +$C = I_2$. + +**Build $K$** (the $4 \times 4$ Kronecker matrix): + +$$K = \begin{pmatrix}-2&1&0&0\\0&-1&0&0\\2&0&-2&1\\0&2&0&-1\end{pmatrix}.$$ + +**Build $\operatorname{vec}(C) = (1, 0, 0, 1)^\top$.** + +**Gaussian elimination** with partial pivoting produces $X_{00} = -0.5$, $X_{10} = 0$, +$X_{01} = -1$, $X_{11} = -1$. + +**Verify:** $AX + XB = I_2$. ✓ + +## Pitfalls & Edge Cases + +- **Unsolvable case:** if $\lambda_i(A) + \lambda_j(B) = 0$ (continuous) or + $\lambda_i(A)\lambda_j(B) = 1$ (discrete) for any pair, $K$ is singular and no unique solution + exists. The solver returns `false` and leaves the solution unchanged. +- **Ill-conditioning:** near-degenerate eigenvalue pairs yield large solution norms and reduced + accuracy. Monitor the residual $\|AX + XB - C\|$ after solving. +- **Symmetry:** the continuous Lyapunov solution is guaranteed symmetric only when $Q$ is + symmetric; numerical symmetry should be enforced by the caller if required downstream. + +## Variants & Generalizations + +- **Bartels–Stewart algorithm:** reduces to quasi-upper-triangular (Schur) form via the double-shift + QR algorithm, then back-substitutes column by column in $O(N^3)$. Preferred for $N \geq 6$. +- **Smith iteration:** fixed-point iteration $X_{k+1} = AX_k A^\top + Q$ for the discrete + Lyapunov equation; simple but slow for $\rho(A)$ close to 1. +- **Cholesky factor updates:** propagate a Cholesky factor of $X$ for guaranteed positive-definiteness + in square-root filtering (see `SquareRootKalmanFilter`). + +## Applications + +- **LQR / LQG design:** the continuous Lyapunov equation $AX + XA^\top + Q = 0$ yields the + closed-loop covariance, and the DARE solution can be cross-validated via its Lyapunov residual. +- **Kalman filter (observability Gramian):** $AW + WA^\top + BB^\top = 0$ gives the + controllability Gramian used to evaluate reachability in balanced truncation. +- **Pole-placement observers (Luenberger):** the Sylvester equation $AX - XF = LC$ arises when + placing observer poles while decoupling the error dynamics. +- **H₂ norm computation:** $\|G\|_2^2 = \operatorname{tr}(B^\top X B)$ where $X$ satisfies the + observability Gramian Lyapunov equation. + +## Connections to Other Algorithms + +- **LU Decomposition:** the Kronecker-vectorisation approach converts any Sylvester equation into + a dense linear system solved by the `LuDecomposition` solver. +- **QR Decomposition:** the Bartels–Stewart and Hammarling algorithms use QR / Schur + decomposition to reduce the problem to quasi-triangular form before back-substitution. +- **DARE:** the discrete algebraic Riccati equation can be warm-started or verified using + one step of the discrete Lyapunov equation. +- **Spectral Radius:** solvability of the continuous equation requires $\lambda_i(A) + \lambda_j(A) \neq 0$; + checking the spectral radius is a fast necessary condition. + +## References & Further Reading + +- Bartels, R. H. & Stewart, G. W. (1972). "Solution of the matrix equation AX + XB = C." + *Communications of the ACM*, 15(9), 820–826. +- Golub, G. H., Nash, S. & Van Loan, C. (1979). "A Hessenberg-Schur method for the problem + AX + XB = C." *IEEE Transactions on Automatic Control*, 24(6), 909–913. +- Horn, R. A. & Johnson, C. R. (1994). *Topics in Matrix Analysis*. Cambridge University Press. + (Chapter 4: Kronecker products and the vec operator.) +- Hammarling, S. J. (1982). "Numerical solution of the stable, non-negative definite Lyapunov + equation." *IMA Journal of Numerical Analysis*, 2(3), 303–323. diff --git a/doc/solvers/README.md b/doc/solvers/README.md index 1eb988ee..e2ef999c 100644 --- a/doc/solvers/README.md +++ b/doc/solvers/README.md @@ -16,3 +16,4 @@ Numerical solvers for linear systems, polynomial roots, and matrix equations. | [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 | +| [Lyapunov / Sylvester Solvers](LyapunovSylvester.md) | Sylvester AX+XB=C and continuous/discrete Lyapunov solvers via Kronecker vectorisation | diff --git a/numerical/solvers/CMakeLists.txt b/numerical/solvers/CMakeLists.txt index faac6736..6abc38d9 100644 --- a/numerical/solvers/CMakeLists.txt +++ b/numerical/solvers/CMakeLists.txt @@ -19,6 +19,7 @@ target_sources(numerical.solver PRIVATE GaussianElimination.hpp LevinsonDurbin.hpp LuDecomposition.hpp + LyapunovSylvester.hpp OdeSystem.hpp QrDecomposition.hpp RungeKuttaIntegrators.hpp @@ -31,6 +32,7 @@ numerical_add_coverage_sources(numerical.solver DurandKerner.cpp GaussianElimination.cpp LuDecomposition.cpp + LyapunovSylvester.cpp QrDecomposition.cpp RungeKuttaIntegrators.cpp SpectralRadius.cpp diff --git a/numerical/solvers/LyapunovSylvester.cpp b/numerical/solvers/LyapunovSylvester.cpp new file mode 100644 index 00000000..ea0ce56f --- /dev/null +++ b/numerical/solvers/LyapunovSylvester.cpp @@ -0,0 +1,7 @@ +#include "numerical/solvers/LyapunovSylvester.hpp" + +namespace solvers +{ + template class SylvesterSolver; + template class SylvesterSolver; +} diff --git a/numerical/solvers/LyapunovSylvester.hpp b/numerical/solvers/LyapunovSylvester.hpp new file mode 100644 index 00000000..5f2cdeed --- /dev/null +++ b/numerical/solvers/LyapunovSylvester.hpp @@ -0,0 +1,148 @@ +#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/LuDecomposition.hpp" +#include +#include + +namespace solvers +{ + template + class SylvesterSolver + { + static_assert(std::is_floating_point_v, "SylvesterSolver supports floating-point types only"); + static_assert(N > 0 && M > 0, "SylvesterSolver requires positive dimensions"); + + static constexpr std::size_t NM = N * M; + + public: + SylvesterSolver() = default; + + OPTIMIZE_FOR_SPEED bool SolveSylvester( + const math::SquareMatrix& a, + const math::SquareMatrix& b, + const math::Matrix& c); + + OPTIMIZE_FOR_SPEED bool SolveContinuousLyapunov( + const math::SquareMatrix& a, + const math::SquareMatrix& q); + + OPTIMIZE_FOR_SPEED bool SolveDiscreteLyapunov( + const math::SquareMatrix& a, + const math::SquareMatrix& q); + + [[nodiscard]] const math::Matrix& Solution() const + { + return solution; + } + + private: + static math::Matrix Vec(const math::Matrix& x); + static math::Matrix UnVec(const math::Matrix& v); + + math::Matrix solution{}; + }; + + template + OPTIMIZE_FOR_SPEED bool SylvesterSolver::SolveSylvester( + const math::SquareMatrix& a, + const math::SquareMatrix& b, + const math::Matrix& c) + { + math::Matrix K{}; + for (std::size_t row = 0; row < NM; ++row) + for (std::size_t col = 0; col < NM; ++col) + { + const std::size_t i = row / N; + const std::size_t k = row % N; + const std::size_t j = col / N; + const std::size_t l = col % N; + T val{}; + if (i == j) + val += a.at(k, l); + if (k == l) + val += b.at(j, i); + K.at(row, col) = val; + } + + math::Matrix vecC = Vec(c); + + LuDecomposition lu{}; + if (!lu.Decompose(K)) + return false; + + solution = UnVec(lu.Solve(vecC)); + return true; + } + + template + OPTIMIZE_FOR_SPEED bool SylvesterSolver::SolveContinuousLyapunov( + const math::SquareMatrix& a, + const math::SquareMatrix& q) + { + static_assert(N == M, "SolveContinuousLyapunov requires square solver (N == M)"); + + return SolveSylvester(a, a.Transpose(), q * T(-1)); + } + + template + OPTIMIZE_FOR_SPEED bool SylvesterSolver::SolveDiscreteLyapunov( + const math::SquareMatrix& a, + const math::SquareMatrix& q) + { + static_assert(N == M, "SolveDiscreteLyapunov requires square solver (N == M)"); + + math::Matrix K{}; + for (std::size_t row = 0; row < NM; ++row) + for (std::size_t col = 0; col < NM; ++col) + { + const std::size_t i = row / N; + const std::size_t k = row % N; + const std::size_t j = col / N; + const std::size_t l = col % N; + T val = a.at(i, j) * a.at(k, l); + if (i == j && k == l) + val -= T(1); + K.at(row, col) = val; + } + + math::Matrix rhs = Vec(q * T(-1)); + + LuDecomposition lu{}; + if (!lu.Decompose(K)) + return false; + + solution = UnVec(lu.Solve(rhs)); + return true; + } + + template + math::Matrix SylvesterSolver::Vec(const math::Matrix& x) + { + math::Matrix v{}; + for (std::size_t i = 0; i < M; ++i) + for (std::size_t k = 0; k < N; ++k) + v.at(i * N + k, 0) = x.at(k, i); + return v; + } + + template + math::Matrix SylvesterSolver::UnVec(const math::Matrix& v) + { + math::Matrix out{}; + for (std::size_t i = 0; i < M; ++i) + for (std::size_t k = 0; k < N; ++k) + out.at(k, i) = v.at(i * N + k, 0); + return out; + } + +#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD + extern template class SylvesterSolver; + extern template class SylvesterSolver; +#endif +} diff --git a/numerical/solvers/test/CMakeLists.txt b/numerical/solvers/test/CMakeLists.txt index 7d8d2754..2a2a029b 100644 --- a/numerical/solvers/test/CMakeLists.txt +++ b/numerical/solvers/test/CMakeLists.txt @@ -14,6 +14,7 @@ target_sources(numerical.solvers_test PRIVATE TestGaussianElimination.cpp TestLevinsonDurbin.cpp TestLuDecomposition.cpp + TestLyapunovSylvester.cpp TestQrDecomposition.cpp TestRungeKuttaIntegrators.cpp TestSpectralRadius.cpp diff --git a/numerical/solvers/test/TestLyapunovSylvester.cpp b/numerical/solvers/test/TestLyapunovSylvester.cpp new file mode 100644 index 00000000..53de5adf --- /dev/null +++ b/numerical/solvers/test/TestLyapunovSylvester.cpp @@ -0,0 +1,138 @@ +#include "numerical/math/Tolerance.hpp" +#include "numerical/solvers/LyapunovSylvester.hpp" +#include + +namespace +{ + class TestSylvesterSolver : public ::testing::Test + { + protected: + solvers::SylvesterSolver sylv{}; + solvers::SylvesterSolver lyap3{}; + }; +} + +TEST_F(TestSylvesterSolver, sylvester_reconstructs_C) +{ + math::SquareMatrix A{ { -3.0f, 1.0f }, { 0.0f, -2.0f } }; + math::SquareMatrix B{ { 1.0f, 2.0f }, { 0.0f, 1.0f } }; + math::SquareMatrix C{ { 1.0f, 0.0f }, { 0.0f, 1.0f } }; + + bool ok = sylv.SolveSylvester(A, B, C); + + EXPECT_TRUE(ok); + auto X = sylv.Solution(); + auto residual = A * X + X * B; + for (std::size_t i = 0; i < 2; ++i) + for (std::size_t j = 0; j < 2; ++j) + EXPECT_NEAR(residual.at(i, j), C.at(i, j), math::Tolerance()); +} + +TEST_F(TestSylvesterSolver, continuous_lyapunov_is_spd_for_stable_A) +{ + math::SquareMatrix A{ { -2.0f, 1.0f }, { 0.0f, -3.0f } }; + math::SquareMatrix Q{ { 1.0f, 0.0f }, { 0.0f, 1.0f } }; + + bool ok = sylv.SolveContinuousLyapunov(A, Q); + + EXPECT_TRUE(ok); + auto X = sylv.Solution(); + EXPECT_GT(X.at(0, 0), 0.0f); + EXPECT_GT(X.at(1, 1), 0.0f); + EXPECT_GT(X.at(0, 0) * X.at(1, 1) - X.at(0, 1) * X.at(1, 0), 0.0f); + EXPECT_NEAR(X.at(0, 1), X.at(1, 0), math::Tolerance()); + auto residual = A * X + X * A.Transpose(); + for (std::size_t i = 0; i < 2; ++i) + for (std::size_t j = 0; j < 2; ++j) + EXPECT_NEAR(residual.at(i, j), -Q.at(i, j), math::Tolerance()); +} + +TEST_F(TestSylvesterSolver, discrete_lyapunov_matches_known) +{ + math::SquareMatrix A{ { 0.5f, 0.0f }, { 0.0f, 0.3f } }; + math::SquareMatrix Q{ { 1.0f, 0.0f }, { 0.0f, 1.0f } }; + + bool ok = sylv.SolveDiscreteLyapunov(A, Q); + + EXPECT_TRUE(ok); + auto X = sylv.Solution(); + auto residual = A * X * A.Transpose() - X; + for (std::size_t i = 0; i < 2; ++i) + for (std::size_t j = 0; j < 2; ++j) + EXPECT_NEAR(residual.at(i, j), -Q.at(i, j), math::Tolerance()); +} + +TEST_F(TestSylvesterSolver, scalar_case_closed_form) +{ + solvers::SylvesterSolver sylv1{}; + math::SquareMatrix a{ { 3.0f } }; + math::SquareMatrix c{ { 6.0f } }; + + bool ok = sylv1.SolveSylvester(a, a, c); + + EXPECT_TRUE(ok); + EXPECT_NEAR(sylv1.Solution().at(0, 0), c.at(0, 0) / (2.0f * a.at(0, 0)), math::Tolerance()); +} + +TEST_F(TestSylvesterSolver, symmetric_solution_for_symmetric_Q) +{ + math::SquareMatrix A{ { -1.0f, 0.0f, 0.0f }, { 0.0f, -2.0f, 0.0f }, { 0.0f, 0.0f, -3.0f } }; + math::SquareMatrix Q{ { 1.0f, 0.0f, 0.0f }, { 0.0f, 1.0f, 0.0f }, { 0.0f, 0.0f, 1.0f } }; + + bool ok = lyap3.SolveContinuousLyapunov(A, Q); + + EXPECT_TRUE(ok); + auto X = lyap3.Solution(); + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + EXPECT_NEAR(X.at(i, j), X.at(j, i), math::Tolerance()); +} + +TEST_F(TestSylvesterSolver, gramian_matches_truncated_sum) +{ + math::SquareMatrix A{ { 0.5f, 0.0f }, { 0.0f, 0.3f } }; + math::SquareMatrix Q{ { 1.0f, 1.0f }, { 1.0f, 1.0f } }; + + bool ok = sylv.SolveDiscreteLyapunov(A, Q); + + EXPECT_TRUE(ok); + auto X = sylv.Solution(); + + math::SquareMatrix truncSum{}; + math::SquareMatrix Ak = math::SquareMatrix::Identity(); + for (std::size_t k = 0; k <= 40; ++k) + { + truncSum += Ak * Q * Ak.Transpose(); + Ak = A * Ak; + } + + for (std::size_t i = 0; i < 2; ++i) + for (std::size_t j = 0; j < 2; ++j) + EXPECT_NEAR(X.at(i, j), truncSum.at(i, j), 1e-3f); +} + +TEST_F(TestSylvesterSolver, matches_hand_computed_2x2) +{ + math::SquareMatrix A{ { -2.0f, 0.0f }, { 0.0f, -3.0f } }; + math::SquareMatrix Q{ { 1.0f, 0.0f }, { 0.0f, 1.0f } }; + + bool ok = sylv.SolveContinuousLyapunov(A, Q); + + EXPECT_TRUE(ok); + auto X = sylv.Solution(); + EXPECT_NEAR(X.at(0, 0), 0.25f, math::Tolerance()); + EXPECT_NEAR(X.at(0, 1), 0.0f, math::Tolerance()); + EXPECT_NEAR(X.at(1, 0), 0.0f, math::Tolerance()); + EXPECT_NEAR(X.at(1, 1), 1.0f / 6.0f, math::Tolerance()); +} + +TEST_F(TestSylvesterSolver, unsolvable_pair_flagged) +{ + math::SquareMatrix A{ { 1.0f, 0.0f }, { 0.0f, 2.0f } }; + math::SquareMatrix B{ { -1.0f, 0.0f }, { 0.0f, -2.0f } }; + math::SquareMatrix C{ { 1.0f, 0.0f }, { 0.0f, 1.0f } }; + + bool ok = sylv.SolveSylvester(A, B, C); + + EXPECT_FALSE(ok); +} diff --git a/roadmap/solvers/LyapunovSylvester/explanation.md b/roadmap/solvers/LyapunovSylvester/explanation.md deleted file mode 100644 index 473fb37e..00000000 --- a/roadmap/solvers/LyapunovSylvester/explanation.md +++ /dev/null @@ -1,31 +0,0 @@ -# Lyapunov / Sylvester Equation Solvers — Overview - -## What it is -Solvers for the linear matrix equations `A X + X B = C` (Sylvester) and its symmetric relatives, the -continuous `A X + X Aᵀ = −Q` and discrete `A X Aᵀ − X = −Q` **Lyapunov** equations. The unknown is a -whole matrix `X`, not a vector. - -## Why it matters (embedded) -These equations are the algebra of stability and energy. A positive-definite Lyapunov solution is a -*certificate* that a system is stable; the same equations produce the **controllability / -observability Gramians** that quantify how well a plant can be steered and seen — prerequisites for -observers, LQR, and model reduction. - -## How it works (intuition) -Written out, `A X + X B = C` is a huge linear system in the entries of `X`. **Bartels-Stewart** avoids -that by first rotating `A` and `B` to (quasi-)triangular form via a Schur/eigen decomposition, which -*decouples* the unknown columns so they can be found one at a time by back-substitution — then the -result is rotated back. For symmetric `A` the eigenbasis diagonalizes everything and each entry is a -simple division by a sum of eigenvalues. - -## Key parameters -- **Equation type** — Sylvester vs continuous vs discrete Lyapunov. -- **Spectrum of `A`, `B`** — determines solvability (no eigenvalues that cancel). - -## Reference -R. H. Bartels, G. W. Stewart, "Solution of the Matrix Equation AX + XB = C," -*Comm. ACM*, 15(9), 1972. - -## See also -`JacobiEigenSolver`, `QrDecomposition`, `ControllabilityObservability` (Gramians), -`DiscreteAlgebraicRiccatiEquation`. diff --git a/roadmap/solvers/LyapunovSylvester/implementation.md b/roadmap/solvers/LyapunovSylvester/implementation.md deleted file mode 100644 index 0e5e8abe..00000000 --- a/roadmap/solvers/LyapunovSylvester/implementation.md +++ /dev/null @@ -1,76 +0,0 @@ -# Lyapunov / Sylvester Equation Solvers — Implementation Pseudocode - -> Roadmap ref: #31 (Tier 4) · 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 SylvesterSolver: # solves A X + X B = C - # continuous Lyapunov: B = Aᵀ, C = −Q, M = N - # discrete Lyapunov: A X Aᵀ − X = −Q - Matrix X # solution - bool solvable -``` - -## Interface - -``` -SylvesterSolver() -bool SolveSylvester(A[N×N], B[M×M], C[N×M]) # hot path -bool SolveContinuousLyapunov(A[N×N], Q[N×N]) # A X + X Aᵀ = −Q -bool SolveDiscreteLyapunov(A[N×N], Q[N×N]) # A X Aᵀ − X = −Q -const Matrix& Solution() const -``` - -## Algorithm (pseudocode) - -``` -function SolveSylvester(A, B, C): # OPTIMIZE_FOR_SPEED (Bartels-Stewart) - # 1. reduce to (quasi-)triangular form via real Schur / eigen-decomposition - (Ua, Ra) = Schur(A) # reuse QrDecomposition sweeps - (Ub, Rb) = Schur(B) - C̃ = Uaᵀ · C · Ub # transform the right-hand side - # 2. back-substitute column by column (columns decouple once triangular) - for j in 0 .. M-1: - rhs = C̃[:,j] − Σ_{k)`; the generic `T` signature keeps a - `Q15`/`Q31` specialisation cheap to add later. - -## Deployment - -- Header: `numerical/solvers/LyapunovSylvester.hpp` — `#pragma once` → - `#pragma GCC optimize("O3","fast-math")`, `OPTIMIZE_FOR_SPEED` on the solvers, and - `extern template class SylvesterSolver;` under `#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD`. -- Coverage: `numerical/solvers/LyapunovSylvester.cpp` → - `template class SylvesterSolver;` and `` (the tested sizes). -- Test: `numerical/solvers/test/TestLyapunovSylvester.cpp` -- Doc: `doc/solvers/LyapunovSylvester.md` -- CMake: `.hpp` → `target_sources`; `.cpp` → `numerical_add_coverage_sources`; - `TestLyapunovSylvester.cpp` → the `_test` target. -- Generic pattern: see `roadmap/DEPLOYMENT.md`. diff --git a/roadmap/solvers/LyapunovSylvester/tests.md b/roadmap/solvers/LyapunovSylvester/tests.md deleted file mode 100644 index 018d4981..00000000 --- a/roadmap/solvers/LyapunovSylvester/tests.md +++ /dev/null @@ -1,58 +0,0 @@ -# Lyapunov / Sylvester Solvers — Unit Test Plan (Pseudocode) - -> GoogleTest · `TEST_F` (`float`) · `StrictMock` only · no heap. - -## Fixture - -``` -class TestSylvesterSolver : public ::testing::Test: - SylvesterSolver sylv - SylvesterSolver lyap3 -# each case below is a TEST_F(TestSylvesterSolver, ) -``` - -## Test cases (Arrange / Act / Assert) - -``` -sylvester_reconstructs_C: - Arrange: fixed A, B, C - Act: SolveSylvester(A, B, C) - Assert: A·X + X·B ≈ C - -continuous_lyapunov_is_spd_for_stable_A: - Arrange: Hurwitz A, Q = I - Assert: X symmetric, positive-definite, and A·X + X·Aᵀ ≈ −Q - -discrete_lyapunov_matches_known: - Arrange: Schur-stable A, Q = I - Assert: A·X·Aᵀ − X ≈ −Q - -scalar_case_closed_form: - Arrange: 1×1 a·x + x·a = c - Assert: x ≈ c / (2a) - -symmetric_solution_for_symmetric_Q: - Assert: X ≈ Xᵀ - -gramian_matches_truncated_sum: - Arrange: controllability Gramian via discrete Lyapunov - Assert: X ≈ Σ_{k=0..K} Aᵏ B Bᵀ (Aᵀ)ᵏ (K large) - -matches_hand_computed_2x2: - Assert: X equals the documented Bartels-Stewart reference - -unsolvable_pair_flagged: - Arrange: A, B with λ_i + λ_j = 0 - Assert: Solve returns false -``` - -## Reference vectors - -- Scalar `a x + x a = c ⇒ x = c/(2a)`. -- 2×2 stable `A`, `Q = I` ⇒ hand-solved symmetric positive-definite `X`. - -## Edge cases - -- Eigenvalue pair summing to zero (continuous) or product = 1 (discrete) ⇒ no unique solution. -- Near-degenerate (tiny denominator) ⇒ ill-conditioned; flag a large solution norm. -- `N == 1` scalar; diagonal `A` closed-form cross-check.