# 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. Commented Jan 26, 2021 at 8:46
• A ballpark estimate is M=50, N=100,000. Commented Jan 26, 2021 at 14:28
• And the matrix size? Commented Jan 26, 2021 at 15:53
• Oh, sorry. Let us say 64x64. Commented Jan 26, 2021 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. Commented Jan 26, 2021 at 17:20

## 1 Answer

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.)

• I know that numpy is templated C that unfortunately calls GESV directly, and that scipy's version can call the positive definite versions, but do not support broadcasting. However, I'm looking for the best purely pythonic solution. Commented Jan 25, 2021 at 21:39
• @Yi What do you mean by "broadcasting"? Commented Jan 26, 2021 at 7:18
• Simultaneously solve for N matrices- For matrix 1 solve for vectors 1...M, for matrix 2 solve for vectors M+1...2M, etc. Is that not the correct term here? Numpy's solve documentation makes me believe that this is a feature they support and scipy does not. Commented Jan 26, 2021 at 8:18
• 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. Commented Jan 26, 2021 at 8:27
• Of course. I'm okay with numpy's way of "cheating", it's just that I want my code to work regardless of operating system/architecture (e.g on both linux/windows) so writing my own C is not an option in this case. Commented Jan 26, 2021 at 13:32