Speed up sparse Hermitian matrix-vector product

Question

I am interested in speeding up as much as possible the dot product in the specific case of a multiplication of a $n\times n$ sparse hermitian matrix $H$ (where $n$ is typically large, let's say $\approx 10^4$) and a vector $v$ with $n$ complex entries. Obviously, H.v does the job using the function Dot[], but I am not sure that this is the fastest way to go (or is it?). For example, can I create a compiled function with the list of indexes and values of $H$ and $v$ as inputs and do better?

EDIT As suggested in the comments, here's a code to generate randomly defined $H$ and $v$:

n = 10000;(* dimension of matrix H *)
m = 50000;(* number of non-vanishing entries *)
pos = RandomInteger[{1, n}, {m, 
    2}];(* list of positions of non-vanishing entries *)
vals = RandomComplex[{-1.0 - 1.0*I, 1.0 + 1.0*I}, 
   m];(* list of values *)
H = SparseArray[
   Thread[pos -> vals]
   , {n, n}];(* make a non-hermitian sparse array *)
H += ConjugateTranspose[H];(* make it hermitian *)
HermitianMatrixQ[H] (* check if it's hermitian *)
v = RandomComplex[{-1.0 - 1.0*I, 1.0 + 1.0*I}, n];
H . v // 
  RepeatedTiming // First (* time taken by the function Dot[] (4 ms \
on my machine) *)

Answer 1

If you happen have to multiply many vectors at once with the same matrix H, then place them as columns into a dense matrix. Afterwards do a single maxrix-matrix multiply:

v = RandomComplex[{-1.0 - 1.0*I, 1.0 + 1.0*I}, n];
w = H . v // RepeatedTiming // First

V = RandomComplex[{-1.0 - 1.0*I, 1.0 + 1.0*I}, {n, 8}];
W = H . V // RepeatedTiming // First

0.00178638

0.00404287

If you would do 16 separate matrix-vector multiplication 0.001761228515625 seconds (which would 0.0281797), it takes only 0.00404287. So this is a speedup of factor 7. This has to do with the scattered read-operations that have to be performed in the CSR format. You face many cache misses. But the number of cache misses are roughly the same for the a single matrix-vector multiply and the matrix-matrix multiply. This is why the the CPU typically cannot work with full throughput: It has to way for the data to arrive. But in the way fetching is done, much of the row of the matrix V is already loaded when the first entry is accessed. This helps the to keep the CPU busy while the next row can be loaded for the multiplication (which can be any row below the current row because of the sparsity of H.

Hence the precise speedup depends very much on the machine you run it on (CPU speed, number of floating point or vector pipelines, cache sizes, memory bandwidth...)

Answer 2

Wow, I just stumbled over something that I cannot explain. Please try to create the matrix H in the following way (note the explicit mention of the background value:

H = SparseArray[pos -> vals, {n, n}, 0. + 0. I];
H = SparseArray[H + ConjugateTranspose[H]];
H . v; // RepeatedTiming // First

0.000205734

This is almost 8 times faster than the original timing (0.00176442) on my machine. Also applying SparseArray after the addition is a good practice to guarantee accordance to the CSR format (correct ordering of the column indices). But apparently, only specifying the background value as complex number made Mathematica pick the most efficient code path...