0
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • $\begingroup$ Could you give ballpark estimates for how large N, M, and the matrix size are? The effectiveness of different solutions may depend on this. $\endgroup$
    – Federico Poloni
    Jan 26, 2021 at 8:46
  • $\begingroup$ A ballpark estimate is M=50, N=100,000. $\endgroup$
    – Yiftach
    Jan 26, 2021 at 14:28
  • $\begingroup$ And the matrix size? $\endgroup$
    – Federico Poloni
    Jan 26, 2021 at 15:53
  • $\begingroup$ Oh, sorry. Let us say 64x64. $\endgroup$
    – Yiftach
    Jan 26, 2021 at 17:12
  • $\begingroup$ 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. $\endgroup$
    – Federico Poloni
    Jan 26, 2021 at 17:20

1 Answer 1

Reset to default
2
$\begingroup$

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

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Yiftach
    Jan 25, 2021 at 21:39
  • $\begingroup$ @Yi What do you mean by "broadcasting"? $\endgroup$
    – Federico Poloni
    Jan 26, 2021 at 7:18
  • $\begingroup$ 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. $\endgroup$
    – Yiftach
    Jan 26, 2021 at 8:18
  • $\begingroup$ 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. $\endgroup$
    – Federico Poloni
    Jan 26, 2021 at 8:27
  • $\begingroup$ 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. $\endgroup$
    – Yiftach
    Jan 26, 2021 at 13:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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