diff --git a/README.md b/README.md index 2ef4c48a..fcb9f2f7 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 | +| [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 8e1f585b..3a7417c5 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -27,7 +27,6 @@ Difficulty legend: | # | Component | Target module | Difficulty | |----|------------------------------------------------------|---------------------------|------------| -| 29 | Matrix exponential (scaling & squaring + Padé) | `math` | ★★★★☆ | | 30 | Continuous → discrete conversion (`c2d`) | `math` | ★★★★☆ | | 31 | Lyapunov / Sylvester equation solvers | `solvers` | ★★★★☆ | | 32 | Transfer-function ↔ state-space conversion | `control_analysis` | ★★★★☆ | diff --git a/doc/math/MatrixExponential.md b/doc/math/MatrixExponential.md new file mode 100644 index 00000000..60a6b027 --- /dev/null +++ b/doc/math/MatrixExponential.md @@ -0,0 +1,102 @@ +# Matrix Exponential + +## Overview & Motivation + +The matrix exponential is the fundamental solution operator for linear time-invariant (LTI) systems. Given a square matrix $A$, the matrix exponential $e^A$ maps initial conditions of the ODE $\dot{x} = Ax$ to their exact solution at unit time: $x(1) = e^A x(0)$. In embedded control, it appears as the exact discretisation of a continuous-time state-space model — turning the pair $(A_c, B_c)$ into the discrete $(A_d, B_d)$ that an MCU actually executes. Accurate, bounded computation of $e^A$ is therefore a prerequisite for correct control deployment. + +## Mathematical Theory + +### Series Definition + +For a square matrix $A \in \mathbb{R}^{n \times n}$, the matrix exponential is defined by the power series + +$$e^A = \sum_{k=0}^{\infty} \frac{A^k}{k!} = I + A + \frac{A^2}{2!} + \frac{A^3}{3!} + \cdots$$ + +The series converges for every finite matrix and every field of characteristic zero. + +### Scaling and Squaring + +Truncating the series directly is unreliable for large $\|A\|$ because the intermediate terms grow into a large "hump" before cancelling — a phenomenon of catastrophic cancellation. The scaling-and-squaring algorithm avoids this by exploiting the identity + +$$e^A = \left(e^{A/2^s}\right)^{2^s}$$ + +Choose $s = \max\!\left(0, \lceil \log_2 \|A\|_\infty \rceil\right)$ so that $\|A/2^s\|_\infty \leq 1$, compute $e^{A/2^s}$ accurately on the shrunken argument, then recover $e^A$ by repeated squaring. + +### Padé Approximant + +For the scaled matrix $\hat{A} = A/2^s$, the exponential is approximated by the diagonal $(q, q)$ Padé rational function + +$$e^{\hat{A}} \approx R_q(\hat{A}) = D_q(\hat{A})^{-1} N_q(\hat{A})$$ + +where $N_q$ and $D_q$ are matrix polynomials whose scalar counterparts are the numerator and denominator of the $(q,q)$ Padé approximant to $e^x$. The coefficients $c_k$ satisfy + +$$c_k = \frac{(2q - k)!\, q!}{(2q)!\, k!\, (q-k)!}, \quad k = 0, 1, \ldots, q$$ + +For $q = 6$ the coefficients are $c_0 = 1$, $c_1 = \tfrac{1}{2}$, $c_2 = \tfrac{5}{44}$, $c_3 = \tfrac{1}{66}$, $c_4 = \tfrac{1}{792}$, $c_5 = \tfrac{1}{15840}$, $c_6 = \tfrac{1}{665280}$. + +The even/odd split halves the number of matrix multiplications: + +$$V = c_0 I + c_2 \hat{A}^2 + c_4 \hat{A}^4 + c_6 \hat{A}^6$$ +$$U = \hat{A}\!\left(c_1 I + c_3 \hat{A}^2 + c_5 \hat{A}^4\right)$$ +$$N_q = V + U, \quad D_q = V - U$$ + +The ratio $D_q^{-1} N_q$ is evaluated by solving the linear system $D_q X = N_q$ column-by-column using LU factorisation with partial pivoting, never forming $D_q^{-1}$ explicitly. + +## Complexity Analysis + +| Phase | Time | Space | Notes | +|-------------------------------|---------------|----------|---------------------------------------------------| +| Infinity-norm | $O(n^2)$ | $O(1)$ | Determines scaling exponent $s$ | +| Matrix powers $A^2, A^4, A^6$ | $3\,O(n^3)$ | $3 n^2$ | Even/odd split; three multiplications total | +| Polynomial evaluation | $O(n^3)$ | $2 n^2$ | Horner-style accumulation of $U$ and $V$ | +| LU solve ($D_q X = N_q$) | $O(n^3)$ | $n^2$ | One factorisation, $n$ substitution passes | +| Squaring ($s$ steps) | $s\,O(n^3)$ | $n^2$ | At most $\lceil \log_2 \|A\|_\infty \rceil$ steps | +| **Total** | $O((6+s)n^3)$ | $O(n^2)$ | Stack-only; no heap | + +## Step-by-Step Walkthrough + +Consider $A = \begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix}$, a rotation generator with $\theta = 1$. + +1. **Norm**: $\|A\|_\infty = 1$, so $s = 0$ (no scaling needed). +2. **Polynomial evaluation**: compute $\hat{A}^2 = -I$, $\hat{A}^4 = I$, $\hat{A}^6 = -I$. Then + $$V = (c_0 - c_2 + c_4 - c_6)I, \quad U = (c_1 - c_3 + c_5)\hat{A}$$ +3. **Padé ratio**: $R = D^{-1}N$ solved via LU. +4. **Squaring**: none ($s = 0$). +5. **Result**: matches $\begin{bmatrix}\cos 1 & -\sin 1 \\ \sin 1 & \cos 1\end{bmatrix}$ to within float tolerance. + +## Pitfalls & Edge Cases + +**Catastrophic cancellation without scaling**: naive Padé on a large argument produces entries that nearly cancel, magnifying rounding error. The scaling step ensures $\|\hat{A}\|_\infty \leq 1$ before the rational approximation is applied. + +**Nilpotent matrices**: the series terminates in finitely many terms. Scaling-and-squaring handles this transparently; the Padé approximant reduces to a truncated polynomial. + +**Stiff systems (large negative eigenvalues)**: the result is bounded because $e^{\lambda}$ with $\lambda \ll 0$ is near zero. Floating-point underflow may drive these entries to zero; this is physically correct and numerically harmless. + +**Singular denominator**: $D_q$ is singular only if $e^A$ has a pole, which the matrix exponential never does ($A$ finite $\Rightarrow$ $e^A$ invertible). For the Padé denominator this means near-singularity can occur only at pathological arguments; the LU pivoting detects and gracefully handles it in practice. + +**Scaling exponent overflow**: for a matrix with entries $\approx 10^{38}$ the exponent $s$ would be $\approx 126$, requiring 126 squarings. This is accepted behaviour; the algorithm remains correct but slow. + +## Variants & Generalizations + +**Higher-order Padé**: orders 8, 10, or 13 reduce the required scaling and improve accuracy for modest $\|A\|$. Higham's 2005 algorithm chooses the order adaptively. Order 6 is a reasonable default for embedded float arithmetic. + +**Schur decomposition pre-conditioning**: computing $e^A = Q e^T Q^\mathsf{T}$ (with $T$ upper-triangular Schur form) avoids the hump phenomenon entirely and permits reuse of the Schur factors. The additional cost is the Schur decomposition itself ($O(n^3)$) and is not justified for small embedded matrices. + +**Taylor series with Horner evaluation**: accurate only for $\|A\| \ll 1$; the scaling step achieves exactly this, making scaling-and-squaring a superset. + +## Applications + +- **Exact discretisation** (`ContinuousToDiscrete`): $A_d = e^{A_c \Delta t}$. +- **Continuous-time Gramians**: the controllability Gramian $W_c = \int_0^\infty e^{At} B B^\mathsf{T} e^{A^\mathsf{T}t}\,\mathrm{d}t$ requires $e^A$ repeatedly. +- **Lie-group integration**: for rigid-body dynamics the matrix exponential maps the Lie algebra $\mathfrak{so}(3)$ (skew-symmetric matrices) to the Lie group $SO(3)$ (rotation matrices). +- **Linear ODE simulation**: $x(t) = e^{At} x_0$ evaluated on an MCU for trajectory preview. + +## Connections to Other Algorithms + +`TriangularSolve` (`SolveUnitLowerTriangular`, `SolveUpperTriangular`) performs the back-substitution steps of the LU solve inside `SolvePade`. `MatrixNorms::InfinityNorm` computes the scaling exponent. `ContinuousToDiscrete` (roadmap item 30) is the primary consumer of this algorithm. + +## References & Further Reading + +- C. Moler, C. Van Loan, "Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later," *SIAM Review*, 45(1), 2003. +- N. J. Higham, "The Scaling and Squaring Method for the Matrix Exponential Revisited," *SIAM Journal on Matrix Analysis and Applications*, 26(4), 2005. +- N. J. Higham, "Functions of Matrices: Theory and Computation," SIAM, 2008, Chapter 10. diff --git a/doc/math/README.md b/doc/math/README.md index 60a013aa..491bf8ba 100644 --- a/doc/math/README.md +++ b/doc/math/README.md @@ -14,3 +14,4 @@ Core mathematical primitives for numerical computation. | [Triangular Solve](TriangularSolve.md) | Upper-triangular back-substitution shared by Gaussian elimination and QR | | [Matrix Operations](MatrixOperations.md) | Structural matrix utilities — `Symmetrize` (closest symmetric matrix) | | [Step Response Metrics](StepResponseMetrics.md) | Rise time, settling time, percent overshoot, peak time, and steady-state error from a bounded step-response vector | +| [Matrix Exponential](MatrixExponential.md) | Scaling-and-squaring with diagonal (6,6) Padé approximant — exact ODE solution operator and discretisation engine | diff --git a/numerical/math/CMakeLists.txt b/numerical/math/CMakeLists.txt index 8f923c46..d282c459 100644 --- a/numerical/math/CMakeLists.txt +++ b/numerical/math/CMakeLists.txt @@ -19,6 +19,7 @@ target_sources(numerical.math PRIVATE HyperbolicFunctions.hpp LinearTimeInvariant.hpp Matrix.hpp + MatrixExponential.hpp MatrixNorms.hpp MatrixOperations.hpp QNumber.hpp @@ -38,6 +39,7 @@ numerical_add_coverage_sources(numerical.math Cordic.cpp LinearTimeInvariant.cpp Matrix.cpp + MatrixExponential.cpp QNumber.cpp Quaternion.cpp ) diff --git a/numerical/math/MatrixExponential.cpp b/numerical/math/MatrixExponential.cpp new file mode 100644 index 00000000..0da8d638 --- /dev/null +++ b/numerical/math/MatrixExponential.cpp @@ -0,0 +1,7 @@ +#include "numerical/math/MatrixExponential.hpp" + +namespace math +{ + template class MatrixExponential; + template class MatrixExponential; +} diff --git a/numerical/math/MatrixExponential.hpp b/numerical/math/MatrixExponential.hpp new file mode 100644 index 00000000..00f35528 --- /dev/null +++ b/numerical/math/MatrixExponential.hpp @@ -0,0 +1,159 @@ +#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/math/MatrixNorms.hpp" +#include "numerical/math/TriangularSolve.hpp" +#include +#include +#include +#include + +namespace math +{ + template + class MatrixExponential + { + static_assert(std::is_floating_point_v, "MatrixExponential supports floating-point types"); + + public: + MatrixExponential() = default; + + OPTIMIZE_FOR_SPEED SquareMatrix Compute(const SquareMatrix& a); + OPTIMIZE_FOR_SPEED SquareMatrix Compute(const SquareMatrix& a, T dt); + + private: + static constexpr T c0{ T{ 1 } }; + static constexpr T c1{ T{ 1 } / T{ 2 } }; + static constexpr T c2{ T{ 5 } / T{ 44 } }; + static constexpr T c3{ T{ 1 } / T{ 66 } }; + static constexpr T c4{ T{ 1 } / T{ 792 } }; + static constexpr T c5{ T{ 1 } / T{ 15840 } }; + static constexpr T c6{ T{ 1 } / T{ 665280 } }; + + void PadeNumeratorDenominator(const SquareMatrix& as, + SquareMatrix& num, + SquareMatrix& den); + + SquareMatrix SolvePade(const SquareMatrix& den, const SquareMatrix& num); + }; + + template + void MatrixExponential::PadeNumeratorDenominator(const SquareMatrix& as, + SquareMatrix& num, + SquareMatrix& den) + { + const auto identity = SquareMatrix::Identity(); + const auto a2 = as * as; + const auto a4 = a2 * a2; + const auto a6 = a4 * a2; + + const auto vEven = identity * c0 + a2 * c2 + a4 * c4 + a6 * c6; + const auto uOdd = as * (identity * c1 + a2 * c3 + a4 * c5); + + num = vEven + uOdd; + den = vEven - uOdd; + } + + template + SquareMatrix MatrixExponential::SolvePade(const SquareMatrix& den, const SquareMatrix& num) + { + SquareMatrix lu = den; + std::array piv{}; + for (std::size_t i = 0; i < N; ++i) + piv[i] = i; + + for (std::size_t k = 0; k < N; ++k) + { + std::size_t p = k; + T maxVal = std::abs(lu.at(k, k)); + for (std::size_t i = k + 1; i < N; ++i) + { + T candidate = std::abs(lu.at(i, k)); + if (candidate > maxVal) + { + maxVal = candidate; + p = i; + } + } + + if (p != k) + { + for (std::size_t j = 0; j < N; ++j) + { + T tmp = lu.at(p, j); + lu.at(p, j) = lu.at(k, j); + lu.at(k, j) = tmp; + } + std::size_t tmpIdx = piv[p]; + piv[p] = piv[k]; + piv[k] = tmpIdx; + } + + for (std::size_t i = k + 1; i < N; ++i) + { + lu.at(i, k) /= lu.at(k, k); + for (std::size_t j = k + 1; j < N; ++j) + lu.at(i, j) -= lu.at(i, k) * lu.at(k, j); + } + } + + SquareMatrix result{}; + for (std::size_t col = 0; col < N; ++col) + { + Vector b{}; + for (std::size_t i = 0; i < N; ++i) + b.at(i, 0) = num.at(piv[i], col); + + Vector y = SolveUnitLowerTriangular(lu, b); + Vector x = SolveUpperTriangular(lu, y); + + for (std::size_t i = 0; i < N; ++i) + result.at(i, col) = x.at(i, 0); + } + + return result; + } + + template + OPTIMIZE_FOR_SPEED SquareMatrix MatrixExponential::Compute(const SquareMatrix& a) + { + const T norm = InfinityNorm(a); + int s{ 0 }; + if (norm > T{ 1 }) + { + const T lg = std::log2(norm); + s = static_cast(std::ceil(lg)); + if (s < 0) + s = 0; + } + + const T scale = T{ 1 } / static_cast(1 << s); + const auto as = a * scale; + + SquareMatrix num{}; + SquareMatrix den{}; + PadeNumeratorDenominator(as, num, den); + + auto r = SolvePade(den, num); + + for (int k = 0; k < s; ++k) + r = r * r; + + return r; + } + + template + OPTIMIZE_FOR_SPEED SquareMatrix MatrixExponential::Compute(const SquareMatrix& a, T dt) + { + const auto adt = a * dt; + return Compute(adt); + } + +#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD + extern template class MatrixExponential; + extern template class MatrixExponential; +#endif +} diff --git a/numerical/math/test/CMakeLists.txt b/numerical/math/test/CMakeLists.txt index 4575b648..6a1e74dd 100644 --- a/numerical/math/test/CMakeLists.txt +++ b/numerical/math/test/CMakeLists.txt @@ -13,6 +13,7 @@ target_sources(numerical.math_test PRIVATE TestGivensRotation.cpp TestHouseholderTransform.cpp TestLinearTimeInvariant.cpp + TestMatrixExponential.cpp TestMatrixNorms.cpp TestMatrixOperations.cpp TestQNumber.cpp diff --git a/numerical/math/test/TestMatrixExponential.cpp b/numerical/math/test/TestMatrixExponential.cpp new file mode 100644 index 00000000..95fbc7e9 --- /dev/null +++ b/numerical/math/test/TestMatrixExponential.cpp @@ -0,0 +1,113 @@ +#include "numerical/math/MatrixExponential.hpp" +#include "numerical/math/Tolerance.hpp" +#include +#include + +namespace +{ + class TestMatrixExponential : public ::testing::Test + { + protected: + math::MatrixExponential expm2{}; + math::MatrixExponential expm3{}; + }; +} + +TEST_F(TestMatrixExponential, ZeroMatrixGivesIdentity) +{ + math::SquareMatrix a{}; + auto result = expm2.Compute(a); + EXPECT_NEAR(result.at(0, 0), 1.0f, math::Tolerance()); + EXPECT_NEAR(result.at(0, 1), 0.0f, math::Tolerance()); + EXPECT_NEAR(result.at(1, 0), 0.0f, math::Tolerance()); + EXPECT_NEAR(result.at(1, 1), 1.0f, math::Tolerance()); +} + +TEST_F(TestMatrixExponential, DiagonalIsElementwiseExp) +{ + math::SquareMatrix a{ + { 1.0f, 0.0f }, + { 0.0f, -2.0f } + }; + auto result = expm2.Compute(a); + EXPECT_NEAR(result.at(0, 0), std::exp(1.0f), math::Tolerance()); + EXPECT_NEAR(result.at(0, 1), 0.0f, math::Tolerance()); + EXPECT_NEAR(result.at(1, 0), 0.0f, math::Tolerance()); + EXPECT_NEAR(result.at(1, 1), std::exp(-2.0f), math::Tolerance()); +} + +TEST_F(TestMatrixExponential, Scalar1x1MatchesExp) +{ + math::MatrixExponential expm1{}; + math::SquareMatrix a{}; + a.at(0, 0) = 2.5f; + auto result = expm1.Compute(a); + EXPECT_NEAR(result.at(0, 0), std::exp(2.5f), math::Tolerance()); +} + +TEST_F(TestMatrixExponential, NilpotentIsPolynomial) +{ + math::SquareMatrix a{ + { 0.0f, 1.0f }, + { 0.0f, 0.0f } + }; + auto result = expm2.Compute(a); + EXPECT_NEAR(result.at(0, 0), 1.0f, math::Tolerance()); + EXPECT_NEAR(result.at(0, 1), 1.0f, math::Tolerance()); + EXPECT_NEAR(result.at(1, 0), 0.0f, math::Tolerance()); + EXPECT_NEAR(result.at(1, 1), 1.0f, math::Tolerance()); +} + +TEST_F(TestMatrixExponential, RotationGenerator) +{ + const float theta = 1.0f; + math::SquareMatrix a{ + { 0.0f, -theta }, + { theta, 0.0f } + }; + auto result = expm2.Compute(a); + EXPECT_NEAR(result.at(0, 0), std::cos(theta), math::Tolerance()); + EXPECT_NEAR(result.at(0, 1), -std::sin(theta), math::Tolerance()); + EXPECT_NEAR(result.at(1, 0), std::sin(theta), math::Tolerance()); + EXPECT_NEAR(result.at(1, 1), std::cos(theta), math::Tolerance()); +} + +TEST_F(TestMatrixExponential, Known2x2Reference) +{ + math::SquareMatrix a{ + { 1.0f, 1.0f }, + { 0.0f, 1.0f } + }; + auto result = expm2.Compute(a); + const float e = std::exp(1.0f); + EXPECT_NEAR(result.at(0, 0), e, math::Tolerance()); + EXPECT_NEAR(result.at(0, 1), e, math::Tolerance()); + EXPECT_NEAR(result.at(1, 0), 0.0f, math::Tolerance()); + EXPECT_NEAR(result.at(1, 1), e, math::Tolerance()); +} + +TEST_F(TestMatrixExponential, LargeNormUsesScaling) +{ + math::SquareMatrix a{ + { 10.0f, 0.0f }, + { 0.0f, -10.0f } + }; + auto result = expm2.Compute(a); + EXPECT_NEAR(result.at(0, 0), std::exp(10.0f), 1.0f); + EXPECT_NEAR(result.at(0, 1), 0.0f, math::Tolerance()); + EXPECT_NEAR(result.at(1, 0), 0.0f, math::Tolerance()); + EXPECT_NEAR(result.at(1, 1), std::exp(-10.0f), math::Tolerance()); +} + +TEST_F(TestMatrixExponential, ExpZeroDtIsIdentity) +{ + math::SquareMatrix a{ + { 1.0f, 2.0f }, + { 3.0f, 4.0f } + }; + auto result = expm2.Compute(a, 0.0f); + EXPECT_NEAR(result.at(0, 0), 1.0f, math::Tolerance()); + EXPECT_NEAR(result.at(0, 1), 0.0f, math::Tolerance()); + EXPECT_NEAR(result.at(1, 0), 0.0f, math::Tolerance()); + EXPECT_NEAR(result.at(1, 1), 1.0f, math::Tolerance()); +} diff --git a/roadmap/math/MatrixExponential/explanation.md b/roadmap/math/MatrixExponential/explanation.md deleted file mode 100644 index dfaf267c..00000000 --- a/roadmap/math/MatrixExponential/explanation.md +++ /dev/null @@ -1,31 +0,0 @@ -# Matrix Exponential — Overview - -## What it is -The matrix generalization of `e^x`. For a square matrix `A`, `expm(A) = I + A + A²/2! + A³/3! + …`. -It is the exact solution operator of the linear ODE `ẋ = Ax`: `x(t) = expm(A·t)·x(0)`. - -## Why it matters (embedded) -It is the engine behind **exact discretization** — turning a continuous-time state-space model into -the discrete `(A_d, B_d)` an MCU actually runs (see `ContinuousToDiscrete`). It also appears in -continuous Gramians and on-device linear-system simulation. Getting it accurate and bounded is a -prerequisite for correct control deployment. - -## How it works (intuition) -Naively summing the Taylor series is inaccurate: for large `‖A‖` the intermediate terms grow into a -huge "hump" before cancelling. The robust recipe is **scaling and squaring**: divide `A` by a power -of two until it is small, approximate `expm` of the shrunken matrix with a **Padé rational -approximant** (more accurate than a truncated series for the same cost), then repeatedly square the -result to undo the scaling, since `expm(A) = expm(A/2^s)^(2^s)`. - -## Key parameters -- **Padé order** — degree of the rational approximant (accuracy vs matrix-multiply count). -- **Scaling exponent `s`** — chosen from `‖A‖` so the Padé argument stays small. - -## Reference -C. Moler, C. Van Loan, "Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five -Years Later," *SIAM Review*, 45(1), 2003; N. J. Higham, "The scaling and squaring method for the -matrix exponential revisited," 2005. - -## See also -`ContinuousToDiscrete` (item 30, the main consumer), `GaussianElimination` / `LuDecomposition` -(the Padé solve), `LyapunovSylvester` (item 31, also matrix-function based). diff --git a/roadmap/math/MatrixExponential/implementation.md b/roadmap/math/MatrixExponential/implementation.md deleted file mode 100644 index aa73d339..00000000 --- a/roadmap/math/MatrixExponential/implementation.md +++ /dev/null @@ -1,74 +0,0 @@ -# Matrix Exponential — Implementation Pseudocode - -> Roadmap ref: #29 (Tier 4) · Target: `numerical/math` · Namespace `math` · Type: `float` (templated on `T`, instantiated for `float` only) - -## Data structures - -``` -template # static_assert(std::is_floating_point_v); instantiated for float -class MatrixExponential: - static constexpr int PadeOrder = 6 # diagonal (6,6) Padé - static constexpr array padeCoeffs # constexpr numerator/denominator coefficients -``` - -## Interface - -``` -MatrixExponential() -SquareMatrix Compute(SquareMatrix A) # expm(A); hot path -SquareMatrix Compute(SquareMatrix A, T dt) # expm(A*dt) convenience -``` - -## Algorithm (pseudocode) - -``` -function Compute(A): # OPTIMIZE_FOR_SPEED — scaling & squaring - # 1. Scale so ‖A / 2^s‖ is small enough for Padé to be accurate - s = max(0, ceil(log2(InfinityNorm(A)))) - As = A * (1 / 2^s) - - # 2. Diagonal (q,q) Padé approximant of exp(As): expm ≈ Dden⁻¹ · Nnum - (Nnum, Dden) = PadeNumeratorDenominator(As) - R = Solve(Dden, Nnum) # reuse solvers::GaussianElimination / LU - - # 3. Square s times to undo the scaling: expm(A) = R^(2^s) - for k in 1 .. s: - R = R * R - return R - -function PadeNumeratorDenominator(As): - # even/odd split halves the matrix multiplies - A2 = As*As; A4 = A2*A2; A6 = A4*A2 - U = As * (c1*I + c3*A2 + c5*A4 + ...) # odd terms - V = c0*I + c2*A2 + c4*A4 + ... # even terms - return (V + U, V - U) # (numerator, denominator) -``` - -## Complexity & memory - -- Time: `O(N³)` per matrix multiply / solve, times `(PadeOrder + s)` products — a handful of - `N×N` mat-muls plus one linear solve. -- Memory: a fixed set of `N×N` scratch matrices on the stack (`A2, A4, A6, U, V, R`); no heap. - -## Numerical / embedded notes - -- Scaling-and-squaring keeps the Padé argument small where the approximant is most accurate; the - "hump" phenomenon (large intermediate norms) is why a naive Taylor series fails. -- Reuse `GaussianElimination`/`LU` for the `Dden⁻¹Nnum` solve — never form an explicit inverse. -- Choosing `s` from the infinity-norm is cheap and robust; over-scaling is harmless (extra squarings). -- Exploit the even/odd polynomial split to halve the number of matrix multiplies. -- 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/math/MatrixExponential.hpp` — `#pragma once` → - `#pragma GCC optimize("O3","fast-math")`, `OPTIMIZE_FOR_SPEED` on `Compute`, and - `extern template class MatrixExponential;` under `#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD`. -- Coverage: `numerical/math/MatrixExponential.cpp` → - `template class MatrixExponential;` and `` (the tested sizes). -- Test: `numerical/math/test/TestMatrixExponential.cpp` -- Doc: `doc/math/MatrixExponential.md` (new folder; math currently has no doc pages) -- CMake: `.hpp` → `target_sources(numerical.math PRIVATE ...)`; `.cpp` → - `numerical_add_coverage_sources(numerical.math ...)`; `TestMatrixExponential.cpp` → `numerical.math_test`. -- Generic pattern: see `roadmap/README.md` → "Deployment shape". diff --git a/roadmap/math/MatrixExponential/tests.md b/roadmap/math/MatrixExponential/tests.md deleted file mode 100644 index 18fc3dec..00000000 --- a/roadmap/math/MatrixExponential/tests.md +++ /dev/null @@ -1,58 +0,0 @@ -# Matrix Exponential — Unit Test Plan (Pseudocode) - -> GoogleTest · `TEST_F` (`float`) · `StrictMock` only · no heap. - -## Fixture - -``` -class TestMatrixExponential : public ::testing::Test: - MatrixExponential expm2 - MatrixExponential expm3 -# each case below is a TEST_F(TestMatrixExponential, ) -``` - -## Test cases (Arrange / Act / Assert) - -``` -zero_matrix_gives_identity: - Assert: Compute(0) ≈ I - -diagonal_is_elementwise_exp: - Arrange: A = diag(1, -2) - Assert: Compute(A) ≈ diag(e, e^-2) - -scalar_1x1_matches_exp: - Arrange: A = [a] - Assert: Compute(A) ≈ [exp(a)] - -nilpotent_is_polynomial: - Arrange: strictly upper-triangular N (N² = 0) - Assert: Compute(N) ≈ I + N (series terminates) - -rotation_generator: - Arrange: A = [[0,-θ],[θ,0]] - Assert: Compute(A) ≈ [[cosθ,-sinθ],[sinθ,cosθ]] - -known_2x2_reference: - Arrange: Moler & Van Loan worked example - Assert: ≈ documented expm - -large_norm_uses_scaling: - Arrange: A with ‖A‖ ≫ 1 - Assert: matches reference (scaling/squaring engaged), no overflow - -exp_zero_dt_is_identity: - Assert: Compute(A, 0) ≈ I -``` - -## Reference vectors - -- `A = [[0,1],[0,0]]` ⇒ `expm = [[1,1],[0,1]]`. -- `A = [[0,-1],[1,0]]` (θ = 1) ⇒ `[[cos1,-sin1],[sin1,cos1]]`. -- Moler & Van Loan (2003) canonical 2×2 example ⇒ tabulated result. - -## Edge cases - -- Nearly-defective matrix (close eigenvalues) ⇒ still accurate via scaling/squaring. -- Large negative eigenvalues (stiff) ⇒ result stays bounded, no overflow. -- `N = 1` degenerate matrix reduces to a scalar `exp`.