Why is inverting a positive definite matrix via Cholesky decomposition slower than regular inversion with numpy?

Question

It's mathematically known that inverting a positive definite matrix via Cholesky decomposition is faster than just using np.linalg.inv(X). However, when I experimented with both and it turns out Cholesky decomposition's performance is worse!

# Inversion through Cholesky
p = X.shape[0]
Ip = np.eye(p)
%timeit scipy.linalg.cho_solve(scipy.linalg.cho_factor(X,lower=True), Ip)

The slowest run took 17.96 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 107 µs per loop


# Simple inversion
%timeit np.linalg.inv(X)

The slowest run took 58.81 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 25.9 µs per loop

The latter took shorter. Why is this? In R, chol2inv(chol(X)) is usually faster than solve(X).

Answer 1

I run a comparison on 1000x1000 matrices, and the Inversion through Cholesky was roughly twice as fast.

Answer 2

Perhaps your matrix is too small. I just tested matrix inversion for a $2\times2$ matrix in Matlab using Cholesky decomposition followed by LU decomposition. 999999 repeats take 5 seconds using Cholesky and only takes 3.4 seconds using LU. Indeed the algorithm of Cholesky followed by back substitution has a smaller big O result, but the result is an asymptotic result and only applies for big matrices.

#LU decomposition
tic
for i=1:999999
    (V_i(:,:,2)+[0 1e-10;0 0])\eye(2);
end
toc
Elapsed time is 3.430676 seconds.

#Cholesky
tic
for i=1:999999
    (V_i(:,:,2))\eye(2);
end
toc
Elapsed time is 4.824175 seconds.