efficiency of inverting a matrix in numpy with Cholesky decomposition

Question

I have a symmetric positive-definite matrix (e.g. Covariance matrix), and I want to calculate its inverse. In math, I know that it is more efficient to use Cholesky decomposition to invert the matrix, especially if your matrix is big. But I was not sure how does "numpy.lianlg.inv()" works. Say I have the following code:

import numpy as np

X = np.arange(10000).reshape(100,100)
X = X + X.T - np.diag(X.diagonal()) #  symmetry 
X = np.dot(X,X.T) # positive-definite

# simple inversion:
inverse1 = np.linalg.inv(X) 

# Cholesky decomposition inversion:
c = np.linalg.inv(np.linalg.cholesky(X))
inverse2 = np.dot(c.T,c)

Which one is more efficient (inverse1 or inverse2)? If the second one is more efficient, why is numpy.linalg.inv() not using this instead?

Regarding your last question - numpy does not know that your matrix is symmetric, so cannot use the latter method. Checking if the matrix is symmetric would be slow. — Eric
– Eric, Commented Jun 3, 2017 at 21:03
Note that also, inv does not take advantage of the fact that cholesky is triangular — it does not use lapack's DTRTRI — Eric
– Eric, Commented Jun 3, 2017 at 21:06
Your code doesn't run for me, claiming that X is not positive definite (presumably due to overflow) — Eric
– Eric, Commented Jun 3, 2017 at 21:10
In general, it's bad idea to invert a matrix. inv is expensive and isn't numerically stable. Usually, you want to multiply the inverse with a vector, i.e., you want to solve a system of equations. In all such cases, it's better to just solve the system using something like linalg.solve (telling solve that the matrix is symmetric and positive definite will make solve use Cholesky). If you want to use the inverse multiple times, compute and store the decomposition for later use. — Praveen
– Praveen, Commented Jun 4, 2017 at 0:09
if you do want to invert the cholesky factor, use scipy.linalg.lapack.dtrtri — Bananach
– Bananach, Commented Sep 23, 2020 at 15:30

Eric · Accepted Answer · 2017-06-04 00:00:00Z

2

With the following setup:

import numpy as np

N = 100
X = np.linspace(0, 1, N*N).reshape(N, N)
X = 0.5*(X + X.T) + np.eye(N) * N

I get the following timings with IPython's %timeit:

In [28]: %timeit np.linalg.inv(X)
255 µs ± 30.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [29]: %timeit c = np.linalg.inv(np.linalg.cholesky(X)); np.dot(c.T,c)
414 µs ± 15.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

edited Jun 4, 2017 at 0:00

answered Jun 3, 2017 at 21:17

Eric

98.1k54 gold badges257 silver badges389 bronze badges

Sign up to request clarification or add additional context in comments.

Collectives™ on Stack Overflow

efficiency of inverting a matrix in numpy with Cholesky decomposition

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related