Does anyone know a good reference for bounds on sparse matrix multiplication? I'm interested in bounds of the number of scalar products required and bounds of the sparsity of the product.
I know that the upper bound of the sparsity of the product $C = AB$ would be $nnz(A)\cdot nnz(B)$, lower bound being zero (if matrix $A$ only has elements in columns, corresponding to rows where $B$ does not have elements), but I suspect that it is not a very tight bound.
On
Let $A$ be $m\times n$ and $B$ be $n\times p$ (assuming both sparse) and let $e_i$ be the $i$th column of the identity matrix of an appropriate size.
If the matrices are stored by rows, the product $C:=AB$ should be computed by rows as well. We have $$\tag{1} e_i^TC=\sum_j\underbrace{(e_i^TAe_j)}_{a_{ij}}e_j^TB $$ so the $i$th row of $C$ is given by linear combinations of the rows of $B$ with the coefficients given by the $i$th row of $A$. Let $$ r_A:=\max_{1\leq i\leq m}\mathrm{nnz}(e_i^TA), \quad r_B:=\max_{1\leq i\leq n}\mathrm{nnz}(e_i^TB) $$ be the maximal numbers of nonzeros per rows of $A$ and $B$, respectively. Since in (1), we take at most $r_A$ combinations of rows of $B$ of the "size" at most $r_B$, it follows that at most $r_Ar_B$ entries will be created in the row $i$ of $C$. Hence $$\tag{2} \mathrm{nnz}(C)\leq m\,r_A\, r_B. $$
Similarly, if $A$ and $B$ are stored by columns, the product is usually computed column by column using $$\tag{3} Ce_j=\sum_iAe_i\underbrace{(e_i^TBe_j)}_{b_{ij}}, $$ that is, the $j$th column of $C$ is given by the combination of the columns of $A$ with the coefficients in the column $j$ of $B$. With $$ c_A:=\max_{1\leq j\leq n}\mathrm{nnz}(Ae_j), \quad c_B:=\max_{1\leq j\leq p}\mathrm{nnz}(Be_j), $$ (3) indicates that $$\tag{4} \mathrm{nnz}(C)\leq p\,c_A\,c_B. $$
We can combine both (2) and (4) to $$ \mathrm{nnz}(C)\leq\min\{m\,r_A\,r_B,\;p\,c_A\,c_B\}. $$
Note that if your matrices are stored by rows or columns, you can use (1) or (3) to improve the estimate of $\mathrm{nnz}(C)$ for the price of one matrix-vector-like procedure (instead of taking maxima $r_A$, $r_B$ or $c_A$, $c_A$, you consider the actual numbers of entries). For example, using (1), for computing the $i$th row we take $\mathrm{nnz}(e_i^TA)$ combinations of sparse vectors of lengths $\mathrm{nnz}(e_j^TB)$, where $j$ runs over the nonzero pattern of the $i$th row of $A$. That is, $$ \mathrm{nnz}(e_i^TC)\leq\sum_{j:\;a_{ij}\neq 0}\mathrm{nnz}(e_j^TB) $$ and $$ \mathrm{nnz}(C)\leq\sum_i\sum_{j:\;a_{ij}\neq 0}\mathrm{nnz}(e_j^TB). $$ Similar estimate can be made if $A$ and $B$ are stored by columns using (3).
However, the exact number of nonzeros cannot be computed without a more expensive analysis of the combinations made in (1) and (3). Some heuristics can be made based on some assumptions of the "overlap" of the entries in $A$ and $B$ but this strongly depends on the applications.
I'm using (2) in my code and have good experience with it. You can also consult other sparse linear algebra packages such as Trilinos or PetSc for other approaches.
The bound $nnz(AB)=nnz(A)nnz(B)$ is tight: Take $$ A = \pmatrix{1 & 0 & \dots & 0 \\ \vdots & \vdots & & \vdots \\1 & 0 & \dots & 0 }\in \mathbb R^{m_a,n_a}, \quad B = \pmatrix{ 1 & \dots & 1 \\0 & \dots & 0 \\\vdots & &\vdots \\0 & \dots & 0 \\}\in \mathbb R^{m_b,n_b}. $$ Then $AB$ is a full matrix with all entries equal to one, $nnz(AB) = m_an_b=nnz(A)nnz(B)$.