# Speed up sparse Hermitian matrix-vector product

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[
, {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) *)

• Please add a sample matrix and vector so people have something to play with as they experiment. Commented Mar 12 at 12:21
• "...can I create a compiled function with the list of indexes and values of 𝐻 and 𝑣 as inputs and do better?" – This is precisely what a sparse-sparse multiplication does and I doubt you can do better. However, if you can calculate/enumerate the sparse matrix elements faster than you can retrieve them from a pre-computed list, then you could write a function $\vec{v}\mapsto H\cdot \vec{v}$ that is faster than a sparse-sparse matrix multiplication. This becomes particularly relevant for extremely large vectors, or if you use ARPACK directly. Commented Mar 12 at 12:36
• Do you happen to have to multiply many vectors at once? Commented Mar 12 at 14:05
• My experience with these special array types is that they are not faster than PackedArrays, and I believe that you can use packed arrays in this case. Also, take a look at the code here, maybe you can copy and paste it to your problem, with some change for matrix-vector multiplication. Commented Mar 12 at 21:12
• There is a matrixMul.cu example for CUDALink that works on real dense matrices. However, I tried this and found it slower than just a CPU A.B in Mathematica. Commented Mar 12 at 21:45

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...)

• Wow! This will be interesting in the future, but for now I have to compute $H\cdot v$ with a single vector, repeating the procedure multiple times as the vector changes iteratively Commented Mar 12 at 22:03

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...