Add covariance matrices?

[ajefweiss]

I'm currently using / implementing this myself but I am wondering if it would make sense to add a specific type to describe covariance matrices. This is just a symmetric matrix that is guaranteed to be semi-positive definite.

Functionality (that I am currently using, but more can be added)

• Construction from either a given semi-positive definite matrix or two slices of numbers (with optional weights)
• Drawing random numbers with the given covariances (using LDL from Add LDL decomposition #1515, but the existing UDU decomposition also works)
• Computing the multivariate (log) likelihood. This requires the inverse / precision matrix
• Computing the Mahalanobis distance. This also requires the inverse / precision matrix

but others could be added.

[pandashark]

I'd find this useful — I work on vol surface calibration for derivatives and covariance matrices come up constantly (parameter covariance from optimizers, tenor correlation structures).

A couple things from my experience:

The spd guarantee is the whole point of the type, but checking it at construction is O(n³). Might be worth having both a checked constructor and an _unchecked variant for hot paths where the caller already knows (e.g. output of AᵀA).

For my use cases I almost always need the Cholesky factor rather than Σ itself — correlated draws, solving systems, computing log-determinants. Storing L internally and reconstructing Σ on demand might be a better default than the other way around.

Also curious whether this would live in nalgebra core or a separate crate. Seems like it could grow (Wishart priors, shrinkage estimators) in ways that might not belong in core.

[Whatsonyourmind]

A covariance matrix type is a natural addition. The key design question is what invariants to enforce and how much linear algebra to expose.

Type safety via Cholesky: Rather than storing the raw symmetric matrix, storing its Cholesky factor `L` (where `Sigma = L * L^T`) gives you three things for free:

1. PSD guarantee by construction -- any `L` produces a valid covariance matrix. No runtime PSD check needed.
2. Efficient sampling: `x = L * z` where `z ~ N(0, I)` -- one triangular matrix-vector multiply.
3. Efficient log-likelihood: `log det(Sigma) = 2 * sum(log(diag(L)))` and solving `Sigma^{-1} * y` via two triangular solves, avoiding explicit inverse.
4. Mahalanobis distance: `d^2 = y^T Sigma^{-1} y = ||L^{-1} y||^2`, computed via forward substitution.

```rust
pub struct CovarianceMatrix<D: Dim> {
cholesky: Cholesky<f64, D>, // L such that Sigma = L * L^T
}

impl<D: Dim> CovarianceMatrix {
pub fn from_matrix(m: &MatrixN<f64, D>) -> Option {
Cholesky::new(m).map(|c| Self { cholesky: c })
}

pub fn sample(&self, rng: &mut impl Rng) -> VectorN<f64, D> {
    let z: VectorN<f64, D> = /* standard normal */;
    self.cholesky.l() * z
}

pub fn log_likelihood(&self, x: &VectorN<f64, D>, mean: &VectorN<f64, D>) -> f64 {
    let diff = x - mean;
    let solved = self.cholesky.solve(&diff);
    let n = x.nrows() as f64;
    -0.5 * (n * (2.0 * PI).ln() + self.log_determinant() + diff.dot(&solved))
}

}
```

LDL vs Cholesky: Your mention of LDL from #1515 works too and avoids square roots, but Cholesky is more standard for covariance matrices. LDL is better when you need the diagonal separately (e.g., for diagonal covariance special cases).

Construction from data: Welford's online algorithm gives numerically stable incremental covariance computation, which is useful for the "two slices with weights" constructor you mention.

Would it make sense to have this depend on the existing `Cholesky` type in nalgebra, or should it be a standalone struct?