39

I do not quite understand why numpy.linalg.solve() gives the more precise answer, whereas numpy.linalg.inv() breaks down somewhat, giving (what I believe are) estimates.

For a concrete example, I am solving the equation C^{-1} * d where C denotes a matrix, and d is a vector-array. For the sake of discussion, the dimensions of C are shape (1000,1000) and d is shape (1,1000).

numpy.linalg.solve(A, b) solves the equation A*x=b for x, i.e. x = A^{-1} * b. Therefore, I could either solve this equation by

(1)

inverse = numpy.linalg.inv(C)
result = inverse * d

or (2)

numpy.linalg.solve(C, d)

Method (2) gives far more precise results. Why is this?

What exactly is happening such that one "works better" than the other?

ali_m
  • 71,714
  • 23
  • 223
  • 298
ShanZhengYang
  • 16,511
  • 49
  • 132
  • 234
  • 12
    Best pick up a textbook. Golub and van Loan would be a good choice. – ev-br Jul 06 '15 at 22:01
  • Just a guess ( possibly from somewhere of Gilbert Strang Books) : I think the inverse is never really calculated during the call of the underlying LAPACK routine. What if you try pinv() instead of inv() ? Look at the reference of pinv() – Moritz Jul 06 '15 at 22:01
  • 3
    In most cases, linear equation solution is faster and produces less round-off error than inversion and multiplication. Fundamentally, these are operations on "floating" number representations, not the real numbers themselves. – Frank M Jul 06 '15 at 22:13
  • Would be interesting to knoe the comparison of pinv wth solve – Moritz Jul 07 '15 at 05:55
  • 3
    You almost *never* want to compute the inverse of a matrix... it's inefficient, and the error is generally higher than other methods. – Bakuriu Oct 03 '16 at 07:14
  • @Bakuriu Let's not go through this again :) Yes, I realize that you usually don't want to, but there are times when you would. Many of the algorithms I can think of come from computational statistics, when there is no other option available. – ShanZhengYang Oct 03 '16 at 07:19

1 Answers1

60

np.linalg.solve(A, b) does not compute the inverse of A. Instead it calls one of the gesv LAPACK routines, which first factorizes A using LU decomposition, then solves for x using forward and backward substitution (see here).

np.linalg.inv uses the same method to compute the inverse of A by solving for A-1 in A·A-1 = I where I is the identity*. The factorization step is exactly the same as above, but it takes more floating point operations to solve for A-1 (an n×n matrix) than for x (an n-long vector). Additionally, if you then wanted to obtain x via the identity A-1·b = x then the extra matrix multiplication would incur yet more floating point operations, and therefore slower performance and more numerical error.

There's no need for the intermediate step of computing A-1 - it is faster and more accurate to obtain x directly.

* The relevant bit of source for inv is here. Unfortunately it's a bit tricky to understand since it's templated C. The important thing to note is that an identity matrix is being passed to the LAPACK solver as parameter B.

ali_m
  • 71,714
  • 23
  • 223
  • 298