# Sparse matrices dot multiplication

I have two huge sparse matrices, and I am interested in dot multiplication them. But I have a memory issue so a notebook is a shutdown. I am using the 11.3 version. Any suggestion on how to avoid memory issues.

Edit:

n = 30000;
m = 22000000;
mat = SparseArray[
RandomInteger[{1, n}, {m, 2}] -> RandomReal[{0, 1}, m], {n, n}];
mat2 = SparseArray[
mat["NonzeroPositions"] ->
RandomReal[{0, 1}, Length@mat["NonzeroPositions"]],
Dimensions[mat]];

• Should the actual matrices be too large, one has to avoid matrix-matrix multiplication to to rephrase everything into matrix-vector mutiplications. For example, u being a vector, instead of res.u, you should use mat.(mat.u). This might bring also a considerable speed-improvement, if only few of results of the form res.u are needed. I gets more complicated when res has to enter a linear system which you want to solve. But if you cannot build res, then you cannot solve the system with direct methods anyways. So you would have to switch to (preconditoned!) iterative methods. – Henrik Schumacher Sep 3 '18 at 21:04
• ... and for these iterative methods (e.g. GMRES), you need only matrix-vector multiplication. So also in that case, you can go with mat.(mat.u). – Henrik Schumacher Sep 3 '18 at 21:05
• Are you sure that both matrices that you dot-multiply, are sparse? If one is dense, then the product will be coerced to a sense matrix. If the dense matrix is of size $k \times n$ and the sparse matrix is of size $n \times m$ with $k n \ll n m$, then the product will be a dense $k \times m$ matrix and may not fit into memory although each factor fits. – Henrik Schumacher Sep 3 '18 at 21:20
• If Mathematica is repeatedly crashing, it might be worth restarting your computer. – mikado Sep 3 '18 at 21:47
• Dot products of two sparse matrices (unless they are very sparse) will result in a SparseArray representation of a dense matrix (this will cause a memory issue). Division on the other hand has a divide by 0 issue, but should not have a memory issue. – Carl Woll Sep 5 '18 at 17:06

More an extended comment than an answer.

Many of us are used to obtain a quite sparse matrix from the matrix product of two sparse matrices. That's is maybe because we are used to obtain sparse arrays as the system matrices of differential equations where each degree of freedom interacts with only few other degrees of freedoms or as adjacency matrices of more or less random graphs. But actually, it suffices to have one degree of freedom in the system (or a vertex that is connected to all other vertices in the graph) in order to obtain a completely dense matrix from a multiplication of a matrix with itself as the following example shows:

n = 10000;
A = SparseArray[{}, {n, n}];
A[[1, All]] = 1;
A[[All, 1]] = 1;
A["Density"]
B = A.A;
B["Density"]


0.00019999

1.

• This is in part because no sparsification is done along the way, e.g. To see if a different background value might work. But your point remains valid, and random values would force a dense matrix there. For many purposes, it can be better to defer matrix-times-matrix altogether, and instead iterate matrix-times-vector. – Daniel Lichtblau Sep 5 '18 at 19:48