Solving large-scale sparse matrix equation

user_13880382 0
Viewed 369

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!

numerical-methods scientific-computing python
2 Answers

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.

user_24528690

It is possible to multiply large matrices using a deep neural network which is designed to predict the result of essentially a matrix factorization-like model. Piecemeal computation is the benefit, making the memory size not a problem. Your labels for supervised will be computed dynamically because you can predict perfectly the exact result of each dot product. Maybe no learning per se is needed, not sure, because your problem is actually solved as soon as you compute the labels. Tensorflow2 is good for this.

There is probably something COTS that already does what you need if you look more, so you probably won't need to go through the effort of making the thing I described.

user_391230
Related