I am working on a problem with very large sparse matrices. I'd like to compute $A^{-1} B$, that is a crucial part of converting DAE to ODE (and there is no workaround). Here size of $A$ is 2E+5 x 2E+5 with 0.7% density (especially all diagonal entries are non-zeros). And $B$ is 2E+5 x 1E+6 whose density is only 1E-6.

I tried to solve this problem in python by scipy.sparse.spsolve, which eventually gave me a memory issue. The same issue happened while using an iterative solver. Unless I save each column to the disk and free the memory before computing the next column in the solution, the memory issue will persist.

So does anybody know a memory-efficient sparse solver or suggest a better way to solve this problem?

  • $\begingroup$ Original question on StackOverflow $\endgroup$
    – Anton Menshov
    Jul 9, 2020 at 14:32
  • 3
    $\begingroup$ I think the question as asked is not really solvable. You have a nice sparse $\mathbf B$ that you are clobbering into a dense mess $\mathbf C = \mathbf A^{-1} \mathbf B$. You must pursue methods that do not require this .. something that can work with $\mathbf C$ in its present/operator form, as a cascade of sparse matvec and sparse backsolve. You might be able to gather ideas by posting a new question that presents the larger context. $\endgroup$ Jul 9, 2020 at 14:44
  • 3
    $\begingroup$ Can you explain why you need to compute this full matrix? I'm speculating you have an ODE/DAE of the form $A y' = f(y) = B y$. Many ODE/DAE solvers are specifically designed to handle mass matrices and maintain sparsity. Explicit methods will require a linear solve with $A$ per stage, which can be done with an iterative method or a precomputed sparse decomposition. Implicit methods often involve linear solves with matrices of the form $A - h \gamma J$, where $J$ is the Jacobian of $f$. As long as $J$ isn't dense, this should avoid these memory issues. $\endgroup$ Jul 9, 2020 at 15:17
  • $\begingroup$ Hi first thanks for all your replies. @StevenRoberts I would select some "coupled dofs" by investigating the non-zero pattern in this $A^{-1}B$ matrix. So that's why I need to know how resultant matrix looks like. $\endgroup$
    – Lambda Z
    Jul 13, 2020 at 9:05

1 Answer 1


$A^{-1}B$ is a $2\cdot 10^5 \times 1\cdot 10^6$ matrix that is likely to be full, so the size of your output is 1.6 TB. I don't think you have other alternatives than writing it to the disk.

  • 3
    $\begingroup$ Except, of course, to just never compute $A^{-1}B$ explicitly as a matrix :-) $\endgroup$ Jul 10, 2020 at 15:21

Your Answer

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

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