20

I have a badly conditioned matrix, whose rcond() is close to zero, and therefore, the inverse of that matrix does not come out to be correct. I have tried using pinv() but that does not solve the problem. This is how I am taking the inverse:

X = (A)\(b);

I looked up for a solution to this problem and found this link (last solution) for improving the matrix. The solution there suggests to use this:

A_new = A_old + c*eye(size(A_old));

Where c > 0. So far employing this technique works in making the matrix A better conditioned and the resultant solution looks better. However, I investigated using different values of c and the resultant solution depends on the value of chosen c.

Other than manually investigating for the value of c, is there an automatic way through which I can find the value of c for which I get the best solution?

Amro
  • 123,847
  • 25
  • 243
  • 454

3 Answers3

16

In discrete inverse theory, adding a small value c to the diagonal of the matrix A about to be inverted, is called damping the inversion and the small value to be added c is called Marquardt-Levenberg coefficient. Sometimes matrix A has zero or close to zero eigenvalues, due to which the matrix becomes singular; adding a small damping coefficient to the diagonal elements makes it stable. Bigger is the value of c, bigger is the damping, your matrix inversion is more stable, but you are further away from true solution. Smaller is the value of c, smaller is the damping, closer is your inverted matrix to true inverted matrix but it might become unstable. One 'adaptive damping' technique sometimes used is - start with a test value of c, invert the matrix A, then decrease the value of c, do the inversion again and so on. stop when you get weird values in inverted matrix due to A becoming singular again, like really large numbers. I think this doesn't exactly answer your question, but it was too long to put it in a comment.

Guddu
  • 2,325
  • 2
  • 18
  • 23
16

As already pointed out in the comments, the answer to your question crucially depends on your application. Maybe adding a small multiple of the identity matrix is the right thing to do, maybe not. In order to determine that you need to tell us: How did this matrix arise? And what do you need the inverse for?

Two common cases are:

  • If you know the matrix A exactly, e.g. because it is the design matrix in a general linear model b = A * X, then modifying it is not a good idea. In this case, the matrix defines a linear system of equations, and if the matrix is singular this means there is no unique solution to this system. To pick one from the infinitely many possible solutions, there are different strategies: X = A \ b picks a solution with as many zero coefficients as possible, while X = pinv(A) * b picks the solution with the minimum L2 norm. See the examples in the documentation of pinv.

  • If the matrix A is estimated from data, e.g. a covariance matrix for an LDA classifier, and you have reason to believe the true value is not singular and the singularity is just due to not having enough data points for the estimate, then applying regularization or "shrinkage" by adding a small multiple of the identity matrix is a common strategy. In this case, Schäfer and Strimmer (2005) describe a way to estimate the optimal regularization coefficient from the data itself.

But I'm sure there are other cases with other answers, too.

Ben Jackson
  • 90,079
  • 9
  • 98
  • 150
A. Donda
  • 8,381
  • 2
  • 20
  • 49
8

Adding small values to the diagonal of A is roughly equivalent to introducing an L2-norm regularization term in the least squares problem Ax=b. That is one seeks to minimize the residual as well the added constraint:

min ||Ax-b||^2 + lambda*||x||^2

where lamdba controls the weight given to minimizing the constraint vs. minimizing the residual norm.

Usually this parameter is chosen using some kind of cross-validation technique.

Amro
  • 123,847
  • 25
  • 243
  • 454
  • So, essentially there's nothing wrong in adding small values to the diagonal of `A`, right? –  Sep 12 '13 at 16:10
  • 1
    yes, to be exact it's adding a multiplier of the identity matrix like you did: `c*eye(.)`. There is a complete [manual](http://www2.imm.dtu.dk/~pcha/Regutools/RTv4manual.pdf) accompanying the [package](http://www.mathworks.com/matlabcentral/fileexchange/52-regtools) I mentioned before. – Amro Sep 12 '13 at 21:21
  • Considering a 2×2 `A` of rank 1, but not otherwise degenerate, visualise a general example. The first term becomes a parabolic-cross-section valley running infinitely far across the plane, and offset from the centre somewhat. The second term is a paraboloid centred on the origin. Before adding the second term, there wasn't a unique solution. But having added it, we've added gradient to the valley floor (unless it crosses the origin), and now instead of getting at least one correct solution out of the infinite set of them, we're getting a unique one that isn't quite correct. – Evgeni Sergeev Mar 08 '16 at 13:50
  • 1
    hmm.. I think @A.Donda made some good points above (+1); It really depends on the application, and nudging the matrix this way only makes sense if it was somehow estimated and came out badly conditioned because of insufficient data during the estimation... When solving a system in general and the matrix is singular, then that's that. IMO If you modify it then you are really solving a slightly different system :) – Amro Mar 09 '16 at 13:32