fn invert_psd_matrix(m: &Array2<f64>) -> Array2<f64>
Cholesky-based inversion of a symmetric PSD matrix.
Computes M^{-1} via L L^T = M then inverting via forward + backward substitution on the identity.
M^{-1}