Inverting non positive definite covariance matrix

I have an expression for a covariance matrix $C$ in terms of the indices $i$ and $j$. In this way I can analytically calculate the elements of my covariance matrix, however when I try to invert $C$ matlab gives a warning about the matrix being close to singular. The inversion therefore doesn't work, by which I mean the multiplying $C$ by the inverse given i do not get the identity. I have tried calculating the pseudo inverse but this also does not work. I have also tried adding a small constant along the diagonal but again the results do not work. In general I'm working with matrices of dimension 1200 but I will give low dimensional example matrix that has the same properties, i.e. the matrices are symmetric about the diagonal and the anti-diagonal

19.9939  19.9954  19.9958  19.9951  19.9933  19.9905  19.9865
19.9954  19.9973  19.9981  19.9978  19.9965  19.9940  19.9905
19.9958  19.9981  19.9993  19.9995  19.9985  19.9965  19.9933
19.9951  19.9978  19.9995  20.0000  19.9995  19.9978  19.9951
19.9933  19.9965  19.9985  19.9995  19.9993  19.9981  19.9958
19.9905  19.9940  19.9965  19.9978  19.9981  19.9973  19.9954
19.9865  19.9905  19.9933  19.9951  19.9958  19.9954  19.9939

As mentioned in the title the matrix isn't positive, however the the negative eigenvalues are very small suggesting that the matrix is not positive only due to machine precision. The negative eigenvalues are $-1.4048e-14$ and $-2.4571e-15$. How can I go about inverting these matrices?

• In many algorithms it is not necessary actually to invert a matrix. Therefore, please tell us why you wish to invert these matrices. It would also be useful to know how such incredibly ill-conditioned matrices arise, because it is possible your original problem could be solved in a different way altogether or it would at least suggest alternative algorithms. – whuber Mar 21 '14 at 15:10
• I have a set of $N$ measurements that are distributed in a multivariate normal distribution and I'm trying to calculate the covariance matrix of these measurements. I believe the reason that the matrix is so poorly conditioned is because the correlations between measurements are very strong. The reason I need to invert the matrix is to calculate the Fisher information and ultimately to perform a scoring algorithm for parameter estimation. – mrkprc1 Mar 21 '14 at 15:22

Assuming your matrices are like the illustrated one, which I characterize has having all its coefficients close to a constant array (with the constant $20$), you can employ the Sherman-Morrison Formula. This formula relates the inverse of a matrix to the inverse of a perturbation of that matrix. It has been used in statistics for iterative and online formulas for performing least squares regressions, because the inclusion of each new row of data to the design matrix $\mathbb{X}$ has the effect of adding a rank-one matrix to the "sum of squares and products matrix" $\mathbb{X}^\prime \mathbb{X}$.

To illustrate, the matrix $\mathbb{A}$ in the question can be written in the form

$$\mathbb{A} = \mathbb{Y} + 20 (1\cdot 1^\prime)$$

where the entries in $\mathbb{Y}$ appear to lie between $0.00$ and $-0.01$ and $1$ is the $7$-vector (a column) whose components are all ones. The rank-one matrix $20(1\cdot 1^\prime)$ (all of whose components equal $20$) is the perturbation of $\mathbb{Y}$ that yields $\mathbb{A}$. The formula asserts

$$\mathbb{A}^{-1} = \mathbb{Y}^{-1} - \frac{20}{1 + 20 (1^\prime \mathbb{Y}^{-1} 1)} (\mathbb{Y}^{-1} 1)\cdot(1^\prime \mathbb{Y}^{-1}).$$

This amounts to inverting $\mathbb{Y}$ and adjusting its coefficients in terms of the row (or column) sums of that inverse. Provided $\mathbb{Y}$ can be stably inverted, this might be an attractive approach.

As I mentioned in a comment, it is probably unnecessary to invert $\mathbb{A}$ at all. For Fisher scoring, for instance, where the neat mathematical formulas have expressions of the form "$x = \mathbb{A}^{-1}b$" (where typically $b$ is the value of the score function), you should be solving the system

$$\mathbb{A}x = b$$

instead of pre-inverting $\mathbb{A}$ and multiplying its inverse by $b$. This tends to lead to more stable and accurate results.

• Is there a book where I can find this sort of thing? – mrkprc1 Mar 24 '14 at 11:42
• There is (but I was unable to find it when writing my answer): see The Matrix Cookbook, section 3.2.4 (p. 18). – whuber Mar 24 '14 at 14:31

Your case is an example where not the diagonal elements, but the distances hiding behind the coefficients are to blame for the not-positive-definiteness of the matrix. I described it here as "Cause 1" (displayed on Fig. 1). For short, your 7 items as points in space do not converge in Euclidean space. That cause is unlikely to be cured by adding a constant to the diagonal of the covariance matrix. You said it did not help you. Without checking your words, I tried to add a constant to all the elements of the matrix; it didn't help either.

Below are the outputs of what I did with your data then. First shown are the eigenvalues of your matrix as it is. I converted the matrix into Euclidean distance matrix by means of the Law of cosines.

I then performed double centering of that distance matrix. Double centering finds geometric centroid of the cloud of points between which the distances are, and converts the distances into new covariate-like scalar products. The resultant matrix describes the distances, but their entries are angle-type similarities, so we are in right to eigen-decompose it the usual way. I don't show the doubly-centered matrix itself, but I show its eigenvalues. Note that there are negative ones, their presence is the evidence of "Cause 1" of not-positive-definiteness.

So, the distances, that are behind the covarinces, should be inflated with a small constant. I added the constant, and the thus corrected Euclidean distances, when doubly-centered, showed no negative eigenvalues. Good.

Now, what we need is to transform distanced back to covariances. We do it again by the Law of cosines, using the diagonal of the input covariance matrix. The resultant covariance matrix has this diagonal and has corrected off-diagonal values, and is now positive definite.

Eigenvalues of the input cov matrix:
139.9695727
.0114858
.0000673
.0000601
-.0000215
-.0000530
-.0001113

Euclidean distances:
.00000  .02000  .04000  .06083  .08124  .10100  .12166
.02000  .00000  .02000  .04123  .06000  .08124  .10100
.04000  .02000  .00000  .01732  .04000  .06000  .08124
.06083  .04123  .01732  .00000  .01732  .04123  .06083
.08124  .06000  .04000  .01732  .00000  .02000  .04000
.10100  .08124  .06000  .04123  .02000  .00000  .02000
.12166  .10100  .08124  .06083  .04000  .02000  .00000

Eigenvalues of the double centered matrix of distances:
10 ** -2   X
1.148575892
.006728980
.006090536
.000000000
-.002108362
-.005304872
-.011125031

Corrected Euclidean distances:
.00000  .02965  .04965  .07048  .09089  .11065  .13131
.02965  .00000  .02965  .05088  .06965  .09089  .11065
.04965  .02965  .00000  .02697  .04965  .06965  .09089
.07048  .05088  .02697  .00000  .02697  .05088  .07048
.09089  .06965  .04965  .02697  .00000  .02965  .04965
.11065  .09089  .06965  .05088  .02965  .00000  .02965
.13131  .11065  .09089  .07048  .04965  .02965  .00000

Eigenvalues of the double centered corrected matrix of distances:
.013491933
.000636840
.000269820
.000222552
.000191635
.000000993
.000000000

Corrected cov matrix:
19.99390    19.99516    19.99537    19.99447    19.99247    19.98948    19.98528
19.99516    19.99730    19.99786    19.99736    19.99587    19.99317    19.98948
19.99537    19.99786    19.99930    19.99929    19.99807    19.99587    19.99247
19.99447    19.99736    19.99929    20.00000    19.99929    19.99736    19.99447
19.99247    19.99587    19.99807    19.99929    19.99930    19.99786    19.99537
19.98948    19.99317    19.99587    19.99736    19.99786    19.99730    19.99516
19.98528    19.98948    19.99247    19.99447    19.99537    19.99516    19.99390

Its eigenvalues:
139.9661876
.0134919
.0006356
.0002698
.0002226
.0001915
.0000010
• So after the double centering could you elaborate on what you mean by "should be inflated by a small distance". Also, how do I get back to the Euclidean distance matrix when I have done this? – mrkprc1 Mar 24 '14 at 15:07
• It simply means: you should add some small positive constant to the distance matrix. Add, then do double centering to check if its eigenvalues has become all nonnegative. If not, add again, etc. repeat. (Please remember that double centering is done on squared distances. But you add a constant to no-squared distances. Also note that eigenvalues of doubly centered distance matrix always has one zero eigenvalue; you check that there is no negative ones.) – ttnphns Mar 24 '14 at 16:11
• Ok great, again is there a book that I might find that covers this technique? – mrkprc1 Mar 24 '14 at 17:38
• No specific book comes to my mind this time, sorry. – ttnphns Mar 25 '14 at 5:57
• Ok, well the method seems to work, in that the output matrix looks similar to the covariance matrix given and is invertable, however with no proof of why it works and without a book where I can verify the method I don't know if I'm comfortable to use this method. – mrkprc1 Mar 25 '14 at 12:14