1
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ If you don't know anything about the structure of the matrices, this is the best you can do. $\endgroup$
    – Thomas
    Commented Nov 28, 2014 at 10:14
  • $\begingroup$ @sonystarmap Wouldn't $n^2$ be better then? :) $\endgroup$ Commented Nov 28, 2014 at 13:52

2 Answers 2

2
$\begingroup$

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)$.

$\endgroup$
4
  • $\begingroup$ +1 Good example, even though you might want to mention that $n_a=m_b$. $\endgroup$
    – Thomas
    Commented Nov 28, 2014 at 10:14
  • $\begingroup$ It is not really tight, as far as I know, you only showed that it is tight in one specific instance. Consider if $A$ and $B$ are completely full, then the bound is pretty loose, $nnz(AB) = m_a n_b$ but your bound is $m_a n_a m_b n_b$! $\endgroup$
    – the swine
    Commented Nov 28, 2014 at 11:23
  • $\begingroup$ I don't think that statement is true. How can you be (formally) proving that there is not a better bound just by showing a single case? If nothing else, $min(nnz(A) nnz(B), m_a n_b)$ is a better bound! $\endgroup$
    – the swine
    Commented Nov 28, 2014 at 12:22
  • $\begingroup$ I meant that there cannot be another function $b(A,B)$ such that $nnz(AB)\le b(A, B)$ and $b(A,B)<nnz(A)nnz(B)$ holds for all $A$ and $B$. There could be improvements of the bound for some $A$ and $B$, of course. Maybe you should specify, what you are looking for. $\endgroup$
    – daw
    Commented Nov 28, 2014 at 13:09
2
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Wow! What a great answer. I won't be able to get the results in Excel, this will have to be done in "C", I'll get back to accept it once I've done so. By any chance, would you have a more specific reference to PETSc (document or function name)? $\endgroup$
    – the swine
    Commented Nov 28, 2014 at 14:21
  • 1
    $\begingroup$ @theswine It seems that in PETSc this is done in the MatMatMult routine (I'm not sure how they do the estimate there or if the allocate data dynamically for each row or column which might be memory inefficient). $\endgroup$ Commented Nov 28, 2014 at 14:29
  • 1
    $\begingroup$ @theswine In the MatrixMatrix class of the EpetraExt package of Trilinos, there is a method Multiply. It uses an estimate from the ML package (an algebraic multigrid package in Trilinos) which looks rather like an alchemy to me. In my experience, it rather underestimates the actual number of nonzeros. $\endgroup$ Commented Nov 28, 2014 at 14:31

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .