# Efficiency of Repeated Sparse Matrix-Vector Products

I recently read Wolfgang's answer to the question found here and found myself wondering about a related followup question.

Assume you have two sparse matrices $A$ and $B$. You need to do the matrix-vector multiplication $a = Xb$ where $X=AB$. He states in his answer that forming X is much more intensive than just doing two multiplies in the form $a=A(Bb)$, which is for the most part intuitive.

My question is, what if you have many $b$'s? Let's assume you have $n$ different $b$'s to multiply. How large does $n$ have to be in relation to the size/sparsity of $A$ and $B$ to make forming $X$ worthwhile or does forming $X$ never become more efficient? If nothing can be said in general, are there specific cases where conclusions can be drawn?

-
You're going to have to say something about the sparsity pattern of the matrices before this question can have any definite answer. On one extreme, consider $A$ and $B$ as diagonal matrices, in which case you should explicitly form $AB$ if there are multiple right-hand sides. On the other hand, consider $A$ and $B$ to be arrowhead matrices (diagonal plus full last row and column), in which case $AB$ will be dense, and you should never form $AB$. A general answer would have to take into account the number of nonzeros in $AB$. – Jack Poulson Dec 13 '12 at 17:16
For argument's sake, let's say it has 10-100 non-zero elements per row, mostly near the diagonal, as if it were formed in some discretization of a 2D mesh. The product would have significant fill in relation to A or B, but not so much as to be considered dense. – Godric Seer Dec 13 '12 at 18:09

With your choice of fuzzy language in your comment, you are carefully skirting the answer to your question. The number of nonzeros in a matrix is a decent proxy for cost to apply it to a vector, since the cost is dominated by the time to load its entries from cache. Roughly speaking, if $$\mathrm{nnz}(AB) > \mathrm{nnz}(A) + \mathrm{nnz}(B),$$ then you should apply them separately. This is typically the case for PDE matrices and (even more so) for power-law graphs such as appear in network analysis.

If you are applying the matrix to many vectors at once, you get to amortize the cost to load the matrix entries. Once again, you can look at the performance model for a matrix of your chosen density and number of vectors to decide whether it is better to explicitly form the product. For square matrices, it rarely is.

One place where explicit matrix products are useful is in "Galerkin" coarse grid operators. There, the coarse grid operator $A_c = P^T A_f P$ resides in a much smaller space than the fine grid operator $A_f$, so the explicit product is much cheaper to apply.

-

As Jed already remarked, it all depends on the matrices. But you can get an idea by considering, for example, the 5-point stencil of the Laplace equation in 2d. There, each row of the stiffness or mass matrix $A$ has five entries. But the product of two such matrices has 13 entries (apply the five point stencil to every location of the five point stencil). In 3d, $A$ has 7 entries per row, but $A^2$ has 33. The situation becomes even worse if you have, say, a $Q_2-Q_1$ discretization of the Stokes equation in 3d.

As I mentioned in the other post, however, a significant cost in forming the product of a matrix is not actually the floating point operations but instead the cost of determining the sparsity pattern and resulting memory allocations. You will need to take that into account when forming a cost model.

-
Even though the the memory allocations are the main cost of the multiplications, it appears from both answers that for many cases the fill in such that even if the $AB$ multiplication was free, the matrix-vector multiplications would be less efficient when measures in raw flops. – Godric Seer Dec 14 '12 at 15:01
I suppose it depends on how many matrix-vector products you want to compute... – Wolfgang Bangerth Dec 15 '12 at 3:29