|
| 1 | +# Matrix Operations |
| 2 | + |
| 3 | +## Overview & Motivation |
| 4 | + |
| 5 | +Small structural matrix utilities that are needed across the library but do not belong to any single algorithm: |
| 6 | + |
| 7 | +- **`Symmetrize`** enforces exact symmetry on a matrix that should be symmetric in theory but drifts under floating-point round-off. |
| 8 | +- **`CongruenceTransform`** computes the quadratic form `A·M·Aᵀ` — the single most common shape in covariance-propagating code (Kalman predict `F·P·Fᵀ`, innovation covariance `H·P·Hᵀ`, Joseph update `(I−KH)·P·(I−KH)ᵀ`, etc.). |
| 9 | + |
| 10 | +Both are recurring needs in Kalman filters, EM parameter updates, and Riccati/Lyapunov solvers. |
| 11 | + |
| 12 | +## Mathematical Theory |
| 13 | + |
| 14 | +Any square matrix decomposes into a symmetric and a skew-symmetric part: |
| 15 | + |
| 16 | +$$M = \underbrace{\tfrac{1}{2}(M + M^\top)}_{\text{symmetric}} + \underbrace{\tfrac{1}{2}(M - M^\top)}_{\text{skew}}$$ |
| 17 | + |
| 18 | +`Symmetrize` returns the symmetric part $\tfrac{1}{2}(M + M^\top)$. It is the orthogonal projection (in the Frobenius inner product) of $M$ onto the subspace of symmetric matrices, so it is the *closest* symmetric matrix to $M$. Applying it to an already-symmetric matrix is a no-op (idempotent). |
| 19 | + |
| 20 | +If $M$ is symmetric, `CongruenceTransform` returns a symmetric result exactly (in exact arithmetic): $(AMA^\top)^\top = A M^\top A^\top = A M A^\top$. This makes it the natural building block for propagating a covariance $P$ through a linear map $A$: $P \mapsto A P A^\top$. |
| 21 | + |
| 22 | +## Complexity Analysis |
| 23 | + |
| 24 | +| Operation | Time | Space | Notes | |
| 25 | +|--------------------|----------|----------|----------------------------------------| |
| 26 | +| Symmetrize | $O(n^2)$ | $O(n^2)$ | One transpose, add, scale | |
| 27 | +| CongruenceTransform| $O(n^2 m + n m^2)$ | $O(nm)$ | For $A \in \mathbb{R}^{n\times m}$, two matrix products | |
| 28 | + |
| 29 | +## Step-by-Step Walkthrough |
| 30 | + |
| 31 | +**Symmetrize** — for $M = \begin{bmatrix}1 & 3 \\ -1 & 2\end{bmatrix}$: $M^\top = \begin{bmatrix}1 & -1 \\ 3 & 2\end{bmatrix}$, so $\tfrac{1}{2}(M + M^\top) = \begin{bmatrix}1 & 1 \\ 1 & 2\end{bmatrix}$ — off-diagonals averaged, diagonal unchanged. |
| 32 | + |
| 33 | +**CongruenceTransform** — with $A \in \mathbb{R}^{n\times m}$ and symmetric $M \in \mathbb{R}^{m\times m}$, the result $A M A^\top \in \mathbb{R}^{n\times n}$ is the covariance of $A x$ when $x$ has covariance $M$. |
| 34 | + |
| 35 | +## Pitfalls & Edge Cases |
| 36 | + |
| 37 | +**Symmetrize is not a fix for indefiniteness** — it removes the skew part but does not make a matrix positive-definite; covariance code typically also adds a small diagonal jitter (`+ εI`) separately. |
| 38 | + |
| 39 | +**CongruenceTransform association** — the implementation evaluates `(A·M)·Aᵀ`, matching the left-associative `operator*`; results are bit-identical to hand-written `A * M * A.Transpose()`. Under `fast-math` the symmetry of the output can still carry tiny round-off asymmetry — follow with `Symmetrize` when exact symmetry is required. |
| 40 | + |
| 41 | +**Float-only** — both are `static_assert(std::is_floating_point_v<T>)`. |
| 42 | + |
| 43 | +## Variants & Generalizations |
| 44 | + |
| 45 | +The skew-symmetric part $\tfrac{1}{2}(M - M^\top)$ is `Symmetrize`'s companion. The transposed congruence $A^\top M A$ (used by `DiscreteAlgebraicRiccatiEquation`) is obtained by passing `A.Transpose()`. |
| 46 | + |
| 47 | +## Applications |
| 48 | + |
| 49 | +- **Kalman predict / update** — `F·P·Fᵀ`, `H·P·Hᵀ`, Joseph form `(I−KH)·P·(I−KH)ᵀ + K·R·Kᵀ` across the KF/EKF/UKF/smoother family. |
| 50 | +- **EM / covariance updates** — re-symmetrize `Q`, `R`, `P` after asymmetric matrix products (`estimators::ExpectationMaximization`). |
| 51 | +- **Riccati / Lyapunov solutions** — propagate and enforce symmetry of the solution matrix. |
| 52 | + |
| 53 | +## Connections to Other Algorithms |
| 54 | + |
| 55 | +Operate on `math::Matrix` / `math::SquareMatrix`. Consumed by the `filters::active` Kalman family and `estimators::ExpectationMaximization`; `CongruenceTransform` pairs naturally with `Symmetrize` for covariance-positivity hygiene. |
| 56 | + |
| 57 | +## References & Further Reading |
| 58 | + |
| 59 | +- Golub, G. H. & Van Loan, C. F., "Matrix Computations", 4th ed., §2 (symmetric/skew decomposition) |
| 60 | +- Higham, N. J., "Accuracy and Stability of Numerical Algorithms", 2nd ed. (symmetry enforcement in covariance recursions) |
0 commit comments