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?
dposvPython 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$