diff --git a/README.md b/README.md
index 212c0dd3..44258e5e 100644
--- a/README.md
+++ b/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.
diff --git a/ROADMAP.md b/ROADMAP.md
index e6452043..c9473cb9 100644
--- a/ROADMAP.md
+++ b/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`                    | ★★★★☆      |
diff --git a/doc/math/GivensRotation.md b/doc/math/GivensRotation.md
new file mode 100644
index 00000000..19ab7ca5
--- /dev/null
+++ b/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)
diff --git a/doc/math/HouseholderTransform.md b/doc/math/HouseholderTransform.md
new file mode 100644
index 00000000..e5896262
--- /dev/null
+++ b/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)
diff --git a/doc/math/MatrixNorms.md b/doc/math/MatrixNorms.md
index 15b0b6a4..9df1653b 100644
--- a/doc/math/MatrixNorms.md
+++ b/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    |
diff --git a/doc/math/MatrixOperations.md b/doc/math/MatrixOperations.md
new file mode 100644
index 00000000..cc7888c7
--- /dev/null
+++ b/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)
diff --git a/doc/math/README.md b/doc/math/README.md
index 456ea5a8..60a013aa 100644
--- a/doc/math/README.md
+++ b/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 |
diff --git a/doc/math/TriangularSolve.md b/doc/math/TriangularSolve.md
new file mode 100644
index 00000000..0df497e9
--- /dev/null
+++ b/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
diff --git a/doc/solvers/QrDecomposition.md b/doc/solvers/QrDecomposition.md
new file mode 100644
index 00000000..c0e49523
--- /dev/null
+++ b/doc/solvers/QrDecomposition.md
@@ -0,0 +1,125 @@
+# QR Decomposition
+
+## Overview & Motivation
+
+The QR decomposition factors a matrix $A$ into an orthogonal factor $Q$ and an upper-triangular
+factor $R$. It is the numerically robust foundation for least-squares regression, eigenvalue
+computation, and square-root filtering. In embedded systems, fitting an overdetermined model
+(more sensor measurements than unknowns) arises constantly in calibration, system identification,
+and sensor fusion. Solving the least-squares problem via QR avoids forming the normal equations
+$A^\top A x = A^\top b$, whose condition number is the square of $A$'s — a critical advantage
+when working with reduced word lengths.
+
+## Mathematical Theory
+
+### Definitions
+
+For an $m \times n$ matrix $A$ with $m \geq n$, the thin QR factorization is
+
+$$A = Q R,$$
+
+where $Q \in \mathbb{R}^{m \times n}$ has orthonormal columns ($Q^\top Q = I_n$) and
+$R \in \mathbb{R}^{n \times n}$ is upper triangular.
+
+### Householder Reflectors
+
+A Householder reflector is a symmetric, orthogonal matrix of the form
+
+$$H = I - \beta v v^\top, \quad \beta = \frac{2}{\|v\|^2},$$
+
+chosen so that $H x = \sigma e_1$ for a target vector $x$, where $\sigma = -\operatorname{sign}(x_1)\|x\|$.
+The sign convention $\sigma = -\operatorname{sign}(x_1)\|x\|$ avoids catastrophic cancellation
+when computing $v_1 = x_1 - \sigma$.
+
+Applying $k$ successive reflectors $H_k, \ldots, H_1$ to $A$ triangularizes it:
+
+$$H_n \cdots H_1 A = R \implies A = H_1 \cdots H_n R = QR.$$
+
+The reflectors are stored implicitly: below the diagonal of the working matrix, alongside a
+scalar $\beta_k$ per column, consuming only $O(mn)$ words.
+
+### Givens Rotations for Streaming Updates
+
+A Givens rotation
+
+$$G(i,j,\theta) = \begin{pmatrix} c & s \\ -s & c \end{pmatrix}, \quad c = \cos\theta,\ s = \sin\theta$$
+
+zeros a single sub-diagonal entry by rotating two rows. When one new row streams in after an
+existing factorization, $n$ plane rotations update $R$ in $O(n^2)$ flops — far cheaper than
+re-factorizing the augmented matrix from scratch.
+
+### Least-Squares Solve
+
+Given the factorization $A = QR$, the least-squares minimum-norm solution is
+
+$$x = R^{-1} Q^\top b,$$
+
+computed as: (1) multiply $b \leftarrow Q^\top b$ by applying the stored reflectors in order
+without forming $Q$ explicitly, then (2) back-substitute $R x = (Q^\top b)_{1:n}$.
+
+## Complexity Analysis
+
+| Case    | Time                          | Space    | Notes                                           |
+|---------|-------------------------------|----------|-------------------------------------------------|
+| Best    | $O(mn^2)$                     | $O(mn)$  | Householder on a well-conditioned dense $A$     |
+| Average | $O(2mn^2 - \tfrac{2}{3}n^3)$  | $O(mn)$  | Standard flop count from Golub & Van Loan       |
+| Worst   | same as average (no pivoting) | $O(mn)$  | Rank-deficient case detected but not re-ordered |
+| Update  | $O(n^2)$                      | $O(n^2)$ | Single Givens row update; no re-factor          |
+
+All memory is stack-allocated via `std::array`; the factor is stored in-place over the input.
+
+## Step-by-Step Walkthrough
+
+Factor $A = \begin{pmatrix} 12 & -51 \\ 6 & 167 \\ -4 & 24 \end{pmatrix}$ (first two columns of the
+Golub & Van Loan example).
+
+**Column 0 ($k=0$):** Extract $x = (12, 6, -4)^\top$, $\|x\| = 14$. Choose $\sigma = -14$
+(sign rule), $v_0 = 12 + 14 = 26$, $v = (26, 6, -4)^\top / 26$, $\beta = 2/\|v\|^2$.
+After applying $H_0$: $R_{00} = 14$, zeros below.
+
+**Column 1 ($k=1$):** Repeat on the $2 \times 1$ trailing sub-column.
+
+The upper-triangular $R$ then has diagonal magnitudes $\{14, 175\}$ for this sub-matrix.
+
+## Pitfalls & Edge Cases
+
+- **Sign convention**: always use $\sigma = -\operatorname{sign}(x_1)\|x\|$ to prevent the
+  subtraction $x_1 - \sigma$ from cancelling when $x_1 > 0$.
+- **Rank deficiency**: a near-zero diagonal entry in $R$ signals a dependent column; the
+  factorization stops and returns false rather than dividing by zero.
+- **Accumulating $Q$**: forming $Q$ explicitly costs an extra $O(mn^2)$ pass and
+  $O(mn)$ storage; `ApplyQtranspose` avoids this for solves.
+- **Streaming updates**: Givens rotations maintain numerical orthogonality one rotation at a time;
+  round-off does not accumulate as severely as re-using a stale Householder factor.
+- **Square matrices**: thin QR with $m = n$ still works; $Q$ is $n \times n$ orthogonal.
+
+## Variants & Generalizations
+
+- **Column pivoting** (QRP): permutes columns to move the largest-norm column to the front,
+  yielding a rank-revealing factorization useful for near-rank-deficient $A$.
+- **Gram-Schmidt (MGS)**: equivalent in exact arithmetic but numerically inferior; only suitable
+  when $A$ is already well-conditioned.
+- **Block Householder**: batches several reflectors into a $WY$ representation for cache-efficient
+  BLAS-3 level operations on large matrices.
+
+## Applications
+
+- Overdetermined least-squares regression: calibration curves, polynomial fits, sensor fusion.
+- Eigenvalue computation: the QR algorithm iterates $A_{k+1} = R_k Q_k$ to convergence.
+- Square-root Kalman filtering: the covariance propagation uses a QR update of the augmented
+  square-root factor rather than an explicit $PFP^\top + Q$ product.
+- Singular Value Decomposition: bidiagonalizes $A$ via two-sided Householder as a first step.
+
+## Connections to Other Algorithms
+
+- Provides the factorization underpinning `SingularValueDecomposition` (Golub-Kahan) and the
+  symmetric eigenvalue solver (Jacobi / tridiagonal QR).
+- The back-substitution step is shared with `GaussianElimination`'s upper-triangular solve.
+- Supersedes the normal equations used by `LinearRegression` for ill-conditioned problems.
+- The square-root Kalman filter stores $P^{1/2}$ and updates it with a Givens-based QR update.
+
+## References & Further Reading
+
+- A. S. Householder, "Unitary Triangularization of a Nonsymmetric Matrix," *JACM*, 5(4), 1958.
+- G. H. Golub & C. F. Van Loan, *Matrix Computations*, 4th ed., Johns Hopkins UP, 2013, Ch. 5.
+- W. H. Press et al., *Numerical Recipes in C*, 3rd ed., Cambridge UP, 2007, §2.10.
diff --git a/doc/solvers/README.md b/doc/solvers/README.md
index 8513d223..4a3a5598 100644
--- a/doc/solvers/README.md
+++ b/doc/solvers/README.md
@@ -14,3 +14,4 @@ Numerical solvers for linear systems, polynomial roots, and matrix equations.
 | [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         |
+| [QR Decomposition](QrDecomposition.md)                                     | Householder factorization and Givens streaming row update for least-squares solves |
diff --git a/numerical/estimators/ConsistencyMetrics.hpp b/numerical/estimators/ConsistencyMetrics.hpp
index 3c56cc9f..59f953d2 100644
--- a/numerical/estimators/ConsistencyMetrics.hpp
+++ b/numerical/estimators/ConsistencyMetrics.hpp
@@ -6,6 +6,7 @@
 
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/math/Matrix.hpp"
+#include "numerical/math/MatrixNorms.hpp"
 #include "numerical/solvers/GaussianElimination.hpp"
 #include <array>
 #include <cmath>
@@ -48,7 +49,6 @@ namespace estimators
 
     private:
         [[nodiscard]] static OPTIMIZE_FOR_SPEED std::optional<StateVector> Solve(const CovarianceMatrix& matrix, const StateVector& rhs);
-        [[nodiscard]] static OPTIMIZE_FOR_SPEED T DotProduct(const StateVector& a, const StateVector& b);
     };
 
     template<typename T, std::size_t Dim>
@@ -57,7 +57,7 @@ namespace estimators
         auto z = Solve(covariance, error);
         if (!z.has_value())
             return std::nullopt;
-        return DotProduct(error, z.value());
+        return math::DotProduct(error, z.value());
     }
 
     template<typename T, std::size_t Dim>
@@ -66,7 +66,7 @@ namespace estimators
         auto z = Solve(innovationCovariance, innovation);
         if (!z.has_value())
             return std::nullopt;
-        return DotProduct(innovation, z.value());
+        return math::DotProduct(innovation, z.value());
     }
 
     template<typename T, std::size_t Dim>
@@ -106,15 +106,6 @@ namespace estimators
         return solver.Solve(matrix, rhs);
     }
 
-    template<typename T, std::size_t Dim>
-    OPTIMIZE_FOR_SPEED T ConsistencyMetrics<T, Dim>::DotProduct(const StateVector& a, const StateVector& b)
-    {
-        T result{};
-        for (std::size_t i = 0; i < Dim; ++i)
-            result += a.at(i, 0) * b.at(i, 0);
-        return result;
-    }
-
 #ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
     extern template class ConsistencyMetrics<float, 1>;
     extern template class ConsistencyMetrics<float, 2>;
diff --git a/numerical/estimators/offline/ExpectationMaximization.hpp b/numerical/estimators/offline/ExpectationMaximization.hpp
index 22bc237a..b93ed3e7 100644
--- a/numerical/estimators/offline/ExpectationMaximization.hpp
+++ b/numerical/estimators/offline/ExpectationMaximization.hpp
@@ -7,6 +7,7 @@
 #include "numerical/filters/active/KalmanSmoother.hpp"
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/math/Matrix.hpp"
+#include "numerical/math/MatrixOperations.hpp"
 #include "numerical/solvers/GaussianElimination.hpp"
 #include <array>
 #include <cmath>
@@ -78,9 +79,6 @@ namespace estimators
             const StateMatrix& fullStateCov,
             const MeasurementCovariance& obsOuterSum,
             std::size_t numSteps) const;
-
-        template<std::size_t N>
-        static math::SquareMatrix<float, N> Symmetrize(const math::SquareMatrix<float, N>& M);
     };
 
     // Implementation //
@@ -191,7 +189,7 @@ namespace estimators
             laggedStateCov, transitionCrossCov.Transpose())
                                .Transpose();
         const float invTm1 = 1.0f / static_cast<float>(numSteps - 1);
-        auto Q_new = Symmetrize<StateSize>((currentStateCov - F_new * transitionCrossCov.Transpose()) * invTm1);
+        auto Q_new = math::Symmetrize((currentStateCov - F_new * transitionCrossCov.Transpose()) * invTm1);
         Q_new += StateMatrix::Identity() * 1e-6f;
         return { F_new, Q_new };
     }
@@ -210,19 +208,11 @@ namespace estimators
             fullStateCov, obsStateCrossCov.Transpose())
                                .Transpose();
         const float invT = 1.0f / static_cast<float>(numSteps);
-        auto R_new = Symmetrize<MeasurementSize>((obsOuterSum - H_new * obsStateCrossCov.Transpose()) * invT);
+        auto R_new = math::Symmetrize((obsOuterSum - H_new * obsStateCrossCov.Transpose()) * invT);
         R_new += MeasurementCovariance::Identity() * 1e-6f;
         return { H_new, R_new };
     }
 
-    template<std::size_t StateSize, std::size_t MeasurementSize, std::size_t MaxSteps>
-    template<std::size_t N>
-    math::SquareMatrix<float, N>
-    ExpectationMaximization<StateSize, MeasurementSize, MaxSteps>::Symmetrize(const math::SquareMatrix<float, N>& M)
-    {
-        return (M + M.Transpose()) * 0.5f;
-    }
-
 #ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
     extern template class ExpectationMaximization<2, 1, 10>;
     extern template class ExpectationMaximization<4, 2, 20>;
diff --git a/numerical/estimators/offline/LinearRegression.hpp b/numerical/estimators/offline/LinearRegression.hpp
index 1d281ad8..c45a35e0 100644
--- a/numerical/estimators/offline/LinearRegression.hpp
+++ b/numerical/estimators/offline/LinearRegression.hpp
@@ -8,7 +8,7 @@
 
 #include "numerical/estimators/Estimator.hpp"
 #include "numerical/math/CompilerOptimizations.hpp"
-#include "numerical/solvers/GaussianElimination.hpp"
+#include "numerical/solvers/QrDecomposition.hpp"
 
 namespace estimators
 {
@@ -40,17 +40,15 @@ namespace estimators
 
         for (size_t i = 0; i < Samples; ++i)
         {
-            X_design.at(i, 0) = T(0.9999f);
+            X_design.at(i, 0) = T{ 1 };
 
             for (size_t j = 0; j < Features; ++j)
                 X_design.at(i, j + 1) = X.at(i, j);
         }
 
-        auto X_transpose = X_design.Transpose();
-        auto XtX = X_transpose * X_design;
-        auto Xty = X_transpose * y;
-
-        coefficients = solvers::SolveSystem<T, Features + 1, 1>(XtX, Xty);
+        solvers::QrDecomposition<T, Samples, Features + 1> qr;
+        qr.Decompose(X_design);
+        coefficients = qr.SolveLeastSquares(y);
     }
 
     template<typename T, std::size_t Samples, std::size_t Features>
diff --git a/numerical/estimators/offline/PolynomialFitting.hpp b/numerical/estimators/offline/PolynomialFitting.hpp
index 052a60e7..491cb655 100644
--- a/numerical/estimators/offline/PolynomialFitting.hpp
+++ b/numerical/estimators/offline/PolynomialFitting.hpp
@@ -6,7 +6,7 @@
 
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/math/Matrix.hpp"
-#include "numerical/solvers/GaussianElimination.hpp"
+#include "numerical/solvers/QrDecomposition.hpp"
 #include <type_traits>
 
 namespace estimators
@@ -43,11 +43,9 @@ namespace estimators
                 v.at(i, j) = v.at(i, j - 1) * x.at(i, 0);
         }
 
-        auto vt = v.Transpose();
-        auto normalMatrix = vt * v;
-        auto rhs = vt * y;
-
-        coefficients = solvers::SolveSystem<T, Degree + 1, 1>(normalMatrix, rhs);
+        solvers::QrDecomposition<T, Samples, Degree + 1> qr;
+        qr.Decompose(v);
+        coefficients = qr.SolveLeastSquares(y);
     }
 
     template<typename T, std::size_t Samples, std::size_t Degree>
diff --git a/numerical/filters/active/ExtendedKalmanFilter.hpp b/numerical/filters/active/ExtendedKalmanFilter.hpp
index ce1b6225..31d17f1c 100644
--- a/numerical/filters/active/ExtendedKalmanFilter.hpp
+++ b/numerical/filters/active/ExtendedKalmanFilter.hpp
@@ -3,6 +3,7 @@
 #include "infra/util/Function.hpp"
 #include "numerical/filters/active/KalmanFilterBase.hpp"
 #include "numerical/math/CompilerOptimizations.hpp"
+#include "numerical/math/MatrixOperations.hpp"
 
 #if defined(__GNUC__) || defined(__clang__)
 #pragma GCC optimize("O3", "fast-math")
@@ -94,7 +95,7 @@ namespace filters
     {
         auto F = stateJacobianFn(this->state());
         this->state() = stateTransitionFn(this->state());
-        this->covariance() = F * this->covariance() * F.Transpose() + this->processNoise();
+        this->covariance() = math::CongruenceTransform(F, this->covariance()) + this->processNoise();
     }
 
     template<typename QNumberType, std::size_t StateSize, std::size_t MeasurementSize, std::size_t ControlSize>
@@ -103,7 +104,7 @@ namespace filters
     {
         auto F = stateJacobianWithControlFn(this->state(), u);
         this->state() = stateTransitionWithControlFn(this->state(), u);
-        this->covariance() = F * this->covariance() * F.Transpose() + this->processNoise();
+        this->covariance() = math::CongruenceTransform(F, this->covariance()) + this->processNoise();
     }
 
     template<typename QNumberType, std::size_t StateSize, std::size_t MeasurementSize, std::size_t ControlSize>
diff --git a/numerical/filters/active/KalmanFilter.hpp b/numerical/filters/active/KalmanFilter.hpp
index 48ad312f..5c5d1d02 100644
--- a/numerical/filters/active/KalmanFilter.hpp
+++ b/numerical/filters/active/KalmanFilter.hpp
@@ -3,6 +3,7 @@
 #include "numerical/filters/active/KalmanFilterBase.hpp"
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/math/LinearTimeInvariant.hpp"
+#include "numerical/math/MatrixOperations.hpp"
 
 #if defined(__GNUC__) || defined(__clang__)
 #pragma GCC optimize("O3", "fast-math")
@@ -108,7 +109,7 @@ namespace filters
     OPTIMIZE_FOR_SPEED void KalmanFilter<QNumberType, StateSize, MeasurementSize, ControlSize>::Predict()
     {
         this->state() = stateTransition * this->state();
-        this->covariance() = stateTransition * this->covariance() * stateTransition.Transpose() + this->processNoise();
+        this->covariance() = math::CongruenceTransform(stateTransition, this->covariance()) + this->processNoise();
     }
 
     template<typename QNumberType, std::size_t StateSize, std::size_t MeasurementSize, std::size_t ControlSize>
@@ -116,7 +117,7 @@ namespace filters
     requires(ControlSize > 0)
     {
         this->state() = stateTransition * this->state() + controlInputMatrix * u;
-        this->covariance() = stateTransition * this->covariance() * stateTransition.Transpose() + this->processNoise();
+        this->covariance() = math::CongruenceTransform(stateTransition, this->covariance()) + this->processNoise();
     }
 
     template<typename QNumberType, std::size_t StateSize, std::size_t MeasurementSize, std::size_t ControlSize>
diff --git a/numerical/filters/active/KalmanFilterBase.hpp b/numerical/filters/active/KalmanFilterBase.hpp
index 4a6968f4..f0636fb5 100644
--- a/numerical/filters/active/KalmanFilterBase.hpp
+++ b/numerical/filters/active/KalmanFilterBase.hpp
@@ -2,6 +2,7 @@
 
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/math/Matrix.hpp"
+#include "numerical/math/MatrixOperations.hpp"
 #include "numerical/solvers/GaussianElimination.hpp"
 
 #if defined(__GNUC__) || defined(__clang__)
@@ -135,7 +136,7 @@ namespace filters
         typename KalmanFilterBase<QNumberType, StateSize, MeasurementSize, ControlSize>::KalmanGain
         KalmanFilterBase<QNumberType, StateSize, MeasurementSize, ControlSize>::ComputeKalmanGain(const MeasurementMatrix& H) const
     {
-        auto S = H * covariance_ * H.Transpose() + measurementNoise_;
+        auto S = math::CongruenceTransform(H, covariance_) + measurementNoise_;
 
         // Solve S Kᵀ = H P instead of explicit K = P Hᵀ S⁻¹ (S is symmetric)
         auto KTranspose = solvers::SolveSystem<QNumberType, MeasurementSize, StateSize>(S, H * covariance_);
@@ -150,7 +151,7 @@ namespace filters
         state_ = state_ + K * innovation;
 
         auto IminusKH = StateMatrix::Identity() - K * H;
-        covariance_ = IminusKH * covariance_ * IminusKH.Transpose() + K * measurementNoise_ * K.Transpose();
+        covariance_ = math::CongruenceTransform(IminusKH, covariance_) + math::CongruenceTransform(K, measurementNoise_);
     }
 
     template<typename QNumberType, std::size_t StateSize, std::size_t MeasurementSize, std::size_t ControlSize>
diff --git a/numerical/filters/active/KalmanSmoother.hpp b/numerical/filters/active/KalmanSmoother.hpp
index 85ec3a13..5998da23 100644
--- a/numerical/filters/active/KalmanSmoother.hpp
+++ b/numerical/filters/active/KalmanSmoother.hpp
@@ -7,6 +7,7 @@
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/math/LinearTimeInvariant.hpp"
 #include "numerical/math/Matrix.hpp"
+#include "numerical/math/MatrixOperations.hpp"
 #include "numerical/solvers/CholeskyDecomposition.hpp"
 #include "numerical/solvers/GaussianElimination.hpp"
 #include <array>
@@ -140,7 +141,7 @@ namespace filters
         for (std::size_t t = 0; t < numSteps; ++t)
         {
             const auto nu = observations[t] - H * predictedMeans_[t];
-            const auto S = H * predictedCovariances_[t] * H.Transpose() + R;
+            const auto S = math::CongruenceTransform(H, predictedCovariances_[t]) + R;
 
             // Kalman gain: K = P H^T S^{-1}  via SolveSystem(S, H*P)^T
             const auto K = solvers::SolveSystem<float, MeasurementSize, StateSize>(
@@ -152,7 +153,7 @@ namespace filters
             // Joseph-form covariance update for numerical PD-preservation
             const auto IminusKH = StateMatrix::Identity() - K * H;
             filteredCovariances_[t] =
-                IminusKH * predictedCovariances_[t] * IminusKH.Transpose() + K * R * K.Transpose();
+                math::CongruenceTransform(IminusKH, predictedCovariances_[t]) + math::CongruenceTransform(K, R);
 
             kalmanGains_[t] = K;
 
@@ -161,7 +162,7 @@ namespace filters
             if (t < numSteps - 1)
             {
                 predictedMeans_[t + 1] = F * filteredMeans_[t];
-                predictedCovariances_[t + 1] = F * filteredCovariances_[t] * F.Transpose() + Q;
+                predictedCovariances_[t + 1] = math::CongruenceTransform(F, filteredCovariances_[t]) + Q;
             }
         }
 
@@ -188,7 +189,7 @@ namespace filters
 
             const auto deltaP = output.smoothedCovariances[t + 1] - predictedCovariances_[t + 1];
             output.smoothedMeans[t] = filteredMeans_[t] + G_t * (output.smoothedMeans[t + 1] - predictedMeans_[t + 1]);
-            output.smoothedCovariances[t] = filteredCovariances_[t] + G_t * deltaP * G_t.Transpose();
+            output.smoothedCovariances[t] = filteredCovariances_[t] + math::CongruenceTransform(G_t, deltaP);
 
             // Lag-1 cross-covariance: P_{t+1, t | T}
             if (t == numSteps - 2)
diff --git a/numerical/filters/active/UnscentedKalmanFilter.hpp b/numerical/filters/active/UnscentedKalmanFilter.hpp
index 39c0cdb3..dac06e68 100644
--- a/numerical/filters/active/UnscentedKalmanFilter.hpp
+++ b/numerical/filters/active/UnscentedKalmanFilter.hpp
@@ -3,6 +3,7 @@
 #include "infra/util/Function.hpp"
 #include "numerical/filters/active/KalmanFilterBase.hpp"
 #include "numerical/math/CompilerOptimizations.hpp"
+#include "numerical/math/MatrixOperations.hpp"
 #include "numerical/solvers/CholeskyDecomposition.hpp"
 #include <algorithm>
 #include <array>
@@ -254,7 +255,7 @@ namespace filters
 
         auto innovation = measurement - predictedMeasurement;
         this->state() = this->state() + K * innovation;
-        this->covariance() = this->covariance() - K * S * K.Transpose();
+        this->covariance() = this->covariance() - math::CongruenceTransform(K, S);
     }
 
 #ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
diff --git a/numerical/math/CMakeLists.txt b/numerical/math/CMakeLists.txt
index 89c1846a..8f923c46 100644
--- a/numerical/math/CMakeLists.txt
+++ b/numerical/math/CMakeLists.txt
@@ -14,10 +14,13 @@ target_sources(numerical.math PRIVATE
     ComplexNumber.hpp
     Cordic.hpp
     Geometry3D.hpp
+    GivensRotation.hpp
+    HouseholderTransform.hpp
     HyperbolicFunctions.hpp
     LinearTimeInvariant.hpp
     Matrix.hpp
     MatrixNorms.hpp
+    MatrixOperations.hpp
     QNumber.hpp
     Quaternion.hpp
     RecursiveBuffer.hpp
@@ -26,6 +29,7 @@ target_sources(numerical.math PRIVATE
     StepResponseMetrics.hpp
     Toeplitz.hpp
     Tolerance.hpp
+    TriangularSolve.hpp
     TrigonometricFunctions.hpp
 )
 
diff --git a/numerical/math/GivensRotation.hpp b/numerical/math/GivensRotation.hpp
new file mode 100644
index 00000000..3a538a06
--- /dev/null
+++ b/numerical/math/GivensRotation.hpp
@@ -0,0 +1,40 @@
+#pragma once
+
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC optimize("O3", "fast-math")
+#endif
+
+#include "numerical/math/CompilerOptimizations.hpp"
+#include <cmath>
+#include <type_traits>
+
+namespace math
+{
+    template<typename T>
+    struct GivensRotation
+    {
+        T c;
+        T s;
+    };
+
+    template<typename T>
+    [[nodiscard]] OPTIMIZE_FOR_SPEED GivensRotation<T> ComputeGivens(T a, T b)
+    {
+        static_assert(std::is_floating_point_v<T>, "ComputeGivens supports floating-point types only");
+
+        T r = std::sqrt(a * a + b * b);
+        if (r < T{ 1e-30f })
+            return { T{ 1 }, T{} };
+
+        return { a / r, b / r };
+    }
+
+    template<typename T>
+    OPTIMIZE_FOR_SPEED void ApplyGivens(const GivensRotation<T>& g, T& x, T& y)
+    {
+        T newX = g.c * x + g.s * y;
+        T newY = -g.s * x + g.c * y;
+        x = newX;
+        y = newY;
+    }
+}
diff --git a/numerical/math/HouseholderTransform.hpp b/numerical/math/HouseholderTransform.hpp
new file mode 100644
index 00000000..4ee2fc2f
--- /dev/null
+++ b/numerical/math/HouseholderTransform.hpp
@@ -0,0 +1,47 @@
+#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 <cmath>
+#include <cstddef>
+#include <type_traits>
+
+namespace math
+{
+    template<typename T, std::size_t N>
+    OPTIMIZE_FOR_SPEED void HouseholderVector(const Vector<T, N>& x, std::size_t start, Vector<T, N>& v, T& beta)
+    {
+        static_assert(std::is_floating_point_v<T>, "HouseholderVector supports floating-point types only");
+
+        T sigma{ T{} };
+        for (std::size_t i = start + 1; i < N; ++i)
+            sigma += x.at(i, 0) * x.at(i, 0);
+
+        for (std::size_t i = 0; i < N; ++i)
+            v.at(i, 0) = (i >= start) ? x.at(i, 0) : T{};
+
+        if (sigma < T{ 1e-30f })
+        {
+            beta = T{};
+            return;
+        }
+
+        T norm = std::sqrt(x.at(start, 0) * x.at(start, 0) + sigma);
+        T v0;
+        if (x.at(start, 0) <= T{})
+            v0 = x.at(start, 0) - norm;
+        else
+            v0 = -sigma / (x.at(start, 0) + norm);
+
+        beta = T{ 2 } * v0 * v0 / (sigma + v0 * v0);
+
+        T invV0 = T{ 1 } / v0;
+        v.at(start, 0) = T{ 1 };
+        for (std::size_t i = start + 1; i < N; ++i)
+            v.at(i, 0) = x.at(i, 0) * invV0;
+    }
+}
diff --git a/numerical/math/Matrix.hpp b/numerical/math/Matrix.hpp
index de64ffc4..01f95723 100644
--- a/numerical/math/Matrix.hpp
+++ b/numerical/math/Matrix.hpp
@@ -285,7 +285,12 @@ namespace math
 
         Matrix result;
         for (size_type i = 0; i < Rows; ++i)
-            result.at(i, i) = T(0.9999f);
+        {
+            if constexpr (std::is_floating_point_v<T>)
+                result.at(i, i) = T(1);
+            else
+                result.at(i, i) = T(0.9999f);
+        }
 
         return result;
     }
diff --git a/numerical/math/MatrixNorms.hpp b/numerical/math/MatrixNorms.hpp
index 628da377..02b53b7e 100644
--- a/numerical/math/MatrixNorms.hpp
+++ b/numerical/math/MatrixNorms.hpp
@@ -55,13 +55,20 @@ namespace math
     }
 
     template<typename T, std::size_t Size>
-    [[nodiscard]] OPTIMIZE_FOR_SPEED T VectorNorm(const Vector<T, Size>& v)
+    [[nodiscard]] OPTIMIZE_FOR_SPEED T DotProduct(const Vector<T, Size>& a, const Vector<T, Size>& b)
     {
         static_assert(std::is_floating_point_v<T>, "MatrixNorms supports floating-point types");
         T sum{};
         for (std::size_t i = 0; i < Size; ++i)
-            sum += v.at(i, 0) * v.at(i, 0);
-        return std::sqrt(sum);
+            sum += a.at(i, 0) * b.at(i, 0);
+        return sum;
+    }
+
+    template<typename T, std::size_t Size>
+    [[nodiscard]] OPTIMIZE_FOR_SPEED T VectorNorm(const Vector<T, Size>& v)
+    {
+        static_assert(std::is_floating_point_v<T>, "MatrixNorms supports floating-point types");
+        return std::sqrt(DotProduct(v, v));
     }
 
     template<typename T, std::size_t Size>
diff --git a/numerical/math/MatrixOperations.hpp b/numerical/math/MatrixOperations.hpp
new file mode 100644
index 00000000..d3c0d148
--- /dev/null
+++ b/numerical/math/MatrixOperations.hpp
@@ -0,0 +1,27 @@
+#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 <cstddef>
+#include <type_traits>
+
+namespace math
+{
+    template<typename T, std::size_t N>
+    [[nodiscard]] OPTIMIZE_FOR_SPEED SquareMatrix<T, N> Symmetrize(const SquareMatrix<T, N>& m)
+    {
+        static_assert(std::is_floating_point_v<T>, "Symmetrize supports floating-point types");
+        return (m + m.Transpose()) * T(0.5f);
+    }
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    [[nodiscard]] OPTIMIZE_FOR_SPEED SquareMatrix<T, Rows> CongruenceTransform(const Matrix<T, Rows, Cols>& a, const SquareMatrix<T, Cols>& m)
+    {
+        static_assert(std::is_floating_point_v<T>, "CongruenceTransform supports floating-point types");
+        return a * m * a.Transpose();
+    }
+}
diff --git a/numerical/math/TriangularSolve.hpp b/numerical/math/TriangularSolve.hpp
new file mode 100644
index 00000000..e2d4d39f
--- /dev/null
+++ b/numerical/math/TriangularSolve.hpp
@@ -0,0 +1,34 @@
+#pragma once
+
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC optimize("O3", "fast-math")
+#endif
+
+#include "infra/util/ReallyAssert.hpp"
+#include "numerical/math/CompilerOptimizations.hpp"
+#include "numerical/math/Matrix.hpp"
+#include "numerical/math/QNumber.hpp"
+#include <cmath>
+#include <cstddef>
+
+namespace math
+{
+    template<typename T, std::size_t N>
+    [[nodiscard]] OPTIMIZE_FOR_SPEED Vector<T, N> SolveUpperTriangular(const Matrix<T, N, N>& r, const Vector<T, N>& c)
+    {
+        Vector<T, N> x{};
+
+        for (std::size_t i = N; i > 0; --i)
+        {
+            std::size_t row = i - 1;
+            T sum = c.at(row, 0);
+            for (std::size_t j = row + 1; j < N; ++j)
+                sum = sum - r.at(row, j) * x.at(j, 0);
+
+            really_assert(std::abs(ToFloat(r.at(row, row))) > 0.0f);
+            x.at(row, 0) = sum / r.at(row, row);
+        }
+
+        return x;
+    }
+}
diff --git a/numerical/math/test/CMakeLists.txt b/numerical/math/test/CMakeLists.txt
index 9c63e2b7..4575b648 100644
--- a/numerical/math/test/CMakeLists.txt
+++ b/numerical/math/test/CMakeLists.txt
@@ -10,8 +10,11 @@ target_link_libraries(numerical.math_test PUBLIC
 target_sources(numerical.math_test PRIVATE
     TestComplexNumber.cpp
     TestCordic.cpp
+    TestGivensRotation.cpp
+    TestHouseholderTransform.cpp
     TestLinearTimeInvariant.cpp
     TestMatrixNorms.cpp
+    TestMatrixOperations.cpp
     TestQNumber.cpp
     TestMatrix.cpp
     TestQuaternion.cpp
@@ -19,4 +22,5 @@ target_sources(numerical.math_test PRIVATE
     TestStatistics.cpp
     TestStepResponseMetrics.cpp
     TestToeplitz.cpp
+    TestTriangularSolve.cpp
 )
diff --git a/numerical/math/test/TestGivensRotation.cpp b/numerical/math/test/TestGivensRotation.cpp
new file mode 100644
index 00000000..b5d18634
--- /dev/null
+++ b/numerical/math/test/TestGivensRotation.cpp
@@ -0,0 +1,43 @@
+#include "numerical/math/GivensRotation.hpp"
+#include "numerical/math/Tolerance.hpp"
+#include <cmath>
+#include <gtest/gtest.h>
+
+namespace
+{
+    class GivensRotationTest : public ::testing::Test
+    {
+    };
+}
+
+TEST_F(GivensRotationTest, ZerosSecondComponent)
+{
+    auto g = math::ComputeGivens(3.0f, 4.0f);
+
+    float x = 3.0f;
+    float y = 4.0f;
+    math::ApplyGivens(g, x, y);
+
+    EXPECT_NEAR(x, 5.0f, math::Tolerance<float>());
+    EXPECT_NEAR(y, 0.0f, math::Tolerance<float>());
+}
+
+TEST_F(GivensRotationTest, PreservesNorm)
+{
+    auto g = math::ComputeGivens(1.0f, 2.0f);
+
+    float x = -3.0f;
+    float y = 7.0f;
+    float before = std::sqrt(x * x + y * y);
+    math::ApplyGivens(g, x, y);
+
+    EXPECT_NEAR(std::sqrt(x * x + y * y), before, math::Tolerance<float>());
+}
+
+TEST_F(GivensRotationTest, DegenerateInputIsIdentity)
+{
+    auto g = math::ComputeGivens(0.0f, 0.0f);
+
+    EXPECT_NEAR(g.c, 1.0f, math::Tolerance<float>());
+    EXPECT_NEAR(g.s, 0.0f, math::Tolerance<float>());
+}
diff --git a/numerical/math/test/TestHouseholderTransform.cpp b/numerical/math/test/TestHouseholderTransform.cpp
new file mode 100644
index 00000000..1fade4f1
--- /dev/null
+++ b/numerical/math/test/TestHouseholderTransform.cpp
@@ -0,0 +1,64 @@
+#include "numerical/math/HouseholderTransform.hpp"
+#include "numerical/math/Tolerance.hpp"
+#include <cmath>
+#include <gtest/gtest.h>
+
+namespace
+{
+    class HouseholderTransformTest : public ::testing::Test
+    {
+    protected:
+        math::Vector<float, 4> Reflect(const math::Vector<float, 4>& x, const math::Vector<float, 4>& v, float beta)
+        {
+            float dot = 0.0f;
+            for (std::size_t i = 0; i < 4; ++i)
+                dot += v.at(i, 0) * x.at(i, 0);
+
+            math::Vector<float, 4> result;
+            for (std::size_t i = 0; i < 4; ++i)
+                result.at(i, 0) = x.at(i, 0) - beta * dot * v.at(i, 0);
+            return result;
+        }
+    };
+}
+
+TEST_F(HouseholderTransformTest, ZerosEntriesBelowPivot)
+{
+    math::Vector<float, 4> x{ { 4.0f }, { 3.0f }, { 0.0f }, { 0.0f } };
+    math::Vector<float, 4> v;
+    float beta{};
+
+    math::HouseholderVector(x, 0, v, beta);
+    auto reflected = Reflect(x, v, beta);
+
+    EXPECT_NEAR(reflected.at(1, 0), 0.0f, math::Tolerance<float>());
+    EXPECT_NEAR(reflected.at(2, 0), 0.0f, math::Tolerance<float>());
+    EXPECT_NEAR(reflected.at(3, 0), 0.0f, math::Tolerance<float>());
+    EXPECT_NEAR(std::abs(reflected.at(0, 0)), 5.0f, math::Tolerance<float>());
+}
+
+TEST_F(HouseholderTransformTest, RespectsStartOffset)
+{
+    math::Vector<float, 4> x{ { 9.0f }, { 0.0f }, { 3.0f }, { 4.0f } };
+    math::Vector<float, 4> v;
+    float beta{};
+
+    math::HouseholderVector(x, 1, v, beta);
+    auto reflected = Reflect(x, v, beta);
+
+    EXPECT_NEAR(reflected.at(0, 0), 9.0f, math::Tolerance<float>());
+    EXPECT_NEAR(std::abs(reflected.at(1, 0)), 5.0f, math::Tolerance<float>());
+    EXPECT_NEAR(reflected.at(2, 0), 0.0f, math::Tolerance<float>());
+    EXPECT_NEAR(reflected.at(3, 0), 0.0f, math::Tolerance<float>());
+}
+
+TEST_F(HouseholderTransformTest, ZeroSubvectorYieldsZeroBeta)
+{
+    math::Vector<float, 4> x{ { 2.0f }, { 0.0f }, { 0.0f }, { 0.0f } };
+    math::Vector<float, 4> v;
+    float beta{ 1.0f };
+
+    math::HouseholderVector(x, 0, v, beta);
+
+    EXPECT_NEAR(beta, 0.0f, math::Tolerance<float>());
+}
diff --git a/numerical/math/test/TestMatrix.cpp b/numerical/math/test/TestMatrix.cpp
index 4972559e..0fc7fe89 100644
--- a/numerical/math/test/TestMatrix.cpp
+++ b/numerical/math/test/TestMatrix.cpp
@@ -139,8 +139,18 @@ TYPED_TEST(MatrixTest, Identity)
 {
     auto identity = TestFixture::MatrixType::Identity();
 
-    auto expected = this->MakeMatrix(0.9999f, 0.0f, 0.0f, 0.9999f);
-    EXPECT_TRUE(AreMatricesNear(identity, expected));
+    float expectedDiag = std::is_same_v<TypeParam, float> ? 1.0f : 0.9999f;
+    for (size_t i = 0; i < 2; ++i)
+    {
+        for (size_t j = 0; j < 2; ++j)
+        {
+            float expected = (i == j) ? expectedDiag : 0.0f;
+            if constexpr (std::is_same_v<TypeParam, float>)
+                EXPECT_NEAR(identity.at(i, j), expected, 1e-4f);
+            else
+                EXPECT_NEAR(identity.at(i, j).ToFloat(), expected, 1e-4f);
+        }
+    }
 }
 
 TYPED_TEST(MatrixTest, RangeLimits)
diff --git a/numerical/math/test/TestMatrixNorms.cpp b/numerical/math/test/TestMatrixNorms.cpp
index 80be27a7..c609f7a9 100644
--- a/numerical/math/test/TestMatrixNorms.cpp
+++ b/numerical/math/test/TestMatrixNorms.cpp
@@ -39,6 +39,20 @@ TEST_F(MatrixNormsTest, VectorNorm)
     EXPECT_NEAR(result, 5.0f, math::Tolerance<float>());
 }
 
+TEST_F(MatrixNormsTest, DotProduct)
+{
+    math::Vector<float, 2> w{ { 1.0f }, { 2.0f } };
+    float result = math::DotProduct(v, w);
+    EXPECT_NEAR(result, 11.0f, math::Tolerance<float>());
+}
+
+TEST_F(MatrixNormsTest, DotProductWithSelfEqualsNormSquared)
+{
+    float dot = math::DotProduct(v, v);
+    float norm = math::VectorNorm(v);
+    EXPECT_NEAR(dot, norm * norm, math::Tolerance<float>());
+}
+
 TEST_F(MatrixNormsTest, NormalizeUnit)
 {
     auto result = math::Normalize(v);
diff --git a/numerical/math/test/TestMatrixOperations.cpp b/numerical/math/test/TestMatrixOperations.cpp
new file mode 100644
index 00000000..482c7761
--- /dev/null
+++ b/numerical/math/test/TestMatrixOperations.cpp
@@ -0,0 +1,71 @@
+#include "numerical/math/MatrixOperations.hpp"
+#include "numerical/math/Tolerance.hpp"
+#include <gtest/gtest.h>
+
+namespace
+{
+    class MatrixOperationsTest : public ::testing::Test
+    {
+    protected:
+        math::SquareMatrix<float, 2> asymmetric{
+            { 1.0f, 3.0f },
+            { -1.0f, 2.0f }
+        };
+    };
+}
+
+TEST_F(MatrixOperationsTest, SymmetrizeAveragesOffDiagonals)
+{
+    auto s = math::Symmetrize(asymmetric);
+
+    EXPECT_NEAR(s.at(0, 0), 1.0f, math::Tolerance<float>());
+    EXPECT_NEAR(s.at(1, 1), 2.0f, math::Tolerance<float>());
+    EXPECT_NEAR(s.at(0, 1), 1.0f, math::Tolerance<float>());
+    EXPECT_NEAR(s.at(1, 0), 1.0f, math::Tolerance<float>());
+}
+
+TEST_F(MatrixOperationsTest, SymmetrizeIsIdempotentOnSymmetricInput)
+{
+    auto once = math::Symmetrize(asymmetric);
+    auto twice = math::Symmetrize(once);
+
+    for (std::size_t i = 0; i < 2; ++i)
+        for (std::size_t j = 0; j < 2; ++j)
+            EXPECT_NEAR(twice.at(i, j), once.at(i, j), math::Tolerance<float>());
+}
+
+TEST_F(MatrixOperationsTest, CongruenceTransformMatchesExplicitProduct)
+{
+    math::Matrix<float, 2, 3> a{
+        { 1.0f, 2.0f, 3.0f },
+        { 4.0f, 5.0f, 6.0f }
+    };
+    math::SquareMatrix<float, 3> m{
+        { 2.0f, 0.0f, 1.0f },
+        { 0.0f, 3.0f, 0.0f },
+        { 1.0f, 0.0f, 4.0f }
+    };
+
+    auto result = math::CongruenceTransform(a, m);
+    auto expected = a * m * a.Transpose();
+
+    for (std::size_t i = 0; i < 2; ++i)
+        for (std::size_t j = 0; j < 2; ++j)
+            EXPECT_NEAR(result.at(i, j), expected.at(i, j), math::Tolerance<float>());
+}
+
+TEST_F(MatrixOperationsTest, CongruenceTransformPreservesSymmetry)
+{
+    math::Matrix<float, 2, 2> a{
+        { 1.0f, 2.0f },
+        { -3.0f, 0.5f }
+    };
+    math::SquareMatrix<float, 2> symmetric{
+        { 4.0f, 1.0f },
+        { 1.0f, 2.0f }
+    };
+
+    auto result = math::CongruenceTransform(a, symmetric);
+
+    EXPECT_NEAR(result.at(0, 1), result.at(1, 0), math::Tolerance<float>());
+}
diff --git a/numerical/math/test/TestTriangularSolve.cpp b/numerical/math/test/TestTriangularSolve.cpp
new file mode 100644
index 00000000..7e8288b7
--- /dev/null
+++ b/numerical/math/test/TestTriangularSolve.cpp
@@ -0,0 +1,45 @@
+#include "numerical/math/Tolerance.hpp"
+#include "numerical/math/TriangularSolve.hpp"
+#include <gtest/gtest.h>
+
+namespace
+{
+    class TriangularSolveTest : public ::testing::Test
+    {
+    protected:
+        math::Matrix<float, 3, 3> r{
+            { 2.0f, -1.0f, 3.0f },
+            { 0.0f, 4.0f, 1.0f },
+            { 0.0f, 0.0f, 5.0f }
+        };
+        math::Vector<float, 3> x{ { 1.0f }, { -2.0f }, { 3.0f } };
+    };
+}
+
+TEST_F(TriangularSolveTest, SolvesUpperTriangularSystem)
+{
+    math::Vector<float, 3> c;
+    for (std::size_t i = 0; i < 3; ++i)
+    {
+        float sum = 0.0f;
+        for (std::size_t j = 0; j < 3; ++j)
+            sum += r.at(i, j) * x.at(j, 0);
+        c.at(i, 0) = sum;
+    }
+
+    auto solution = math::SolveUpperTriangular(r, c);
+
+    for (std::size_t i = 0; i < 3; ++i)
+        EXPECT_NEAR(solution.at(i, 0), x.at(i, 0), math::Tolerance<float>());
+}
+
+TEST_F(TriangularSolveTest, IdentityReturnsRightHandSide)
+{
+    auto identity = math::SquareMatrix<float, 3>::Identity();
+    math::Vector<float, 3> b{ { 7.0f }, { -3.0f }, { 2.0f } };
+
+    auto solution = math::SolveUpperTriangular(identity, b);
+
+    for (std::size_t i = 0; i < 3; ++i)
+        EXPECT_NEAR(solution.at(i, 0), b.at(i, 0), math::Tolerance<float>());
+}
diff --git a/numerical/solvers/CMakeLists.txt b/numerical/solvers/CMakeLists.txt
index 7c08d33e..697a8007 100644
--- a/numerical/solvers/CMakeLists.txt
+++ b/numerical/solvers/CMakeLists.txt
@@ -19,6 +19,7 @@ target_sources(numerical.solver PRIVATE
     GaussianElimination.hpp
     LevinsonDurbin.hpp
     OdeSystem.hpp
+    QrDecomposition.hpp
     RungeKuttaIntegrators.hpp
     Solver.hpp
     SpectralRadius.hpp
@@ -28,6 +29,7 @@ numerical_add_coverage_sources(numerical.solver
     DiscreteAlgebraicRiccatiEquation.cpp
     DurandKerner.cpp
     GaussianElimination.cpp
+    QrDecomposition.cpp
     RungeKuttaIntegrators.cpp
     SpectralRadius.cpp
 )
diff --git a/numerical/solvers/GaussianElimination.hpp b/numerical/solvers/GaussianElimination.hpp
index 66c7fee1..f5769549 100644
--- a/numerical/solvers/GaussianElimination.hpp
+++ b/numerical/solvers/GaussianElimination.hpp
@@ -8,6 +8,7 @@
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/math/Matrix.hpp"
 #include "numerical/math/QNumber.hpp"
+#include "numerical/math/TriangularSolve.hpp"
 #include "numerical/solvers/Solver.hpp"
 #include <cmath>
 
@@ -29,7 +30,6 @@ namespace solvers
         std::size_t FindPivotRow(const InputMatrix& matrix, std::size_t col) const;
         void SwapRows(InputMatrix& matrix, SolutionVector& vector, std::size_t row1, std::size_t row2) const;
         void EliminateBelow(InputMatrix& matrix, SolutionVector& vector, std::size_t col) const;
-        SolutionVector BackSubstitute(const InputMatrix& matrix, const SolutionVector& vector) const;
     };
 
     template<typename T, std::size_t N, std::size_t Cols>
@@ -86,27 +86,6 @@ namespace solvers
         }
     }
 
-    template<typename T, std::size_t N>
-    typename GaussianElimination<T, N>::SolutionVector
-    GaussianElimination<T, N>::BackSubstitute(const InputMatrix& matrix, const SolutionVector& vector) const
-    {
-        SolutionVector x;
-
-        for (std::size_t i = N; i > 0; --i)
-        {
-            std::size_t row = i - 1;
-            T sum = vector.at(row, 0);
-
-            for (std::size_t j = row + 1; j < N; ++j)
-                sum = sum - matrix.at(row, j) * x.at(j, 0);
-
-            really_assert(std::abs(math::ToFloat(matrix.at(row, row))) > 0.0f);
-            x.at(row, 0) = sum / matrix.at(row, row);
-        }
-
-        return x;
-    }
-
     template<typename T, std::size_t N>
     OPTIMIZE_FOR_SPEED
         typename GaussianElimination<T, N>::SolutionVector
@@ -125,7 +104,7 @@ namespace solvers
             EliminateBelow(augA, augB, col);
         }
 
-        return BackSubstitute(augA, augB);
+        return math::SolveUpperTriangular(augA, augB);
     }
 
     template<typename T, std::size_t N, std::size_t Cols>
diff --git a/numerical/solvers/QrDecomposition.cpp b/numerical/solvers/QrDecomposition.cpp
new file mode 100644
index 00000000..c9ccce3c
--- /dev/null
+++ b/numerical/solvers/QrDecomposition.cpp
@@ -0,0 +1,7 @@
+#include "numerical/solvers/QrDecomposition.hpp"
+
+namespace solvers
+{
+    template class QrDecomposition<float, 4, 3>;
+    template class QrDecomposition<float, 3, 3>;
+}
diff --git a/numerical/solvers/QrDecomposition.hpp b/numerical/solvers/QrDecomposition.hpp
new file mode 100644
index 00000000..92715f53
--- /dev/null
+++ b/numerical/solvers/QrDecomposition.hpp
@@ -0,0 +1,192 @@
+#pragma once
+
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC optimize("O3", "fast-math")
+#endif
+
+#include "numerical/math/CompilerOptimizations.hpp"
+#include "numerical/math/GivensRotation.hpp"
+#include "numerical/math/HouseholderTransform.hpp"
+#include "numerical/math/Matrix.hpp"
+#include "numerical/math/TriangularSolve.hpp"
+#include <cmath>
+#include <cstddef>
+#include <type_traits>
+
+namespace solvers
+{
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    class QrDecomposition
+    {
+        static_assert(std::is_floating_point_v<T>, "QrDecomposition supports floating-point types only");
+        static_assert(Rows >= Cols, "QrDecomposition requires Rows >= Cols");
+
+    public:
+        QrDecomposition() = default;
+
+        OPTIMIZE_FOR_SPEED bool Decompose(const math::Matrix<T, Rows, Cols>& a);
+        math::Matrix<T, Rows, Cols> Q() const;
+        math::Matrix<T, Cols, Cols> R() const;
+        OPTIMIZE_FOR_SPEED math::Vector<T, Cols> SolveLeastSquares(const math::Vector<T, Rows>& b) const;
+        void ApplyQtranspose(math::Vector<T, Rows>& v) const;
+        void GivensUpdateRow(const math::Matrix<T, 1, Cols>& newRow);
+
+    private:
+        static void ApplyReflector(math::Matrix<T, Rows, Cols>& block, const math::Vector<T, Rows>& v, T beta, std::size_t start);
+
+        math::Matrix<T, Rows, Cols> qr{};
+        math::Vector<T, Cols> betas{};
+        bool factored{ false };
+    };
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    OPTIMIZE_FOR_SPEED void QrDecomposition<T, Rows, Cols>::ApplyReflector(
+        math::Matrix<T, Rows, Cols>& block, const math::Vector<T, Rows>& v, T beta, std::size_t start)
+    {
+        for (std::size_t j = start; j < Cols; ++j)
+        {
+            T dot{ T{} };
+            for (std::size_t i = start; i < Rows; ++i)
+                dot += v.at(i, 0) * block.at(i, j);
+
+            T scale = beta * dot;
+            for (std::size_t i = start; i < Rows; ++i)
+                block.at(i, j) -= scale * v.at(i, 0);
+        }
+    }
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    OPTIMIZE_FOR_SPEED bool QrDecomposition<T, Rows, Cols>::Decompose(const math::Matrix<T, Rows, Cols>& a)
+    {
+        qr = a;
+        betas = math::Vector<T, Cols>{};
+        factored = false;
+
+        math::Vector<T, Rows> colVec{};
+        math::Vector<T, Rows> v{};
+
+        T maxPivot{ T{} };
+
+        for (std::size_t k = 0; k < Cols; ++k)
+        {
+            for (std::size_t i = 0; i < Rows; ++i)
+                colVec.at(i, 0) = (i >= k) ? qr.at(i, k) : T{};
+
+            T beta{};
+            math::HouseholderVector(colVec, k, v, beta);
+            betas.at(k, 0) = beta;
+
+            if (beta != T{})
+            {
+                ApplyReflector(qr, v, beta, k);
+
+                for (std::size_t i = k + 1; i < Rows; ++i)
+                    qr.at(i, k) = v.at(i, 0);
+            }
+
+            T pivot = std::abs(qr.at(k, k));
+            if (pivot > maxPivot)
+                maxPivot = pivot;
+
+            if (pivot <= maxPivot * T{ 1e-6f })
+                return false;
+        }
+
+        factored = true;
+        return true;
+    }
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    math::Matrix<T, Rows, Cols> QrDecomposition<T, Rows, Cols>::Q() const
+    {
+        math::Matrix<T, Rows, Cols> result{};
+        for (std::size_t j = 0; j < Cols; ++j)
+            result.at(j, j) = T{ 1 };
+
+        math::Vector<T, Rows> v{};
+
+        for (std::size_t k = Cols; k > 0; --k)
+        {
+            std::size_t col = k - 1;
+            T beta = betas.at(col, 0);
+            if (beta == T{})
+                continue;
+
+            v.at(col, 0) = T{ 1 };
+            for (std::size_t i = col + 1; i < Rows; ++i)
+                v.at(i, 0) = qr.at(i, col);
+
+            ApplyReflector(result, v, beta, col);
+        }
+
+        return result;
+    }
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    math::Matrix<T, Cols, Cols> QrDecomposition<T, Rows, Cols>::R() const
+    {
+        math::Matrix<T, Cols, Cols> r{};
+        for (std::size_t i = 0; i < Cols; ++i)
+            for (std::size_t j = i; j < Cols; ++j)
+                r.at(i, j) = qr.at(i, j);
+        return r;
+    }
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    void QrDecomposition<T, Rows, Cols>::ApplyQtranspose(math::Vector<T, Rows>& v) const
+    {
+        math::Vector<T, Rows> hv{};
+
+        for (std::size_t k = 0; k < Cols; ++k)
+        {
+            T beta = betas.at(k, 0);
+            if (beta == T{})
+                continue;
+
+            hv.at(k, 0) = T{ 1 };
+            for (std::size_t i = k + 1; i < Rows; ++i)
+                hv.at(i, 0) = qr.at(i, k);
+
+            T dot{ T{} };
+            for (std::size_t i = k; i < Rows; ++i)
+                dot += hv.at(i, 0) * v.at(i, 0);
+
+            T scale = beta * dot;
+            for (std::size_t i = k; i < Rows; ++i)
+                v.at(i, 0) -= scale * hv.at(i, 0);
+        }
+    }
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    OPTIMIZE_FOR_SPEED math::Vector<T, Cols>
+    QrDecomposition<T, Rows, Cols>::SolveLeastSquares(const math::Vector<T, Rows>& b) const
+    {
+        math::Vector<T, Rows> c = b;
+        ApplyQtranspose(c);
+
+        math::Vector<T, Cols> c1{};
+        for (std::size_t i = 0; i < Cols; ++i)
+            c1.at(i, 0) = c.at(i, 0);
+
+        return math::SolveUpperTriangular(R(), c1);
+    }
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    void QrDecomposition<T, Rows, Cols>::GivensUpdateRow(const math::Matrix<T, 1, Cols>& newRow)
+    {
+        math::Matrix<T, 1, Cols> row = newRow;
+
+        for (std::size_t k = 0; k < Cols; ++k)
+        {
+            math::GivensRotation<T> g = math::ComputeGivens(qr.at(k, k), row.at(0, k));
+
+            for (std::size_t j = k; j < Cols; ++j)
+                math::ApplyGivens(g, qr.at(k, j), row.at(0, j));
+        }
+    }
+
+#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
+    extern template class QrDecomposition<float, 4, 3>;
+    extern template class QrDecomposition<float, 3, 3>;
+#endif
+}
diff --git a/numerical/solvers/test/CMakeLists.txt b/numerical/solvers/test/CMakeLists.txt
index 5fe2d859..dc8d0182 100644
--- a/numerical/solvers/test/CMakeLists.txt
+++ b/numerical/solvers/test/CMakeLists.txt
@@ -13,6 +13,7 @@ target_sources(numerical.solvers_test PRIVATE
     TestDurandKerner.cpp
     TestGaussianElimination.cpp
     TestLevinsonDurbin.cpp
+    TestQrDecomposition.cpp
     TestRungeKuttaIntegrators.cpp
     TestSpectralRadius.cpp
 )
diff --git a/numerical/solvers/test/TestQrDecomposition.cpp b/numerical/solvers/test/TestQrDecomposition.cpp
new file mode 100644
index 00000000..00584eaf
--- /dev/null
+++ b/numerical/solvers/test/TestQrDecomposition.cpp
@@ -0,0 +1,188 @@
+#include "numerical/math/Tolerance.hpp"
+#include "numerical/solvers/QrDecomposition.hpp"
+#include <gtest/gtest.h>
+
+namespace
+{
+    class TestQrDecomposition : public ::testing::Test
+    {
+    protected:
+        solvers::QrDecomposition<float, 4, 3> qr{};
+        solvers::QrDecomposition<float, 3, 3> qrSquare{};
+    };
+}
+
+TEST_F(TestQrDecomposition, reconstructs_A)
+{
+    math::Matrix<float, 4, 3> a{
+        { 1.0f, 2.0f, 3.0f },
+        { 4.0f, 5.0f, 6.0f },
+        { 7.0f, 8.0f, 10.0f },
+        { 1.0f, 0.0f, 2.0f }
+    };
+
+    qr.Decompose(a);
+    auto reconstructed = qr.Q() * qr.R();
+
+    for (std::size_t i = 0; i < 4; ++i)
+        for (std::size_t j = 0; j < 3; ++j)
+            EXPECT_NEAR(reconstructed.at(i, j), a.at(i, j), 1e-5f);
+}
+
+TEST_F(TestQrDecomposition, Q_columns_are_orthonormal)
+{
+    math::Matrix<float, 4, 3> a{
+        { 1.0f, 2.0f, 3.0f },
+        { 4.0f, 5.0f, 6.0f },
+        { 7.0f, 8.0f, 10.0f },
+        { 1.0f, 0.0f, 2.0f }
+    };
+
+    qr.Decompose(a);
+    auto qMat = qr.Q();
+    auto qtq = qMat.Transpose() * qMat;
+
+    for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t j = 0; j < 3; ++j)
+            EXPECT_NEAR(qtq.at(i, j), (i == j) ? 1.0f : 0.0f, 1e-5f);
+}
+
+TEST_F(TestQrDecomposition, R_is_upper_triangular)
+{
+    math::Matrix<float, 4, 3> a{
+        { 1.0f, 2.0f, 3.0f },
+        { 4.0f, 5.0f, 6.0f },
+        { 7.0f, 8.0f, 10.0f },
+        { 1.0f, 0.0f, 2.0f }
+    };
+
+    qr.Decompose(a);
+    auto r = qr.R();
+
+    for (std::size_t i = 1; i < 3; ++i)
+        for (std::size_t j = 0; j < i; ++j)
+            EXPECT_NEAR(r.at(i, j), 0.0f, 1e-5f);
+}
+
+TEST_F(TestQrDecomposition, solves_overdetermined_least_squares)
+{
+    math::Matrix<float, 4, 3> a{
+        { 1.0f, 0.0f, 0.0f },
+        { 0.0f, 1.0f, 0.0f },
+        { 0.0f, 0.0f, 1.0f },
+        { 1.0f, 1.0f, 1.0f }
+    };
+
+    math::Vector<float, 4> b{ { 1.0f }, { 2.0f }, { 3.0f }, { 6.5f } };
+
+    qr.Decompose(a);
+    auto x = qr.SolveLeastSquares(b);
+
+    EXPECT_NEAR(x.at(0, 0), 1.125f, 1e-4f);
+    EXPECT_NEAR(x.at(1, 0), 2.125f, 1e-4f);
+    EXPECT_NEAR(x.at(2, 0), 3.125f, 1e-4f);
+}
+
+TEST_F(TestQrDecomposition, matches_known_3x3)
+{
+    math::Matrix<float, 3, 3> a{
+        { 12.0f, -51.0f, 4.0f },
+        { 6.0f, 167.0f, -68.0f },
+        { -4.0f, 24.0f, -41.0f }
+    };
+
+    qrSquare.Decompose(a);
+    auto r = qrSquare.R();
+
+    EXPECT_NEAR(std::abs(r.at(0, 0)), 14.0f, 1e-3f);
+    EXPECT_NEAR(std::abs(r.at(1, 1)), 175.0f, 1e-3f);
+    EXPECT_NEAR(std::abs(r.at(2, 2)), 35.0f, 1e-3f);
+    EXPECT_NEAR(r.at(1, 0), 0.0f, 1e-3f);
+    EXPECT_NEAR(r.at(2, 0), 0.0f, 1e-3f);
+    EXPECT_NEAR(r.at(2, 1), 0.0f, 1e-3f);
+}
+
+TEST_F(TestQrDecomposition, apply_qtranspose_equals_explicit)
+{
+    math::Matrix<float, 4, 3> a{
+        { 1.0f, 2.0f, 3.0f },
+        { 4.0f, 5.0f, 6.0f },
+        { 7.0f, 8.0f, 10.0f },
+        { 1.0f, 0.0f, 2.0f }
+    };
+
+    qr.Decompose(a);
+
+    math::Vector<float, 4> b{ { 1.0f }, { 2.0f }, { 3.0f }, { 4.0f } };
+    math::Vector<float, 4> bCopy = b;
+
+    qr.ApplyQtranspose(bCopy);
+
+    auto qMat = qr.Q();
+    auto explicit_qt_b = qMat.Transpose() * b;
+
+    for (std::size_t i = 0; i < 3; ++i)
+        EXPECT_NEAR(bCopy.at(i, 0), explicit_qt_b.at(i, 0), 1e-5f);
+}
+
+TEST_F(TestQrDecomposition, rank_deficient_returns_false)
+{
+    math::Matrix<float, 4, 3> a{
+        { 1.0f, 2.0f, 3.0f },
+        { 2.0f, 4.0f, 6.0f },
+        { 3.0f, 6.0f, 9.0f },
+        { 4.0f, 8.0f, 12.0f }
+    };
+
+    bool ok = qr.Decompose(a);
+
+    EXPECT_FALSE(ok);
+}
+
+TEST_F(TestQrDecomposition, givens_update_adds_row)
+{
+    math::Matrix<float, 3, 3> a3{
+        { 1.0f, 2.0f, 3.0f },
+        { 4.0f, 5.0f, 6.0f },
+        { 7.0f, 8.0f, 10.0f }
+    };
+
+    math::Matrix<float, 1, 3> newRow{
+        { 1.0f, 0.0f, 2.0f }
+    };
+
+    math::Matrix<float, 4, 3> a4{
+        { 1.0f, 2.0f, 3.0f },
+        { 4.0f, 5.0f, 6.0f },
+        { 7.0f, 8.0f, 10.0f },
+        { 1.0f, 0.0f, 2.0f }
+    };
+
+    qrSquare.Decompose(a3);
+    qrSquare.GivensUpdateRow(newRow);
+    auto rUpdated = qrSquare.R();
+
+    qr.Decompose(a4);
+    auto rFull = qr.R();
+
+    for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t j = i; j < 3; ++j)
+            EXPECT_NEAR(std::abs(rUpdated.at(i, j)), std::abs(rFull.at(i, j)), 1e-3f);
+}
+
+TEST_F(TestQrDecomposition, identity_factors_to_identity)
+{
+    auto a = math::SquareMatrix<float, 3>::Identity();
+
+    qrSquare.Decompose(a);
+    auto q = qrSquare.Q();
+    auto r = qrSquare.R();
+
+    for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t j = 0; j < 3; ++j)
+            EXPECT_NEAR(q.at(i, j), a.at(i, j), math::Tolerance<float>());
+
+    for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t j = 0; j < 3; ++j)
+            EXPECT_NEAR(r.at(i, j), a.at(i, j), math::Tolerance<float>());
+}
diff --git a/roadmap/solvers/QrDecomposition/explanation.md b/roadmap/solvers/QrDecomposition/explanation.md
deleted file mode 100644
index 08bf2bc5..00000000
--- a/roadmap/solvers/QrDecomposition/explanation.md
+++ /dev/null
@@ -1,29 +0,0 @@
-# QR Decomposition — Overview
-
-## What it is
-A factorization `A = Q·R` where `Q` has orthonormal columns and `R` is upper triangular. It is the
-numerically robust way to solve least-squares problems and the structural backbone of the eigenvalue,
-SVD, and square-root filtering algorithms.
-
-## Why it matters (embedded)
-Fitting an overdetermined model (more measurements than unknowns) is everywhere — calibration,
-sensor fusion, system identification. Solving it through `QR` avoids forming `AᵀA`, whose condition
-number is the *square* of `A`'s, so it preserves precision on reduced-word-length hardware.
-
-## How it works (intuition)
-A **Householder reflector** mirrors a whole column onto a coordinate axis, zeroing everything below
-the diagonal in one shot; applying `Cols` reflectors triangularizes `A`, and the accumulated
-reflections form the orthogonal `Q`. **Givens rotations** achieve the same by zeroing one entry at a
-time, which is ideal when a single new row streams in and only a rank-1 update is needed.
-
-## Key parameters
-- **Shape (`Rows ≥ Cols`)** — tall matrices for least squares; square for general use.
-- **Reflector vs rotation** — Householder for dense factorization, Givens for streaming updates.
-
-## Reference
-A. S. Householder, "Unitary Triangularization of a Nonsymmetric Matrix," *JACM*, 5(4), 1958;
-Golub & Van Loan, *Matrix Computations*, Ch. 5.
-
-## See also
-`LuDecomposition` (general solve), `CholeskyDecomposition` (SPD), `SingularValueDecomposition`,
-`JacobiEigenSolver`.
diff --git a/roadmap/solvers/QrDecomposition/implementation.md b/roadmap/solvers/QrDecomposition/implementation.md
deleted file mode 100644
index 24e9d62a..00000000
--- a/roadmap/solvers/QrDecomposition/implementation.md
+++ /dev/null
@@ -1,77 +0,0 @@
-# QR Decomposition (Householder / Givens) — Implementation Pseudocode
-
-> Roadmap ref: #27 (Tier 4) · Target: `numerical/solvers` · Namespace `solvers` · Type: `float` (templated on `T`, instantiated for `float` only)
-
-## Data structures
-
-```
-template<typename T, std::size_t Rows, std::size_t Cols>   # static_assert(std::is_floating_point_v<T>); instantiated for float
-class QrDecomposition:
-    Matrix<T, Rows, Cols> qr        # reflectors below diagonal, R on/above diagonal
-    Vector<T, Cols>       betas     # Householder scaling per column
-    static_assert Rows >= Cols
-```
-
-## Interface
-
-```
-QrDecomposition()
-bool Decompose(const Matrix<T,Rows,Cols>& A)      # hot path; false if rank-deficient
-Matrix<T,Rows,Cols> Q() const                     # thin Q (accumulated on request)
-Matrix<T,Cols,Cols> R() const                     # upper triangular
-Vector<T,Cols> SolveLeastSquares(const Vector<T,Rows>& b) const
-void ApplyQtranspose(Vector<T,Rows>& v) const     # v ← Qᵀ v  (no explicit Q)
-```
-
-## Algorithm (pseudocode)
-
-```
-function Decompose(A):                                 # OPTIMIZE_FOR_SPEED
-    qr = A
-    for k in 0 .. Cols-1:
-        x       = qr[k.., k]                           # sub-column below diagonal
-        v, beta = Householder(x)                       # reflector zeroing x below x[0]
-        for j in k .. Cols-1:                          # apply H = I − beta·v·vᵀ to trailing block
-            qr[k.., j] -= beta · v · (vᵀ · qr[k.., j])
-        store v (below diag), betas[k]
-        if |qr[k,k]| < eps: return false               # rank-deficient
-    return true
-
-function SolveLeastSquares(b):                          # OPTIMIZE_FOR_SPEED
-    c = ApplyQtranspose(b)                              # c = Qᵀ b
-    return BackSubstitute(R, c[0..Cols-1])             # R x = c₁  (upper-triangular solve)
-
-# Givens option (streaming): rotate one new row into R with Cols plane rotations,
-# zeroing each sub-diagonal entry — an O(Cols²) rank-1 update, no full re-factor.
-```
-
-## Complexity & memory
-
-- Householder factorization: `O(Rows·Cols²)` (≈ `2·Rows·Cols² − (2/3)Cols³`).
-- Accumulating `Q` is optional — `ApplyQtranspose` avoids forming it for a solve.
-- Memory: `O(Rows·Cols)` for the packed factors — static, no heap.
-
-## Numerical / embedded notes
-
-- **Householder is backward stable** and preferred for dense `A`; **Givens** wins for sparse or
-  *streaming* updates (adding one row/column) because each rotation touches only two rows.
-- Least squares via QR (`R x = Qᵀb`) is far better conditioned than the normal equations
-  `AᵀA x = Aᵀb` (whose condition number is squared).
-- Foundational: reused by `SingularValueDecomposition`, the symmetric-eigen building blocks, and the
-  square-root Kalman filter; reuse the triangular back-substitution shared with `GaussianElimination`.
-- Fix the sign of each Householder pivot (`−sign(x₀)·‖x‖`) to avoid cancellation.
-- Float-only: `static_assert(std::is_floating_point_v<T>)`; the generic `T` signature keeps a
-  `Q15`/`Q31` specialisation cheap to add later.
-
-## Deployment
-
-- Header: `numerical/solvers/QrDecomposition.hpp` — `#pragma once` →
-  `#pragma GCC optimize("O3","fast-math")`, `OPTIMIZE_FOR_SPEED` on `Decompose`/`SolveLeastSquares`, and
-  `extern template class QrDecomposition<float, 4, 3>;` under `#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD`.
-- Coverage: `numerical/solvers/QrDecomposition.cpp` →
-  `template class QrDecomposition<float, 4, 3>;` and `<float, 3, 3>` (the tested shapes).
-- Test: `numerical/solvers/test/TestQrDecomposition.cpp`
-- Doc: `doc/solvers/QrDecomposition.md`
-- CMake: `.hpp` → `target_sources`; `.cpp` → `numerical_add_coverage_sources`;
-  `TestQrDecomposition.cpp` → the `_test` target.
-- Generic pattern: see `roadmap/DEPLOYMENT.md`.
diff --git a/roadmap/solvers/QrDecomposition/tests.md b/roadmap/solvers/QrDecomposition/tests.md
deleted file mode 100644
index 50b4350d..00000000
--- a/roadmap/solvers/QrDecomposition/tests.md
+++ /dev/null
@@ -1,60 +0,0 @@
-# QR Decomposition — Unit Test Plan (Pseudocode)
-
-> GoogleTest · `TEST_F` (`float`) · `StrictMock` only · no heap.
-
-## Fixture
-
-```
-class TestQrDecomposition : public ::testing::Test:
-    QrDecomposition<float, 4, 3> qr        # tall, overdetermined
-    QrDecomposition<float, 3, 3> qrSquare
-# each case below is a TEST_F(TestQrDecomposition, <name>)
-```
-
-## Test cases (Arrange / Act / Assert)
-
-```
-reconstructs_A:
-    Arrange: fixed 4×3 A
-    Act:     Decompose(A)
-    Assert:  Q·R ≈ A  (within 1e-5)
-
-Q_columns_are_orthonormal:
-    Assert:  Qᵀ·Q ≈ I  (Cols×Cols)
-
-R_is_upper_triangular:
-    Assert:  entries below the diagonal of R are ≈ 0
-
-solves_overdetermined_least_squares:
-    Arrange: line-fit system with more rows than unknowns
-    Assert:  solution matches the normal-equation reference (within tol)
-
-matches_known_3x3:
-    Arrange: golden 3×3 from Golub & Van Loan
-    Assert:  R (up to column sign) equals the documented factor
-
-apply_qtranspose_equals_explicit:
-    Assert:  ApplyQtranspose(b) == Q()ᵀ · b
-
-rank_deficient_returns_false:
-    Arrange: A with a dependent column
-    Assert:  Decompose returns false (near-zero R pivot)
-
-givens_update_adds_row:
-    Arrange: factor A, then stream one extra row via Givens
-    Assert:  updated R matches a full re-factorization of the augmented A
-
-identity_factors_to_identity:
-    Assert:  Q ≈ I, R ≈ I
-```
-
-## Reference vectors
-
-- Golden 3×3 with documented `Q`,`R` (sign convention fixed by `−sign(x₀)‖x‖`).
-- Overdetermined fit ⇒ least-squares solution equal to the `(AᵀA)⁻¹Aᵀb` reference.
-
-## Edge cases
-
-- Rank-deficient / dependent columns ⇒ zero pivot flagged.
-- Square vs tall shapes; `Cols == 1` (a single reflector).
-- Sign ambiguity of columns — assertions compare up to per-column sign.