title: Inverting ill-conditioned matrix in NumPy
body: <p>I am trying to invert a covariance matrix to calculate a Mahalanobis-distance. However, my covariance matrix is very ill-conditioned (condition number > 10^17). This apparently leads to large rounding errors, resulting in a "matrix inverse" that actually does not have much to do with the real inverse.</p>
  
  <p>My code is:</p>
  
  <pre><code>sample = np.array([[523.77, 468.3, 491.55], [234.96, 221.17, 228.43], [342.02, 337.87, 351.31]], np.float64)
  
  print "Sample:"
  print sample
  print
  print "Sample Inverse:"
  print np.linalg.inv(sample)
  print
  print "Sample x Sample Inverse:"
  print np.dot(sample, np.linalg.inv(sample))
  print
  print
  print "Cov(Sample):"
  print np.cov(sample)
  print "Determinant", np.linalg.det(np.cov(sample))
  print "Condition Number", np.linalg.cond(np.cov(sample))
  print
  print "Cov(Sample) Inverse:"
  print np.linalg.inv(np.cov(sample))
  print
  print "Cov(Sample) x Cov(Sample) Inverse:"
  print np.dot(np.cov(sample), np.linalg.inv(np.cov(sample)))
  </code></pre>
  
  <p>Which gives the following result:</p>
  
  <pre><code>Sample:
  [[ 523.77  468.3   491.55]
   [ 234.96  221.17  228.43]
   [ 342.02  337.87  351.31]]
  
  Sample Inverse:
  [[ 0.01204863  0.03620989 -0.04040286]
   [-0.10240561  0.36837033 -0.09623795]
   [ 0.08675788 -0.38953001  0.13473714]]
  
  Sample x Sample Inverse:
  [[  1.00000000e+00   0.00000000e+00   0.00000000e+00]
   [  3.55271368e-15   1.00000000e+00   3.55271368e-15]
   [  0.00000000e+00   0.00000000e+00   1.00000000e+00]]
  # This is obviously fine
  
  
  Cov(Sample):
  [[ 775.9353      190.68715      40.55945   ]
   [ 190.68715      47.58543333   15.68986667]
   [  40.55945      15.68986667   47.36003333]]
  Determinant -1.3802775949e-10
  Condition Number 3.45634882629e+17
  
  Cov(Sample) Inverse:
  [[ -1.45440004e+13   6.08180372e+13  -7.69277831e+12]
   [  6.08180372e+13  -2.54320238e+14   3.21685688e+13]
   [ -7.69277831e+12   3.21685688e+13  -4.06895190e+12]]
  
  Cov(Sample) x Cov(Sample) Inverse:
  [[-0.9375    -1.         0.09375  ]
   [-0.4375    -0.125      0.1953125]
   [ 0.        -0.25       1.0625   ]]
  # This is wrong
  </code></pre>
  
  <p>which is obviously not what I want...</p>
  
  <p>In reality, my covariance matrix is also not a 3x3, but rather a 145x145 matrix, which makes things worse: The condition number is now larger than 10^18, and numpy.linalg.det() returns 0 for the determinant (which would mean that no inverse exists).</p>
  
  <p>Any suggestions how this problem might be fixed are very much appreciated. Thanks in advance for all your efforts! :)</p>