Singular covariance matrix in Mahalanobis distance in Matlab

Question

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?

Answer 1

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?

Answer 2

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.

Answer 3

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

Answer 4

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.