5
$\begingroup$

I have a symmetric positive semidefinite covariance matrix $A$, which is approximately computed as the output of a quadratic regression. I then need to invert $A$, but often it is close to singular. I've reduced the problem by using scaling. That is, I create a diagonal matrix $D$, with elements $D_{ii} = 1/\sqrt{A_{ii}}$. Then

$A^{-1} = D(DAD)^{-1}D$

Where $DAD$ has a lower conditioning number then $A$. Unfortunately in some iterations this is not enough. The size of $A$ is quite small, say maximum $50 \times 50$.

I need the inverse of $A$ because I have to use it in a long calculation, where terms such as $x^TA^{-1}x$ and $A^{-1}B$ appear lots for time. Also: $A^{-1}$ represents a covariance matrix, so it has to be symmetric and positive definite.

Is there some better way to make $A$ invertible?

|cite|improve this question
$\endgroup$
  • $\begingroup$ Computing the pseudo-inverse instead is not an option? $\endgroup$ – Gabriel Landi Nov 18 '12 at 21:29
  • $\begingroup$ I didn't think about that. If I use the pseudo-inverse do I have some measure of the distance between the pseudo-inverse and the "true" inverse? $\endgroup$ – Jugurtha Nov 18 '12 at 23:14
  • $\begingroup$ When the matrix is invertible, the pseudo-inverse is the inverse. $\endgroup$ – Wolfgang Bangerth Nov 19 '12 at 0:25
5
$\begingroup$

Multiply the diagonal elements by a factor $q>1$ but close to 1. It will usually do the job. (I'd first try $q=1.0001$, but one can experiment with the number of zeros in this expression; e.g., use $q=1.01$ for very noisy data.) This is justified under certain conditions, as it is a specific form of regularization. For more on regularization, see my tutorial
http://mat.univie.ac.at/~neum/ms/regtutorial.pdf

An important exception is when some diagonal element is tiny or negative, in which case your data were insufficent for the attempted estimation.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thanks, could I do this while also using preconditioning? $\endgroup$ – Jugurtha Nov 19 '12 at 12:30
  • $\begingroup$ @Jugurtha: I didn't understand your two programming lines, so I can't see what you mean. Please rephrase these lines in ordinary mathematical terms. $\endgroup$ – Arnold Neumaier Nov 19 '12 at 15:19
  • $\begingroup$ I edited the question, I hope it is clearer now. $\endgroup$ – Jugurtha Nov 19 '12 at 16:15
  • $\begingroup$ @Jugurtha: Yes. But one would call what you do ''scaling'', not preconditioning. You can do the scaling after the diagonal regularization; but it should not need to be necessary, as inverting positive definite matrices is stable even without scaling. $\endgroup$ – Arnold Neumaier Nov 19 '12 at 18:53
4
$\begingroup$

The question to ask is why you need to invert the matrix. If a matrix is near-singular, it's true that you can define something like a pseudo-inverse in some stable way but it's nevertheless true that for near-singular matrices, solving linear systems in any way is unlikely to result in anything useful because the result is so strongly dependent on small changes in the right hand side.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ I added some motivation about why I want to invert $A$. $\endgroup$ – Jugurtha Nov 19 '12 at 12:31
  • $\begingroup$ But you didn't understand what I was saying. If $A$ is ill-conditioned or near singular, then expressions like $x^TA^{-1}x$ or $A^{-1}B$ will not make much sense because the values you get will be very sensitive to small perturbations in $x$ or $B$. The only exception would be if you can show that $x$ or $B$ have no components that lie in the directions of the eigenvectors corresponding to the very small eigenvalues of $A$. $\endgroup$ – Wolfgang Bangerth Nov 19 '12 at 12:45
  • $\begingroup$ I understand that if $A$ is singular we cannot invert it. In my case $A^{-1}$ is the estimate of a covariance matrix, which I know to exist. So when I get a badly conditioned $A$ I would like to perturb it so that I get an invertible matrix which is not too different from the original one. Which is more or less what Arnold is suggesting. The alternative it to increase my dataset until $A$ becomes invertible, but this is quite expensive. $\endgroup$ – Jugurtha Nov 19 '12 at 13:27
  • $\begingroup$ I understand what you're trying to do. But what I'm telling you is that if you have a near-singular $A$, then trying to compute $x^TA^{-1}x$ doesn't make any sense. You can paper over that by computing some $x^T\widetilde{A^{-1}}x$ but you're just hiding the problem: Yes, you can compute it; but no, it's not going to be a useful number you compute this way. $\endgroup$ – Wolfgang Bangerth Nov 19 '12 at 14:28
  • $\begingroup$ Ok, so in your opinion the only sensible thing to do if $A$ is almost singular is to get a larger dataset/repeat the simulation until I get an invertible $A$? $\endgroup$ – Jugurtha Nov 19 '12 at 14:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.