feat: add qr decomposition

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, 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 |
+| [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) |
 | [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

@@ -27,7 +27,6 @@ Difficulty legend:
 
 | #  | Component                                            | Target module             | Difficulty |
 |----|------------------------------------------------------|---------------------------|------------|
-| 27 | QR decomposition (Householder / Givens)              | `solvers`                 | ★★★★☆      |
 | 28 | LU decomposition with partial pivoting               | `solvers`                 | ★★★★☆      |
 | 29 | Matrix exponential (scaling & squaring + Padé)       | `math`                    | ★★★★☆      |
 | 30 | Continuous → discrete conversion (`c2d`)             | `math`                    | ★★★★☆      |

doc/math/GivensRotation.md

@@ -0,0 +1,54 @@
+# Givens Rotation
+
+## Overview & Motivation
+
+A Givens rotation zeros a single matrix entry by rotating two rows in a plane, leaving all others untouched. Because it touches only two rows, it is the tool of choice for *incremental* linear algebra: streaming a new row into an existing QR factor, sparse triangularization, and the implicit-shift sweeps of QR/SVD eigen-iterations. `ComputeGivens` derives the rotation coefficients; `ApplyGivens` applies them to a scalar pair.
+
+## Mathematical Theory
+
+Given two scalars $a$ (the value to keep) and $b$ (the value to zero), the rotation
+
+$$G = \begin{bmatrix} c & s \\ -s & c \end{bmatrix}, \quad c = \frac{a}{r},\; s = \frac{b}{r},\; r = \sqrt{a^2 + b^2}$$
+
+is orthogonal ($c^2 + s^2 = 1$) and satisfies
+
+$$G \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix}$$
+
+Applied across a pair of rows, it zeros the target entry while preserving Euclidean length.
+
+## Complexity Analysis
+
+| Operation       | Time   | Space  | Notes                            |
+|-----------------|--------|--------|----------------------------------|
+| ComputeGivens   | $O(1)$ | $O(1)$ | One `sqrt`, two divides          |
+| ApplyGivens     | $O(1)$ | $O(1)$ | Rotates one scalar pair          |
+| Rotate two rows | $O(n)$ | $O(1)$ | `ApplyGivens` across `n` columns |
+
+## Step-by-Step Walkthrough
+
+`ComputeGivens(3, 4)`: $r = 5$, $c = 0.6$, $s = 0.8$. Applying to $(x, y) = (3, 4)$: $x' = 0.6\cdot 3 + 0.8\cdot 4 = 5$, $y' = -0.8\cdot 3 + 0.6\cdot 4 = 0$ — the second component is annihilated and the norm is preserved.
+
+## Pitfalls & Edge Cases
+
+**Degenerate pair** — when $a = b = 0$ the rotation is undefined; `ComputeGivens` returns the identity $(c, s) = (1, 0)$ so applying it is a safe no-op.
+
+**Float-only** — `static_assert(std::is_floating_point_v<T>)`.
+
+## Variants & Generalizations
+
+Householder reflectors (`math::HouseholderTransform`) zero an entire sub-column at once and are cheaper for dense factorization; Givens wins when only one entry (or one streamed row) changes. Fast/"square-root-free" Givens variants trade the `sqrt` for extra bookkeeping.
+
+## Applications
+
+- **Streaming QR update** — rotate a new row into an existing `R` (`QrDecomposition::GivensUpdateRow`).
+- **QR / SVD iterations** — implicit-shift bulge chasing.
+- **Sparse triangularization** — zero isolated entries without touching the rest.
+
+## Connections to Other Algorithms
+
+Used by `solvers::QrDecomposition`; complements `math::HouseholderTransform`.
+
+## References & Further Reading
+
+- Givens, W., "Computation of Plane Unitary Rotations Transforming a General Matrix to Triangular Form", SIAM J. Appl. Math. 6(1), 1958
+- Golub, G. H. & Van Loan, C. F., "Matrix Computations", 4th ed., §5.1 (Givens rotations)

doc/math/HouseholderTransform.md

@@ -0,0 +1,55 @@
+# Householder Transform
+
+## Overview & Motivation
+
+A Householder reflector is the workhorse of numerically stable dense linear algebra. A single reflector mirrors a vector onto a coordinate axis, zeroing every entry below a chosen pivot in one orthogonal step. Chaining reflectors triangularizes a matrix (QR), bidiagonalizes it (SVD), or tridiagonalizes a symmetric matrix (eigensolvers). `HouseholderVector` computes the reflector for one sub-column; it is the shared primitive those factorizations call.
+
+## Mathematical Theory
+
+For a vector $x$, the Householder reflector is the orthogonal matrix
+
+$$H = I - \beta\, v v^\top$$
+
+chosen so that $Hx$ is zero below the pivot. With $\sigma = \sum_{i>\text{start}} x_i^2$ and $\|x\| = \sqrt{x_{\text{start}}^2 + \sigma}$, the reflector maps $x_{\text{start}} \mapsto \mp\|x\|$. The pivot sign is chosen as $-\operatorname{sign}(x_{\text{start}})\|x\|$ to avoid cancellation:
+
+$$v_{\text{start}} = 1, \quad v_i = x_i / v_0, \quad \beta = \frac{2 v_0^2}{\sigma + v_0^2}$$
+
+$H$ is symmetric and orthogonal ($H = H^\top = H^{-1}$), so applying it is backward stable.
+
+## Complexity Analysis
+
+| Operation             | Time   | Space  | Notes                                   |
+|-----------------------|--------|--------|-----------------------------------------|
+| HouseholderVector     | $O(n)$ | $O(1)$ | Builds `v` (stored implicit unit pivot) |
+| Apply $H$ to a vector | $O(n)$ | $O(1)$ | `x - β·v·(vᵀx)` — never form `H`        |
+
+## Step-by-Step Walkthrough
+
+For $x = [4, 3, 0, 0]^\top$, pivot `start = 0`: $\sigma = 9$, $\|x\| = 5$. Since $x_0 > 0$, $v_0 = -\sigma/(x_0 + \|x\|) = -1$, giving $\beta = 1$ and $v = [1, -3, 0, 0]^\top$. Reflecting yields $Hx = [-5, 0, 0, 0]^\top$ — the sub-column collapsed onto the axis.
+
+## Pitfalls & Edge Cases
+
+**Already-zero sub-column** — when $\sigma \approx 0$ there is nothing to zero; the routine returns $\beta = 0$ (identity reflector) and callers skip the update.
+
+**Never form `H` explicitly** — apply it as `x − β·v·(vᵀx)` to keep the cost $O(n)$ per column instead of $O(n^2)$.
+
+**Float-only** — `static_assert(std::is_floating_point_v<T>)`; the sign-fixing and normalisation assume real floating-point arithmetic.
+
+## Variants & Generalizations
+
+Givens rotations (`math::GivensRotation`) achieve the same zeroing one entry at a time, preferable for sparse or streaming updates. Complex Householder reflectors extend this to unitary triangularization.
+
+## Applications
+
+- **QR decomposition** — one reflector per column triangularizes `A`.
+- **SVD / eigensolvers** — bidiagonalization and tridiagonalization.
+- **Square-root Kalman filtering** — covariance factor updates.
+
+## Connections to Other Algorithms
+
+Used by `solvers::QrDecomposition`; complements `math::GivensRotation` and `math::SolveUpperTriangular`. Operates on `math::Vector`.
+
+## References & Further Reading
+
+- Householder, A. S., "Unitary Triangularization of a Nonsymmetric Matrix", JACM 5(4), 1958
+- Golub, G. H. & Van Loan, C. F., "Matrix Computations", 4th ed., §5.1 (Householder reflections)

doc/math/MatrixNorms.md

@@ -8,13 +8,13 @@ Vector normalisation, closely related, produces the unit-length direction of a v
 
 ## Mathematical Theory
 
-### Vector Norm
+### Dot Product & Vector Norm
 
-For a vector $\mathbf{v} \in \mathbb{R}^n$, the Euclidean (L2) norm is
+For vectors $\mathbf{a}, \mathbf{b} \in \mathbb{R}^n$, the dot product is $\mathbf{a}\cdot\mathbf{b} = \sum_{i=1}^{n} a_i b_i$. The Euclidean (L2) norm is the square root of the self dot product:
 
-$$\|\mathbf{v}\|_2 = \sqrt{\sum_{i=1}^{n} v_i^2}$$
+$$\|\mathbf{v}\|_2 = \sqrt{\mathbf{v}\cdot\mathbf{v}} = \sqrt{\sum_{i=1}^{n} v_i^2}$$
 
-The unit vector $\hat{\mathbf{v}} = \mathbf{v} / \|\mathbf{v}\|_2$ satisfies $\|\hat{\mathbf{v}}\|_2 = 1$. Normalisation is undefined when $\|\mathbf{v}\|_2 = 0$ and must be guarded.
+`VectorNorm` is implemented in terms of `DotProduct`. The unit vector $\hat{\mathbf{v}} = \mathbf{v} / \|\mathbf{v}\|_2$ satisfies $\|\hat{\mathbf{v}}\|_2 = 1$. Normalisation is undefined when $\|\mathbf{v}\|_2 = 0$ and must be guarded.
 
 ### Matrix Norms
 
@@ -36,8 +36,9 @@ The 1-norm and infinity-norm are dual: $\|A\|_\infty = \|A^\top\|_1$.
 
 ## Complexity Analysis
 
-| Operation     | Time    | Space  | Notes                        |
-|---------------|---------|--------|------------------------------|
+| Operation     | Time    | Space  | Notes                                  |
+|---------------|---------|--------|----------------------------------------|
+| DotProduct    | $O(n)$  | $O(1)$ | Single pass; `VectorNorm` builds on it |
 | FrobeniusNorm | $O(mn)$ | $O(1)$ | Single pass, no allocation   |
 | OneNorm       | $O(mn)$ | $O(1)$ | Column-wise sum, running max |
 | InfinityNorm  | $O(mn)$ | $O(1)$ | Row-wise sum, running max    |

doc/math/MatrixOperations.md

@@ -0,0 +1,60 @@
+# Matrix Operations
+
+## Overview & Motivation
+
+Small structural matrix utilities that are needed across the library but do not belong to any single algorithm:
+
+- **`Symmetrize`** enforces exact symmetry on a matrix that should be symmetric in theory but drifts under floating-point round-off.
+- **`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.).
+
+Both are recurring needs in Kalman filters, EM parameter updates, and Riccati/Lyapunov solvers.
+
+## Mathematical Theory
+
+Any square matrix decomposes into a symmetric and a skew-symmetric part:
+
+$$M = \underbrace{\tfrac{1}{2}(M + M^\top)}_{\text{symmetric}} + \underbrace{\tfrac{1}{2}(M - M^\top)}_{\text{skew}}$$
+
+`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).
+
+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$.
+
+## Complexity Analysis
+
+| Operation           | Time               | Space    | Notes                                                   |
+|---------------------|--------------------|----------|---------------------------------------------------------|
+| Symmetrize          | $O(n^2)$           | $O(n^2)$ | One transpose, add, scale                               |
+| CongruenceTransform | $O(n^2 m + n m^2)$ | $O(nm)$  | For $A \in \mathbb{R}^{n\times m}$, two matrix products |
+
+## Step-by-Step Walkthrough
+
+**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.
+
+**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$.
+
+## Pitfalls & Edge Cases
+
+**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.
+
+**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.
+
+**Float-only** — both are `static_assert(std::is_floating_point_v<T>)`.
+
+## Variants & Generalizations
+
+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()`.
+
+## Applications
+
+- **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.
+- **EM / covariance updates** — re-symmetrize `Q`, `R`, `P` after asymmetric matrix products (`estimators::ExpectationMaximization`).
+- **Riccati / Lyapunov solutions** — propagate and enforce symmetry of the solution matrix.
+
+## Connections to Other Algorithms
+
+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.
+
+## References & Further Reading
+
+- Golub, G. H. & Van Loan, C. F., "Matrix Computations", 4th ed., §2 (symmetric/skew decomposition)
+- Higham, N. J., "Accuracy and Stability of Numerical Algorithms", 2nd ed. (symmetry enforcement in covariance recursions)

doc/math/README.md

@@ -9,4 +9,8 @@ Core mathematical primitives for numerical computation.
 | [CORDIC](Cordic.md)                                 | Iterative shift-add engine for sin/cos, atan2, magnitude, and vector rotation — no multiplier |
 | [Quaternion](Quaternion.md)                         | Unit-quaternion rotation type: Hamilton product, SLERP, rotation-matrix and Euler conversions |
 | [MatrixNorms](MatrixNorms.md)                       | Frobenius, 1-norm, infinity-norm on matrices; vector L2 norm/normalize                        |
+| [Householder Transform](HouseholderTransform.md)    | Householder reflector for a sub-column — orthogonal, backward-stable factorization primitive  |
+| [Givens Rotation](GivensRotation.md)                | Plane rotation zeroing one entry — streaming/sparse factorization primitive                   |
+| [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 |

doc/math/TriangularSolve.md

@@ -0,0 +1,52 @@
+# Triangular Solve
+
+## Overview & Motivation
+
+Nearly every dense linear solver ends with a triangular system. Gaussian elimination reduces `Ax = b` to an upper-triangular `Ux = c`; QR factorization solves least squares through `Rx = Qᵀb`; Cholesky solves two triangular systems back-to-back. `SolveUpperTriangular` is the shared back-substitution kernel these routines call so the algorithm lives in exactly one place.
+
+## Mathematical Theory
+
+Given an upper-triangular matrix $R \in \mathbb{R}^{n \times n}$ (entries below the diagonal ignored) and a right-hand side $c$, back-substitution solves $Rx = c$ bottom-up:
+
+$$x_i = \frac{1}{r_{ii}}\left(c_i - \sum_{j=i+1}^{n} r_{ij}\, x_j\right), \quad i = n, n-1, \dots, 1$$
+
+Each unknown depends only on those already computed, so a single reverse sweep suffices.
+
+## Complexity Analysis
+
+| Operation            | Time     | Space  | Notes                           |
+|----------------------|----------|--------|---------------------------------|
+| SolveUpperTriangular | $O(n^2)$ | $O(n)$ | One reverse sweep, stack output |
+
+## Step-by-Step Walkthrough
+
+For $R = \begin{bmatrix}2 & -1 & 3 \\ 0 & 4 & 1 \\ 0 & 0 & 5\end{bmatrix}$, $c = [10,\, -5,\, 15]^\top$:
+
+1. $x_3 = 15 / 5 = 3$
+2. $x_2 = (-5 - 1\cdot 3)/4 = -2$
+3. $x_1 = (10 - (-1)(-2) - 3\cdot 3)/2 = 1$ ⇒ $x = [1,\, -2,\, 3]^\top$.
+
+## Pitfalls & Edge Cases
+
+**Singular / near-zero pivot** — a zero diagonal entry makes the system unsolvable. The routine asserts `|r_{ii}| > 0` via `really_assert`; callers (`GaussianElimination`, `QrDecomposition`) detect rank deficiency before reaching it.
+
+**Generic on `T`** — the kernel is templated on any supported numeric type (`float`, `Q15`, `Q31`) because Gaussian elimination is instantiated for all three; it uses `math::ToFloat` only for the pivot assertion.
+
+## Variants & Generalizations
+
+A lower-triangular forward-substitution is the mirror image (top-down sweep) and can be added when Cholesky/LU forward solves need it. Block triangular solves generalise this to matrix right-hand sides.
+
+## Applications
+
+- **Gaussian elimination** — final back-substitution after forward elimination.
+- **QR least squares** — solving `Rx = Qᵀb`.
+- Any factor-then-solve routine producing a triangular factor.
+
+## Connections to Other Algorithms
+
+Shared by `solvers::GaussianElimination` and `solvers::QrDecomposition`. Operates on `math::Matrix` / `math::Vector`.
+
+## References & Further Reading
+
+- Golub, G. H. & Van Loan, C. F., "Matrix Computations", 4th ed., §3.1 (triangular systems)
+- Trefethen, L. N. & Bau, D., "Numerical Linear Algebra", Lecture 17