Solving large-scale sparse matrix equation
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.
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!