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 agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.