# Singular covariance matrix in Mahalanobis distance in Matlab

I am using the Mahalanobis distance to classify an unknown 64-dimensional vector into one of 75 classes. There are n samples of 64-dimensional vectors for each class, arranged into an Nx64 matrix format. Reading this into Matlab with no problem.

However when the 64x64 covariance matrix is calculated using the cov() function it is not being positive-definite or non-singular. Hence the inverse of the covariance matrix required to calculate the Mahal. dist. is undefined.

Now the data that is used to get the feature vector is basically co-ordinates of points taken from handwriting samples. In such a case it seems that the distribution is naturally giving rise to a singular covariance matrix. Is there any way to apply the Mahalanobis distance, or some way of modifying the vectors to yeild a non-singular covar. matrix?

• Can you tell us the condition number of your covariance matrix using cond(A)? – Zen Sep 21 '12 at 13:37
• 6.8396e+052 for the covariance matrix. the data is from file uploadmb.com/dw.php?id=1348237753 – AruniRC Sep 21 '12 at 14:32
• Your matrix is horribly ill conditioned. Read Michael's comment and try to reduce the dimensionality of your problem. You may also try some "hack" to regularize your matrix, like adding a very small number to each element of the diagonal. But I encourage you to think about why the matrix is so ill conditioned. Some info: en.wikipedia.org/wiki/Condition_number – Zen Sep 21 '12 at 19:34

To add on BGreene 's answer, you can use the Moore-Penrose inverse. When you use the Mahalanobis distance modified in this way, on the same data used to estimate the covariance matrix, you lose nothing. The covariance matrix is singular because your data happen to live in a linear subspace, and your modified Mahalanobis distance is identical to the Mahalanobis distance you could compute by first transforming your feature set to get a lesser set generating the same space. If you later use it on new data, which happens to have components outside that linear space, that components will just be zeroed in the distance. But if your data are reasonably representative of what you will see later, that should not be a problem.

• I am glad to see that some people probably read my comments and agree. But why does the Moore-Penrose inverse do the right thing? +1 @kjetil for making some very astute observations! – Michael R. Chernick Sep 22 '12 at 15:42
• thanks. echoing what @MichaelChernick said, could you provide some links/theory about how and why the modified Mahal. dist. using pseudo-inverse would be same as reducing dimensionality of the feature set. – AruniRC Sep 22 '12 at 17:28
• Can you look at my new answer and see if that helps you to work it out? The facts I am referencing must be very well known, but I have never seen it being written out (except in my own notes), maybe somebody else have a reference. – kjetil b halvorsen Sep 22 '12 at 18:41

You could try using the Moore-Penrose pseudoinverse (pinv function in Matlab), to invert your covariance matrix

• ah yes! thanks. the question is whether the pseudo-inverse can be used successfully in the Mahalanobis distance calculation instead of the usual inverted covariance matrix. trying it out asap. if it works, will accept this. – AruniRC Sep 21 '12 at 14:43
• Before you decide that the covariance matrix "should be singular" and that you can just go ahead and use a pseudoinverse think about what you are doing. You don't have the covariance matrix. What you have is a sample estimate. How big is N. 64 dimensions for the vector is pretty large. It is possible that N is too small for the sample covariance matrix to be nonsingular even if the population covariance is nonsingular. – Michael R. Chernick Sep 21 '12 at 16:44
• Also consider what it means if the population covariance matrix is singular. Then the data sit in a proper subspace of the 64 dimensional Euclidean space. So there is no 64-dimensional distance that can be defined! If all the classes have the same singular covariance matrix then you can use a pooled estimate of the covariance and you can think of the problem in a lower dimensional space where everything lives and the lower dimensional covariance matric is nonsingular. – Michael R. Chernick Sep 21 '12 at 16:51
• So perhaps a Mahalanobis distance can be defined there and maybe there is a pseudoinverse to the population covariance matrix that would work. That is some theory but what about the practical problem where you use an estimate. I think the sample covariance matrix will converge to the population covariance matrix. But for a given sample size N will the sample matrix be singular or nonsingular? If it is singular will it be guaranteed to lie in the same subspace as the population covariance matrix. I think it is possible to be nonsingular for a given data set and a given sample size N. – Michael R. Chernick Sep 21 '12 at 16:55
• I think it is also possible that it could be singular for some N. Eventually it would I think because it is converging to a singular matrix. But I think that at a given N when it is singular it could be that the sample covariance is in a different subspace. It is not an easy problem and I do not see a straightforward solution. – Michael R. Chernick Sep 21 '12 at 16:58

The typical handwriting dataset everybody seems to use is quite ill formed for PCA and many mathematical methods. It is pixels, and obviously some of these (bottom right corner, whatever) are never really painted in any of the samples you have.

Even worse, the values are discrete. This can cause all kind of artifacts.

PCA now comes from the perfect world where all dimensions are continuous and just have a different amount of variance to them. And it assumes the variance is caused by importance, not by natural scale of the axes ...

Now throw in that we are (probably) talking about pixel data here. For obvious reasons, we can expect neighboring pixels to be strongly correlated. Guess what PCA will do... you could as well just downsample it to a lower resolution image anyway.

Consider using something else. Just working around the matrix inversion will probably not save your analysis. And maybe you can do without the dimensionality reduction, too?

• fair point. making a better set of features is an option that seems viable to consider. would that make a difference? – AruniRC Sep 22 '12 at 1:14
• Look into methods that do not require dimensionality reduction or PCA. Given that each pixel represents essentially a comparable thing, I'm not convinced that Mahalanobis will provide a substantial benefit over plain Euclidean, actually. – Has QUIT--Anony-Mousse Sep 22 '12 at 1:38

This is really a comment but to long: Thanks Michael! The way to think about the Moore-Penrose inverse is as follows: Any Matrix represents a linear operator. First some notation: Let $A$ be an $n \times m$-matrix, where we assume for simplicity that $n \ge m$. This matrix really represents a linear transformation: $A \colon {\mathbb R^m } \mapsto {\mathbb R^n}$. Let ${\mathcal N}(A)$ be the nullspace of $A$. Then we can decompose $A$ as a direct sum of two operators, on acting on the nullspace (sending it to the zero vector in $\mathbb R^n$), the other acting on the ortogonal complement of the nullspace. Now the singular value decomposition is giving us an coordinate sytem adapted to this situation. Suppose the rank of $A$ is $r \le \min(n,m)$. Then we can write the SVD as $$A = U \Lambda V^T = [U_1\colon U_0] \begin{pmatrix} \Lambda_1 & 0 \\ 0 & \Lambda_0\end{pmatrix} = U_1 \Lambda_1 V_1^T$$ Where $\Lambda_0=0$, $\Lambda_1$ is $r \times r$ and $U_1$ has $r$ columns, $V_1$ is $m\times r$. Note that the the coumns of $V$ gives an orthogonal basis for $\mathbb R^m$, the first $r$ of which, that is, the columns of $V_1$ gives a basis for ${\mathcal N}(A)^{\perp}$. Now, writing a general point in this basis, using $x$ for the coordinates, we get $$x_1 v_1 + \dots + x_r v_r + x_{r+1} v_{r+1} \dots x_m v_m$$ which we can write as $V_1 x_1 + V_0 x_0$ where now $X_1, x_0$ are subvectors of the vector $x$. Now letting $A$ act we find that $$A (V_1 x_1 + V_0 x_0) = U_1 \Lambda_1 V_1^T V_1 x_1 + 0 = U_1 \begin{pmatrix} \lambda_1 x_1 \\ \vdots \\ \lambda_r x_r \end{pmatrix}$$ Observe that the $r$ columns of $U_1$ forms an orthogonal basis of the image space of $A$ in $\mathbb R^n$.

Call now $A$ reduced to acting on ${\mathcal N}(A)$ for the nonsingulatr part of $A$. A then consists of the direct sum of its nonsingular part and the zero operator acting on the nullspace. We get btyhne Moore-penrose inverse by taking a direct sum of the usual inverse of the nonsingular part and a zero operator.

All other generalized inverses can be otained in this way, as a direct sum of the usual inverse of the nonsingular part, and some arbitrary operator $B$ replacing the zero operator in the case of Moore-Penrose. This explains the special role of the Moore-Penrose generalized inverse.