0
$\begingroup$

So I have this line of code that involves a matrix inversion

X = A @ B @ np.linalg.pinv(S)

$A$ is an $n$ by $n$ matrix, $B$ is an $n$ by $m$ matrix and $S$ is an $m$ by $m$ matrix. $m$ is smaller than $n$ but usually not orders of magnitude smaller. Usually $m \sim \frac{1}{2} n$. $S$ is a symmetrical positive definite matrix. How do I make this line of code run faster in python?

I realize I can do X = np.linalg.solve(S.T, (A@B).T).T

but is there more I can do given $S$ is a symmetrical positive definite matrix?

Improve this question
$\endgroup$
7
  • 2
    $\begingroup$ Does $S$ change between executions of this line of code, or is it a fixed matrix? If $S$ is symmetric and positive definite, then pinv(S) is simply $S^{-1}$. $\endgroup$
    – Brian Borchers
    Commented Jul 4 at 22:28
  • $\begingroup$ If the same matrix $S$ is going to be used repeatedly, then you can significantly speed up the solution process. $\endgroup$
    – Brian Borchers
    Commented Jul 5 at 0:13
  • 1
    $\begingroup$ What implementation of BLAS/LAPACK are you using with Python? np.show_config() can be helpful in figuring this out. A high performance BLAS/LAPACK (such as Intel's MKL) can be much faster than the default libraries. $\endgroup$
    – Brian Borchers
    Commented Jul 5 at 0:16
  • 1
    $\begingroup$ Are you matrices dense or sparse? $\endgroup$
    – nicoguaro ♦
    Commented Jul 5 at 8:14
  • 1
    $\begingroup$ OpenBLAS is the fastest of the open source BLAS libraries. MKL is somewhat faster but only runs on Intel (not AMD) processors. You could install MKL (licenses are free but it is not open source) and configure numpy to use it. $\endgroup$
    – Brian Borchers
    Commented Jul 5 at 14:24

1 Answer 1

Reset to default
2
$\begingroup$

Since S is SPD, you can (and should) use np.linalg.cholesky to get a cholesky factorization which will be ~2x faster than the LU factorization that np.linalg.solve uses.

Improve this answer
$\endgroup$
1
  • 3
    $\begingroup$ Good answer. Also, one can get the benefits of cholesky without the manual work by using the wrapper scipy.linalg.solve with assume_a='pos' (note that this is scipy and not numpy). $\endgroup$
    – Federico Poloni
    Commented Jul 5 at 10:32

Not the answer you're looking for? Browse other questions tagged or ask your own question.