[Numpy-discussion] Solving Ax = b: inverse vs cholesky factorization

Bruce Southey bsouthey@gmail....
Mon Nov 8 14:06:03 CST 2010


On 11/08/2010 01:38 PM, Joon wrote:
> On Mon, 08 Nov 2010 13:23:46 -0600, Pauli Virtanen <pav@iki.fi> wrote:
>
> > Mon, 08 Nov 2010 13:17:11 -0600, Joon wrote:
> >> I was wondering when it is better to store cholesky factor and use 
> it to
> >> solve Ax = b, instead of storing the inverse of A. (A is a symmetric,
> >> positive-definite matrix.)
> >>
> >> Even in the repeated case, if I have the inverse of A (invA) stored,
> >> then I can solve Ax = b_i, i = 1, ... , n, by x = dot(invA, b_i). Is
> >> dot(invA, b_i) slower than cho_solve(cho_factor, b_i)?
> >
> > Not necessarily slower, but it contains more numerical error.
> >
> > http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/
> >
> >> I heard calculating the inverse is not recommended, but my 
> understanding
> >> is that numpy.linalg.inv actually solves Ax = I instead of literally
> >> calculating the inverse of A. It would be great if I can get some
> >> intuition about this.
> >
> > That's the same thing as computing the inverse matrix.
> >
>
> Oh I see. So I guess in invA = solve(Ax, I) and then x = dot(invA, b) 
> case, there are more places where numerical errors occur, than just 
> x = solve(Ax, b) case.
>
> Thank you,
> Joon
>
>
>
Numpy uses SVD to get the (pseudo) inverse, which is usually very 
accurate at getting (pseudo) inverse.

There are a lot of misconceptions involved but ultimately it comes down 
to two options:
If you need the inverse (like standard errors) then everything else is 
rather moot.
If you are just solving a system then there are better numerical solvers 
available in speed and accuracy than using inverse or similar approaches.

Bruce




More information about the NumPy-Discussion mailing list