# Numerical Instability of calculating inverse covariance matrix

I have a 65 samples of 21-dimensional data (pasted here) and I am constructing the covariance matrix from it. When computed in C++ I get the covariance matrix pasted here. And when computed in matlab from the data (as shown below) I get the covariance matrix pasted here

Matlab code for computing cov from data:

data = csvread('path/to/data');
matlab_cov = cov(data);


As you can see the difference in covariance matrices are minute (~e-07), which is probably due to numerical problems in the compiler using floating point arithmetic.

However, when I compute the pseudo-inverse covariance matrix from the covariance matrix produced by matlab and the one produced by my C++ code, I get widely different results. I am computing them in the same way i.e.:

data = csvread('path/to/data');
matlab_cov = cov(data);
matlab_inv = pinv(matlab_cov);
my_inv = pinv(my_cov);


The difference is so huge that when I am computing the mahalanobis distance from a sample (pasted here) to the distribution of the 65 samples by:

$(65/64^2) \times ((sample-mean)\times {\sum}^{-1} \times (sample-mean)')$

using the different inverse covariance matrices (${\sum}^{-1}$) I get widely different results i.e.:

 (65/(64^2))*((sample-sample_mean)*my_inv*(sample-sample_mean)')
ans =

1.0167e+05

(65/(64^2))*((sample-sample_mean)*matlab_inv*(sample-sample_mean)')
ans =

109.9612


Is it normal for the small (e-7) differences in covariance matrix to have such an effect on the computation of the pseudo-inverse matrix? And if so, what can I do to mitigate this effect?

Failing this, are there any other distance metrics I can use that do not involve the inverse covariance? I use the Mahalanobis distance as we know for n samples it follows a beta distribution, which I use for hypothesis testing

EDIT: Adding C++ code for calculating covariance matrix below: The vector<vector<double> > represents the collection of rows from the file pasted.

Mat covariance_matrix = Mat(21, 21, CV_32FC1, cv::Scalar(0));
for(int j = 0; j < 21; j++){
for(int k = 0; k < 21; k++){
for(std::vector<vector<double> >::iterator it = data.begin(); it!= data.end(); it++){
covariance_matrix.at<float>(j,k) += (it->at(j) - mean.at(j)) * (it->at(k) - mean[k]);
}
covariance_matrix.at<float>(j,k) /= 64;
}
}

• Inverting matrixes..... That's a dangerous thing! Ususally it's preferable to find alternatives to that (e.g. pseudoinverse) Feb 27, 2013 at 14:54
• @Aly: the matrices you are looking to invert are not "valid" covariances matrices because they are not positive definite; numerically they even have some eigenvalues that are negative (but close to zero). I would probably just add some very small constant along the diagonal; it is a form of Tikhonov correction really ($Χ + \lambda I$). Also don't use floats, use doubles to store your covariance matrix. (And besides you already use OpenCV, you might as well use Eigen or Armadillo..) Feb 27, 2013 at 15:28
• @Aly: Wikipedia, really. (it is the lemma : Tikhonov Regularization). The method that whuber mentioned using the SVD will give you a non-negative definite matrix if you set small eigenvalues to zero; you will still need to add some small constant to all your eigenvalues to make them positive definite. Practically both methods do the same. Just I resorted in not using the SVD but directly affect the samples eigenvalues by adding $\lambda$ to all of them. I haven't come across any references, both methods are quite intuitive I believe. Feb 27, 2013 at 15:57
• @user11852 Please can you make your comments an answer, I am still experimenting, but if promising will accept. Also, if others make their suggestions answers I will up vote as they have been very helpful/useful to my understanding of the problem
– Aly
Feb 27, 2013 at 16:22
• I commented in your other thread that having variables that sum to 1, as you data set does, encourages instability and contains a redundant variable. Please try dropping one column. You do not even need the pinv: the covariance matrix is no longer singular. Feb 27, 2013 at 17:58

The matrices you are looking to invert are not "valid" covariances matrices because they are not positive definite; numerically they even have some eigenvalues that are negative (but close to zero). This most probably due to machine zeros, for example the last eigenvalue of your "matlab_covariance" matrix is -0.000000016313723. To correct to positive definite you can do two things :

1. Just add some very small constant along the diagonal; a form of Tikhonov correction really ($Χ+\lambda I$ with $\lambda \rightarrow 0$).
2. (Based on what whuber proposed) Use SVD, set the "problematic" eigenvalues to a fixed small value (not zero), reconstruct your covariance matrix and then invert that. Clearly if you set some of those eigenvalues to zero you will end up with a non-negative (or semi-positive) matrix, that will not be invertible still.

A non-negative matrix does not have a inverse but it does have a pseudo inverse (all matrices with real or complex entries have a pseudo-inverse), nevertheless the Moore–Penrose pseudo-inverse is more computationally expensive than a true inverse and if the inverse exists it is equal to the pseudo-inverse. So just go for the inverse :)

Both methods practically try to handle the eigenvalues that evaluate to zero (or below zero). The first method is bit hand-wavy but probably much faster to implement. For something slightly more stable you might want to compute the SVD and then set the $\lambda$ equal to absolute of the smallest eigenvalue (so you get non-negative) plus something very small (so you get positive). Just be careful not to enforce positivity to a matrix that is obviously negative (or already positive). Both methods will alter the conditioning number of your matrix.

In statistical terms what you do by adding that $\lambda$ across the diagonal of your covariance matrix you add noise to your measurements. (Because the diagonal of the covariance matrix is the variance of each point and by adding something to those values you just say "the variance at the points I have readings for is actually a bit bigger than what I thought originally".)

A fast test for the positive-definiteness of a matrix is the existence (or not) of the Cholesky decomposition of it.

Also as a computational note:

1. Don't use floats, use doubles to store your covariance matrix.
2. Use numerical linear algebra libraries in C++ (like Eigen or Armadillo) to get inverses of matrices, matrix products, etc. It's faster, safer and more concise.

EDIT: Given you have a Cholesky decomposition of your matrix $K$ such that $LL^T$ (you have to do that to check you are having a Pos.Def. matrix) you should be able to immediately solve the system $Kx=b$. You just solve Ly = b for y by forward substitution, and then L^Tx = y for x by back substitution. (In eigen just use the .solve(x) method of your Cholesky object) Thanks to bnaul and Zen for pointing out that I focused so much on getting the $K$ be Pos.Def. that I forgot why we cared about that in the first place :)

• +1. Using Mathematica and applying it to the data (rather than the posted covariance matrix, which may have been presented with too little precision) I find no negative eigenvalues. That is as it should be: when a covariance matrix is computed exactly, it's guaranteed positive semi-definite, so any negative eigenvalues must be attributed to imprecision in the calculations. Any decent generalized inverse procedure ought to "recognize" those tiny negative values as zeros and treat them accordingly.
– whuber
Feb 27, 2013 at 17:52
• Thanks guys for the effort, as stated I have up voted and will try these out and either comment or accept accordingly.
– Aly
Feb 27, 2013 at 18:11
• Sorry, I'm a bit confused, how does solving the Cholesky give use the Mahalanobis distance?
– Aly
Feb 28, 2013 at 14:13
• Check the link in original bnaul's post. But don't use $LU$ but Cholesky (that's what they mean by LDL*). Feb 28, 2013 at 14:28
• Is there a good way to estimate how large an $\lambda$ I need? I used$\lambda=1e-6$, and it was fine until suddenly it wasn't and my program crashed unable to make a Cholesky factorization. May 6, 2022 at 3:09

The posted answers and comments all make good points about the dangers of inverting nearly singular matrices. However, as far as I can tell, no one has mentioned that computing Mahalanobis distance doesn't actually require inverting the sample covariance. See this StackOverflow question for a description of how to do so using the $LU$ decomposition.

The principle is the same as solving a linear system: when trying to solve for $x$ such that $Ax=b$, there are much more efficient and numerically stable methods than taking $x=A^{-1}b$.

Edit: it probably goes without saying, but this method produces the exact distance value, whereas adding $\lambda I$ to $S$ and inverting yields only an approximation.

• You're right, @bnaul. However, without some kind of regularization, the LU decomposition will not work either. I will add a comment about this in my answer.
– Zen
Feb 28, 2013 at 4:24
• @bnaul: why do the LU when you do to the Cholesky to check positive definiteness? Assuming you have a valid covariance matrix $K=LL^T$ solving $Ly = b$ for y by forward substitution, and then $L^Tx = y$ for x by back substitution will be faster. Good point though, I definitely focus on get positive-definiteness that I forgot why I cared about it originally! :D Feb 28, 2013 at 8:46

(Years later) a tiny example: with $A$ rank-deficient, $r < n, \ n - r$ eigenvalues of $A^T A$ will be 0 to within machine precision -- and about half of these "zeros" may be $< 0$ :

#!/usr/bin/env python2
""" many eigenvalues of A'A are tiny but < 0 """
# e.g. A 1 x 10: [-1.4e-15 -6.3e-17 -4e-17 -2.7e-19 -8.8e-21  1e-18 1.5e-17 5.3e-17 1.4e-15  7.7]

from __future__ import division
import numpy as np
from numpy.linalg import eigvalsh  # -> lapack_lite
# from scipy.linalg import eigvalsh  # ~ same
from scipy import __version__

np.set_printoptions( threshold=20, edgeitems=10, linewidth=140,
formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
print "versions: numpy %s  scipy %s \n" % (
np.__version__, __version__  )

np.random.seed( 3 )

rank = 1
n = 10
A = np.random.normal( size=(rank, n) )
print "A: \n", A
AA = A.T.dot(A)
evals = eigvalsh( AA )
print "eigenvalues of A'A:", evals