Computing sparse matrix products into a dense result

Question

I need to assemble a matrix (in dense form, of moderate size, say dimension 1000) which is most easily expressed as the product of several (4) sparse matrices. These matrices are most easily expressed by their action (a matrix-vector multiplication routine), although I suppose it's not too hard to write them out explicitly.

My question is what is the best way for computing the product, keeping in mind that the result will wind up being stored as a dense matrix anyways. Clearly, expressing everything as a dense matrix and computing the products seems ridiculously inefficient, since the intermediates are still quite sparse.

More details: The actual form of the matrix is like $$ A = DBD^HC $$ where $D$ is a mesh adjacency matrix (only two nonzeros per row), and $B$ and $C$ are banded matrices with bandwidth under 10.

Answer 1

If you really want all of the entries of the dense matrix product, you could compute it column by column, by computing the action on columns of the identity matrix.

Answer 2

Since the question was asked a long time ago, I would give a detailed overview answer.

The first question is how the assembled matrix $A$ will be used later, as the construction of the matrix itself is rarely a goal. In the question, it is clearly defined that one needs to have a dense matrix in the end of the day. However, for some purposes, having a dense matrix instead of its factorization into several matrices might be detrimental and should be avoided. The simplest example would be a solution of the system of linear equations using an iterative solver.

Now, if the construction of the dense matrix is unavoidable and\or required by the further algorithms you have several options.

• As clipper pointed out, the entries of the dense matrix $A$ can be manually computed column-by-column by applying matrix-vector products to the columns of the identity matrix: $D(B(D^H(CI_n)))$, where $I_n$ is the $n$th column of the identity matrix. In such a way, you will be able to avoid forming dense intermediate results and take advantage of the sparsity of $D$, $B$, and $C$.
• Depending on the tools you are using, matrix library can already support sparse matrix-matrix products. For example, a commonly used template-based library Eigen will naturally support such products.
• Explicitly use some SparseBLAS implementation. For example, Intel MKL now supports subroutines mkl_?csrmultcsr and mkl_?csrmultd that compute the product of two sparse matrices and store the results in a sparse CSR or dense format, respectively.

Using some SparseBLAS implementation has an advantage of a well-optimized and debugged code for matrix operations, already built-in potential for parallelization. However, benchmarking it against an MVP-based solution would be required for the sparsity patterns and problem sizes you are interested in. Using Eigen-like library is the easiest options, though you will lack a lot of control and can lose some performance, especially if for those particular functions Eigen does not use external BLAS/LAPACK engines and relies only on its own implementations.