# Solve for large array of PD matrices

I have N matrices that are positive definite, and I have to solve for a M vectors. As M is large in my case, doing all solves simultaneously using np.linalg.solve burdens my RAM and sometimes not possible. However, splitting to batches using and solving on them unnecessarily performs the factorization step multiple times and does not cache it. Both options do not leverage the fact that the matrices are positive definite.

What is the best course of action, in python, for solving for all vectors?

• Could you give ballpark estimates for how large N, M, and the matrix size are? The effectiveness of different solutions may depend on this. – Federico Poloni Jan 26 at 8:46
• A ballpark estimate is M=50, N=100,000. – Yiftach Jan 26 at 14:28
• And the matrix size? – Federico Poloni Jan 26 at 15:53
• Oh, sorry. Let us say 64x64. – Yiftach Jan 26 at 17:12
• OK, so the only large dimension is N; the other sizes look very manageable. I am afraid this boils down to the well-known "for loops are slow in Python" problem. Maybe calling linalg.solve or the dposv Python binding in a Numba or Cython loop would work. Note that you can solve the M systems with one single Lapack call by sticking them as columns in a matrix. – Federico Poloni Jan 26 at 17:20

I would just compute the Cholesky factorization and then solve in batches using it. This will get technical, though: you will need to call Lapack functions by hand, I am afraid (*potrf and *potrs), since Python doesn't help you here, so to use the exact same algorithm you may want to check how it is done in the source of linalg.solve and dposv.f (good luck with the Fortran).
Also, your go-to function in these cases is scipy.linalg.solve; it has options to exploit symmetry and positive definiteness, unlike its numpy counterpart. (Both numpy and scipy have a linalg.solve function, which accept different arguments and yes, I agree that it's confusing.)
• Oh I see now, thanks. You are correct, that seems a valid reason to use np.linalg.solve rather than its scipy counterpart. If you are in a setting where you gain a lot by 'vectorizing' then I don't know of a pure-Python solution. Note that even np.linalg.solve "cheats" and uses C internally. – Federico Poloni Jan 26 at 8:27