Solving large-scale sparse matrix equation

Question

I am doing a sparse-matrix multiplication A^-1@B in Python within scipy.sparse. Both A and B are stored in csc format. The size of A and B is around 1E+6 x 2E+5 and 1E+6 x 1E+6. My intuition tells me I shouldn't invert A but use scipy.sparse.spsolve alternatively. But spsolve(A,B) gives me memory issue. I also try iterative solver like gmres and cg. But since the result from such an iterative solver is usually dense, the memory issue still shows up. Does anybody know if there is an appropriate way to solve this large-scale problem? Many thanks to you guys!

Answer 1

Several things to note:

1. An inverse of the sparse matrix is not necessarily sparse, and generally does not retain the same sparsity pattern, unless you enforce it (with a degree of approximation). See this or this as sample references on the topic of sparse approximate inverse.

2. Multiplication of two sparse matrices is also not the best operation, the sparsity pattern generally also is not preserved.

3. When I mention that the sparsity pattern is not preserved, it usually implies that the sparsity of the problem decreases, thus requiring more and more memory.

So, having only the information you provided, these are the following possibilities I can think of:

1. The resultant multiplication of A^{-1}*B in the exact computations lost most of its sparsity due to the aforementioned reasons. Therefore, you might be trying to compute something that can no longer be stored in memory (which is quite a possibility, given the row and column dimensions of the matrices).

2. You are mentioning solving a linear system Ax=B, so, I assume that is your original goal. Instead of the usual single-RHS solve of Ax=b, where b is a general vector, you want to solve a multi-RHS system AX=B, where B is a matrix, luckily, B is sparse.

   You might consider solving for each nth column of B (denoted as B_n) separately using an iterative solver Ax_n=B_n, and then trying to clean each x_n based on some tolerance (with 0 being an extreme exact tolerance). If each of your x_n is dense, then you are back to square 1. The solution to your system AX=B denoted as X is not sparse, given the tolerance/sparsification criterion of your choice.

In general, I would suggest you to evaluate if the result you are looking for can/should/must/is sparse. If not, will you be satisfied with a certain approximate result.

If you are interested in a more detailed discussion of your problem, consider Computational Science SE as a potential community. However, it would be imperative to post some details about your problem, sparsity patterns, maybe the physical problem where A and B come from.