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

[Nathaniel Smith]

On Mon, Nov 8, 2010 at 12:00 PM, Joon <groups.and.lists at gmail.com> wrote:
> Another question is, is it better to do cho_solve(cho_factor(A), b) than
> solve(A, b)?

If A is symmetric positive definite, then using the cholesky
decomposition should be somewhat faster than using a more general
solver. (Because, basically, the cholesky decomposition routine
"knows" that your matrix is symmetric, so it only has to "look at"
half of it, while a generic solver routine has to "look at" your whole
matrix regardless). And indeed, that seems to be the case in numpy:

In [18]: A = np.random.normal(size=(500, 500))
In [19]: A = np.dot(A, A.T)
In [20]: b = np.random.normal(size=(500, 1))

In [21]: %timeit solve(A, b)
1 loops, best of 3: 147 ms per loop

In [22]: %timeit cho_solve(cho_factor(A), b)
10 loops, best of 3: 92.6 ms per loop

Also of note -- going via the inverse is much slower:

In [23]: %timeit dot(inv(A), b)
1 loops, best of 3: 419 ms per loop

(I didn't check what happens if you have to solve the same set of
equations many times with A fixed and b varying, but I would still use
cholesky for that case. Also, note that for solve(), cho_solve(),
etc., b can be a matrix, which lets you solve for many different b
vectors simultaneously.)

-- Nathaniel