# Fast calculation of $A^T B$

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?

• 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. – slek120 Nov 17 '18 at 16:29
• 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. – Mathnoob Nov 18 '18 at 17:32