2

When you're multiplying sparse matrices against other sparse matrices or dense matrices, what is the conventional approach for each? How are the sparse matrices stored? What does matrix multiplication work?

My understanding is there are several ways it can be stored:

(1) Store the non-zero elements in a hash table where the key is the 2D indices (mapped to 1D usually) and the value is the value of the non-zero element at that index.

(2) Store the non-zero elements in a $k$ sized array where $k$ is the number of non-zero elements. Each element stores the indices and the value at that index.

(3) Similar to (2), but you have a $n$ sized array where $n$ is the size of the original matrix. For each element of the array you store a list of non-zero elements.

I'm not sure which of these or some other formulation is used in practice and under what circumstances.

An additional question I have is do we know how this is done in linear algebra packages like BLAS and LAPACK?

anonuser01
  • 291
  • 1
  • 7
  • BLAS and LAPACK don't support sparse matrix operations. It is difficult to get good performance in sparse matrix-matrix multiplication (SPMM)- you can find many papers on the topic. – Brian Borchers Aug 23 '21 at 02:25
  • @BrianBorchers I'll do some literature review on this. Is there a quick 1-3 sentence summary you could give on why good performance is difficult with SPMM? Does BLAS and LAPACK just default to everything as if it were a dense matrix operation – anonuser01 Aug 23 '21 at 02:42
  • 1
    @anonuser01, there is no interface to input a sparse matrix into BLAS/LAPACK. Sparse matrix operations are simply not implemented in those. – Abdullah Ali Sivas Aug 23 '21 at 02:49
  • BLAS and LAPACK only perform dense matrix operations. – Brian Borchers Aug 23 '21 at 04:18
  • The conventional approach for sparse matrix-matrix multiplication is: don't do it. :p Reframe all your operations as matrix-vector multiplications using associativity, or Schur complements of larger matrices. – Federico Poloni Aug 23 '21 at 07:38
  • @FedericoPoloni So for a matrix-matrix multiplication, I slice one matrix to vectors, make each one dense and then multiply it? – Uwe.Schneider Feb 09 '24 at 20:59
  • 1
    @Uwe.Schneider There are various ways if you really have to do it, but my advice was to use associativity or block matrices to avoid matrix-matrix multiplications if possible. – Federico Poloni Feb 09 '24 at 22:31

1 Answers1

4

Depending on the application, there are many ways to store sparse matrices. The most common ones are CSC, CSR and COO, but there are over 20 types I know of, and probably there are more than those.

Usually, you don't want to multiply sparse matrices together, because that is expensive. If you have to, it is easy to multiply a CSR matrix with a CSC matrix and vice versa. Multiplying two COO matrices together is not that hard, either. Tim Davis' SuiteSparse has some functions capable of doing these operations in parallel.

Multiplying a sparse matrix with a dense matrix is easy. The difficult part is cache management. If you don't care about it, you can take SpMat-Vec operation from Yousef Saad's Sparskit2 and rewrite it for dense matrices. Generalizes easily.

BLAS and LAPACK are not capable of sparse matrix operations. You need to use sparseBLAS/graphBLAS. I don't know any feature complete sparse LAPACK-like library, because many of those operations are pretty difficult to implement for sparse matrices. Also, we know how these libraries work because their development is well-documented, and the results are published. If you have access to a university library, you can read all the papers written since 70s on these topics (from BLAS/LAPACK to thread-parallel and distributed version and sparse/graph linear algebra).

Abdullah Ali Sivas
  • 2,636
  • 1
  • 6
  • 20
  • Could you expand on why multiplying sprase matrices together is expensive? – anonuser01 Aug 23 '21 at 02:44
  • 3
    There are many technical details, like hundreds -if not thousands- of papers full of issues related to SpM-SpM. But the simplest one is related to how sparse matrices are stored. You just store the nonzeros, so when you try, say, $C_{ik} = \sum_j A_{ij}B_{jk}$ you have to go and "search" for the entries $A_{ij}$ and $B_{jk}$. That is the first issue, say you solved it. Now you have a new problem, allocate new memory to store each nonzero $C_{ik}$. Because ideally, you want to store $C$ as a sparse matrix too. But you don't know which entries are nonzero before you do the multiplication. – Abdullah Ali Sivas Aug 23 '21 at 02:57
  • @AbdullahAliSivas Thank you for your answer! Could you tell me where the CSC - CSR algorithm is explained in greater detail? – Uwe.Schneider Feb 09 '24 at 20:59
  • 1
    Maybe you ought to start here https://www.gnu.org/software/gsl/doc/html/spmatrix.html. – A rural reader Feb 10 '24 at 00:49
  • @Aruralreader Thank you for the link! I saw, that the possible operations are described, but not really, how it is actually implemented. – Uwe.Schneider Feb 10 '24 at 10:34