4
$\begingroup$

I want to compute the matrix $$ A = \sum_{i=1}^N v_i v_i^T $$ where each $v_i$ is a given vector of length $2500$, so that $A$ is $2500 \times 2500$, and my $N$ is about 2 million. Rather than call DSYR 2 million times, I found I can speed up the calculation significantly, by storing the $v_i$ in block matrices, and then calling DSYRK when the block matrix is full. I define a $\textrm{block_size} \times 2500$ matrix $W_k$ as: $$ W_k = \pmatrix{ v_1^T \\ \vdots \\ v_b^T } $$ so that $$ A = \sum_{k=1}^{\textrm{nblocks}} W_k^T W_k $$ The reason I set the rows of $W_k$ equal to the vectors $v_i$ and not the columns is because I use C and this way the incX of the rows will be 1 which I think helps to speed up the BLAS operations.

I am finding this procedure is still slower than I would like, and my question is how do I choose the $\textrm{block_size}$ (i.e. number of rows of $W_k$) to get the best performance?

Currently, I am choosing $\textrm{block_size} = 50000$, which I only chose so that the total size of the $W_k$ storage matrix is 1GB ($50000 \times 2500 \times 8 = 1 GB$). Is there a better choice of $\textrm{block_size}$ (number of rows of $W_k$) to get optimal performance?

In case it matters I am using the ATLAS CBLAS library on a machine with 24 cores.

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ You’re going to have to access the matrix both by rows and columns no matter whether your start with your $W$ matrix or its transpose. It’s possible that your BLAS might be more efficient one way or the other. You could also try other BLAS implementations (MKL, OpenBLAS, etc.) $\endgroup$
    – Brian Borchers
    Commented Jan 11, 2019 at 19:38
  • $\begingroup$ Of these ~2 million vectors $v_i$, at most 2500 can be linearly independent. Can't you use that somehow? $\endgroup$
    – Christoph
    Commented Jan 14, 2019 at 7:24

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.