embedded-pro · gabrielfrasantos · Aug 6, 2026 · Aug 6, 2026 · Aug 6, 2026 · Aug 6, 2026
diff --git a/doc/estimators/README.md b/doc/estimators/README.md
@@ -17,9 +17,3 @@ Statistical estimation algorithms for fitting models to observed data and making
 | Algorithm                                           | Description                                                              |
 |-----------------------------------------------------|--------------------------------------------------------------------------|
 | [Recursive Least Squares](RecursiveLeastSquares.md) | Sample-by-sample parameter estimation with exponential forgetting factor |
-
-## Consistency Metrics
-
-| Algorithm                                               | Description                                                                         |
-|---------------------------------------------------------|-------------------------------------------------------------------------------------|
-| [Consistency Metrics (NEES/NIS)](ConsistencyMetrics.md) | Normalised Estimation Error Squared and Normalised Innovation Squared with χ² gates |
diff --git a/doc/filters/active/UnscentedKalmanFilter.md b/doc/filters/active/UnscentedKalmanFilter.md
@@ -146,7 +146,7 @@ graph LR
 |------------------------------------------------------------------|---------------------------------------------------------------------------------------|
 | [Kalman Filter](KalmanFilter.md)                                 | The UKF reduces to the standard KF when $f$ and $h$ are linear                        |
 | [Extended Kalman Filter](ExtendedKalmanFilter.md)                | Uses Jacobians instead of sigma points; simpler but less accurate for nonlinear cases |
-| [Cholesky Decomposition](../../solvers/CholeskyDecomposition.md) | Used internally to generate sigma points from the covariance matrix                   |
+| [Cholesky Decomposition](../../math/CholeskyDecomposition.md) | Used internally to generate sigma points from the covariance matrix                   |
 
 ## References & Further Reading
 

diff --git a/doc/solvers/CholeskyDecomposition.md → doc/math/CholeskyDecomposition.md b/doc/solvers/CholeskyDecomposition.md → doc/math/CholeskyDecomposition.md
@@ -64,8 +64,9 @@ $$L = \begin{bmatrix} 2 & 0 \\ 1 & 2 \end{bmatrix}$$
 
 | Algorithm                                                             | Relationship                                                      |
 |-----------------------------------------------------------------------|-------------------------------------------------------------------|
-| [Unscented Kalman Filter](../filters/active/UnscentedKalmanFilter.md) | Uses Cholesky to generate sigma points from the covariance matrix |
-| [Gaussian Elimination](GaussianElimination.md)                        | General-purpose alternative; does not exploit symmetry            |
+| [Unscented Kalman Filter](../filters/active/UnscentedKalmanFilter.md) | Uses `Factor` to generate sigma points from the covariance matrix |
+| [Consistency Metrics](ConsistencyMetrics.md)                          | Uses `Solve` for the SPD covariance solve behind NEES/NIS         |
+| [Gaussian Elimination](../solvers/GaussianElimination.md)             | General-purpose alternative; does not exploit symmetry            |
 
 ## References & Further Reading
 

diff --git a/doc/estimators/ConsistencyMetrics.md → doc/math/ConsistencyMetrics.md b/doc/estimators/ConsistencyMetrics.md → doc/math/ConsistencyMetrics.md
@@ -86,7 +86,7 @@ The average NEES over a Monte-Carlo ensemble of $M$ runs and $K$ time steps yiel
 
 ## Connections to Other Algorithms
 
-NEES and NIS are statistical companions to the standard Kalman filter update step. They depend on the covariance propagation produced by the `filters/active` family (KF, EKF, UKF). The linear solve reuses `solvers::GaussianElimination`, and the state-error representation aligns with `estimators::EstimationMetrics`.
+NEES and NIS are statistical companions to the standard Kalman filter update step. They depend on the covariance propagation produced by the `filters/active` family (KF, EKF, UKF). The covariance being symmetric positive-definite, the linear solve reuses `math::CholeskyDecomposition::Solve` (factorisation plus forward/back substitution), keeping the metric a dependency-free `math` primitive.
 
 ## References & Further Reading
 

diff --git a/doc/math/README.md b/doc/math/README.md
@@ -11,7 +11,9 @@ Core mathematical primitives for numerical computation.
 | [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                       |
+| [Triangular Solve](TriangularSolve.md)              | Lower/upper-triangular forward and back-substitution shared by Gaussian elimination, Cholesky, and QR |
+| [Cholesky Decomposition](CholeskyDecomposition.md)  | SPD factorization $A=LL^T$ with `Factor` (L) and `Solve` (SPD system) — twice as fast as LU     |
 | [Matrix Operations](MatrixOperations.md)            | Structural matrix utilities — `Symmetrize` (closest symmetric matrix)                          |
 | [Step Response Metrics](StepResponseMetrics.md)     | Rise time, settling time, percent overshoot, peak time, and steady-state error from a bounded step-response vector |
 | [Matrix Exponential](MatrixExponential.md)          | Scaling-and-squaring with diagonal (6,6) Padé approximant — exact ODE solution operator and discretisation engine |
+| [Consistency Metrics](ConsistencyMetrics.md)        | NEES/NIS estimator consistency with χ² gates — normalised estimation error squared and normalised innovation squared |
diff --git a/doc/solvers/JacobiEigenSolver.md b/doc/solvers/JacobiEigenSolver.md
@@ -119,11 +119,11 @@ graph LR
     JAC -.->|"eigenvalue floor keeps SPD"| CHOL
 ```
 
-| Algorithm                                          | Relationship                                                                                                          |
-|----------------------------------------------------|-----------------------------------------------------------------------------------------------------------------------|
-| [QR Decomposition](QrDecomposition.md)             | Both are built from orthogonal (Givens/Householder) transforms; QR underlies the alternative tridiagonal eigen-method |
-| [Cholesky Decomposition](CholeskyDecomposition.md) | Requires symmetric positive-definite input; Jacobi eigenvalues certify or restore definiteness                        |
-| [Spectral Radius](SpectralRadius.md)               | Returns only the dominant eigenvalue magnitude; Jacobi returns the full spectrum and vectors                          |
+| Algorithm                                                  | Relationship                                                                                                          |
+|------------------------------------------------------------|-----------------------------------------------------------------------------------------------------------------------|
+| [QR Decomposition](QrDecomposition.md)                     | Both are built from orthogonal (Givens/Householder) transforms; QR underlies the alternative tridiagonal eigen-method |
+| [Cholesky Decomposition](../math/CholeskyDecomposition.md) | Requires symmetric positive-definite input; Jacobi eigenvalues certify or restore definiteness                        |
+| [Spectral Radius](SpectralRadius.md)                       | Returns only the dominant eigenvalue magnitude; Jacobi returns the full spectrum and vectors                          |
 
 ## References & Further Reading
 

diff --git a/doc/solvers/README.md b/doc/solvers/README.md
@@ -6,7 +6,6 @@ Numerical solvers for linear systems, polynomial roots, and matrix equations.
 
 | Algorithm                                                                  | Description                                                            |
 |----------------------------------------------------------------------------|------------------------------------------------------------------------|
-| [Cholesky Decomposition](CholeskyDecomposition.md)                         | Fast factorization for symmetric positive-definite matrices            |
 | [Gaussian Elimination](GaussianElimination.md)                             | Direct solver for dense linear systems using partial pivoting          |
 | [Levinson-Durbin](LevinsonDurbin.md)                                       | Fast solver for Toeplitz linear systems exploiting structural symmetry |
 | [Durand-Kerner](DurandKerner.md)                                           | Simultaneous iterative root-finder for polynomials                     |

diff --git a/numerical/estimators/CMakeLists.txt b/numerical/estimators/CMakeLists.txt
@@ -12,14 +12,8 @@ target_link_libraries(numerical.estimators ${NUMERICAL_VISIBILITY}
 )
 
 target_sources(numerical.estimators PRIVATE
-    ConsistencyMetrics.hpp
     Estimator.hpp
 )
 
-numerical_add_coverage_sources(numerical.estimators
-    ConsistencyMetrics.cpp
-)
-
 add_subdirectory(offline)
 add_subdirectory(online)
-add_subdirectory(test)
diff --git a/numerical/estimators/offline/LinearRegression.hpp b/numerical/estimators/offline/LinearRegression.hpp
@@ -4,8 +4,6 @@
 #pragma GCC optimize("O3", "fast-math")
 #endif
 
-#include "numerical/math/CompilerOptimizations.hpp"
-
 #include "numerical/estimators/Estimator.hpp"
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/solvers/QrDecomposition.hpp"

diff --git a/numerical/estimators/offline/YuleWalker.hpp b/numerical/estimators/offline/YuleWalker.hpp
@@ -6,6 +6,7 @@
 
 #include "numerical/math/CompilerOptimizations.hpp"
 #include "numerical/math/Matrix.hpp"
+#include "numerical/math/Statistics.hpp"
 #include "numerical/math/Toeplitz.hpp"
 #include "numerical/solvers/GaussianElimination.hpp"
 
@@ -44,10 +45,7 @@ namespace estimators
     template<typename T, std::size_t Samples, std::size_t Order>
     OPTIMIZE_FOR_SPEED void YuleWalker<T, Samples, Order>::Fit(const InputVector& x)
     {
-        mean = T(0.0f);
-        for (std::size_t i = 0; i < Samples; ++i)
-            mean += x[i];
-        mean = mean / T(static_cast<float>(Samples));
+        mean = math::Mean(x);
 
         auto r = ComputeAutocovariance(x, mean);