# More efficient matrix-vector product

Dear mathematica users,

In my present research I am faced with a real dense $n\times n$ matrix $A$ where $n \geq 3000$ (hopefully even more). The coefficients of this matrix are fixed, but I will have to repeatedly multiply it by a variable vector: $Ax$.

I am not complaining about Mathematica's speed to do the job, which seems quite nice, but since I will need to do this repeatedly for very many times, I was wondering if there was a way to optimise the process. Perhaps, declaring it in a With will help, but besides that I am out of ideas.

Another alternative would be a low-rank approximation using SVD. My thoughts were: with the SVD I can write $$A = \sigma_1 u_1 v_1^\text{T} + \sigma_2 u_2 v_2^\text{T} + \ldots$$ so $$Ax = (\sigma_1 v_1^\text{T}x) u_1 + (\sigma_2 v_2^\text{T})u_2 + \ldots$$ As an example, using a rank 100 approximation to a $3000\times 3000$ matrix (which yields a Frobenius error of $\sim5-10\;\%$) I was able to reduce the computation time by a factor of roughly 3 or 4.

I thank in advance for any positive feedback.

Best regards,

Gabriel Landi

Edit: Forgot to say that $A$ is symmetric and has zero diagonal.

-
Do you have a graphics card that supports CUDA or OpenCL? I've found CUDADot to be very fast even with a pretty weak graphics card. – ssch Jan 15 '13 at 13:44
@ssch the CUDA part of CUDADot is very fast, but depending on the problem the overhead for a call can make it slower than a standard Dot. – Yves Klett Jan 15 '13 at 14:02
@YvesKlett You are right. I did some testing and I was mistaken, only with matrix-matrix multiplication was there a real advantage. Do you know more than one of the x vectors at a time so you can put them in a matrix? – ssch Jan 15 '13 at 14:14
Improving the speed of this product will be extremely difficult since it is already highly optimized. The only way I can think of apart from CUDA is a direct call to LinearAlgebraBLASGEMV. – Oleksandr R. Jan 15 '13 at 14:34
Does A have a full rank? If so, can you work in its eigenspace? Alternatively, can you calculate $A X$ where $X = (\begin{matrix}\vec{x}_1&\vec{x}_2&\cdots\end{matrix})$, instead? That would allow you to use CUDA with full matrix-matrix computations, and the eigenspace idea is still applicable. – rcollyer Jan 15 '13 at 15:30