I need to compute a matrix-matrix product, $A^T B$, where $A$ is $n \times r$ sparse, and $B$ is $n \times q$ dense. The number of rows $n$ is far larger than both $r$ and $q$. In fact $n$ is so large I cannot store the entire matrix $B$ in memory, so I build it one row at a time and update the matrix product $A^T B$ accordingly. The way I currently do this is to store $A$ in CSR format and then apply the following algorithm:

loop i = 1, ..., n     // loop over rows of B
  bi = build_row(B, i) // construct row i of B
  loop p = A->row_ptr[i], ..., A->row_ptr[i+1] - 1 // loop over non-zero entries in row i of A
    DAXPY(A->data[p], bi, ATB(A->col_ind[p], :))

This algorithm seems to have complexity $O(nkq)$, where $k$ is the number of non-zero entries in a row of $A$ (this is the same for all rows of $A$). This algorithm works fine, but is quite slow so I am looking for ideas on how to speed it up.

Dense BLAS algorithms often take advantage of blocks of data instead of operating on one row at a time, so my thinking was to do the same (instead of building B one row at a time, build a block of rows, where the block size is still small enough to fit in memory, and then apply a matrix multiplication algorithm). For this type of method, I don't think CSR format would work, since we would need the columns of $A$ in sequential memory locations, so CSC format might work better. A possible algorithm (with A in CSC format) is:

loop i = 1, ..., num_blocks
  row0 = (i-1) * block_size           // row of start of block in B
  row1 = row0 + block_size            // row of end of block in B
  blockB = build_block(B, row0, row1) // build B(row0:row1,:)
  loop j = 1, ..., r // loop over rows of A^T
    loop p = A->col_ptr[j], ..., A->col_ptr[j+1]-1 // loop over non-zero entries in row j of A^T
      if (row0 <= A->row_ind[p]) AND (A->row_ind[p] <= row1) // if this row of A is in the right range, perform DAXPY operation
        DAXPY(A->data[p], blockB(A->row_ind[p]-row0, :), ATB(j, :))

I think the above algorithm has complexity $O(k r q \times NB)$, where $k$ is the number of non-zero entries in a column of $A$, and $NB$ is the number of blocks (num_blocks).

I haven't yet coded up this method, since I'm not sure it will actually be faster. Does anyone have any suggestions on how to make a fast algorithm to do this? Perhaps there is another sparse matrix storage format which would be better to use?

  • $\begingroup$ Have you tried calculating column wise? I would store A in CSC format and use a CSR matrix-vector product (which calculates A^T) on columns of B. $\endgroup$ – slek120 Nov 17 '18 at 16:29
  • $\begingroup$ Unfortunately, for your case, I think it is difficult to predict which format would be better without trying it out. You could maybe try get a smaller block of the matrix and run some benchmarks on that. Additionally, I think it also depends on what architecture you are trying to perform this operation on. But block operations is definitely the way to go, in any case. $\endgroup$ – Mathnoob Nov 18 '18 at 17:32

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.