nnz-preserving sparse matrix multiplication

Question

Let A and B be sparse matrices in $ \mathbb{R}^{m\times m}$ with (roughly) the same density $p$. I want to efficiently compute a matrix $C$ that in some sense is "closest" to $AB$ while still having density $p$. "Closest" should probably mean "minimizes $||C - AB||_F$" but I'm open to other interpretations.

If $\textrm{nnz}$ is the number of non-zero entries in A (= the number of non-zero entries in B), the main computational requirement is that the operation should only need $O(\textrm{nnz})$ memory.

Question 1: this seems like a problem that would have a lot of research already, but I failed to find much. Are there any references/keywords I could use when googling?

---

My naive idea:

1. do a sparse matmul without actually writing the resulting entries in memory. Instead, feed them into an on-line order selection algorithm to choose the $\textrm{nnz}$'th biggest-norm element in the output stream.
2. Re-do the sparse matmul, but only record entries with norm greater or equal than the threshold found in step 1.

This seems to straightforwardly minimize the Frobenius norm but is also kind of inefficient: asymptotically the runtime can still be $O(m^2)$ for bad sparsity patterns, and practically, we're essentially doing the multiplication twice.

Also, it's unclear what to do with a matrix whose non-zero entries are all identical. Maybe just take the first $nnz$ entries? But that makes analysis much harder.

Question 2: are there any obvious improvements to this?

Question 3: what's the dynamic behavior of this algorithm, under the assumption that all non-zero entries are distinct? E.g. let $\odot$ be this density-preserving multiplication and $M_i$ be sparse matrices. Can we say anything about $||(M_1 M_2 \ldots M_n) - (M_1 \odot M_2 \odot \ldots \odot M_n)||_F$ ? Probably not for pathological cases, but what about if $M_i$ are iid random sparse matrices?

Answer 1

I do not know anything about the specific problem, but you might want to look into the techniques the "SParse Approximate Inverse (SPAI)" community has come up with over the past two decades. There, one is looking for a matrix $B \approx A^{-1}$ so that $\|BA-I\|$ is minimized (with regard to some norm, typically the Frobenius norm) and requiring that $B$ is sparse. The question then is how one should choose the sparsity pattern of $B$ so that $BA$ is approximately equal to the identity -- a problem not so far removed from what you are looking for.

I can't point to any specific paper (not knowing the area well), but if you search for SPAI, you should find literature on the topic.

Answer 2

One reason why there might not be much research on this is that one usually avoids sparse matrix multiplications as much as possible in the first place. When applying the product to a vector $ABv$, you can associate from the right $A(Bv)$, and when solving linear systems you can add auxiliary variables to avoid forming the product; for instance, turn $c = ABx$ into $y=Bx,\, c=Ay$, and solve for $\begin{bmatrix}x\\y\end{bmatrix}$.

In fact, my first thought after reading your question is that this is probably an "XY problem", and there might be a better way to do the task you need this approximate matmul for.

Anyway, a trick to avoid doing the multiplication twice is storing the computed $(i,j,C_{ij})$ into a priority queue, and deleting from it the entry with the smallest norm as soon as a larger one arrives.