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, 2013 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. Jan 15, 2013 at 14:02
• 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. Jan 15, 2013 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. Jan 15, 2013 at 15:30
• – Jens
Jan 15, 2013 at 17:56

SVD is an $$O(n^3)$$ operation (and truncated rank-$$r$$ factorization is probably still of complexity $$O(r \, n^2)$$, so pretty slow. If your matrix A has indeed low rank $$r \ll n$$, then Adaptive Cross Approximation (ACA) can provide a rank-$$r$$ factorization in $$O(r^2 \, n)$$ time and $$O(r \, n)$$ additional memory. This algorithm is greedy, so you can prescibe a given relative accucary and it will find also the required rank $$r$$. The resulting factors usually have a size that is only insignificantly larger than the factors provided by truncated SVD. The method is also very powerful because it needs only a few columns and rows of the matrix $$A$$ so that $$A$$ does not need to be fully assembled in memory. If course, it is not 100% reliable: Think of a sparse matrix with few nonzero rows and columns; these will quite certainly not be detected. But it works greate provided that the matrix $$A$$ is sufficiently smooth, (i.e., if it originates from sampling a smooth integral kernel $$K(x,y)$$),

Here is an example with a 400-fold speedup. It features a squared distance matrix A which really has (numerically) low rank ($$r \leq 6$$). Please find the code for ACACompression here.

n = 5000;
A = DistanceMatrix[RandomReal[{-1, 1}, {n, 3}]]^2;
u = RandomReal[{-1, 1}, {n}];

{U, V} = ACACompression[A, "MaxRank" -> 150, "Tolerance" -> 1. 10^-12]
U = Transpose[U];

u = RandomReal[{-1, 1}, {n}];
btrue = A.u; // RepeatedTiming // First
bproxy = U.(V.u); // RepeatedTiming // First

Max[Abs[btrue - bproxy]]


0.0080

0.000019

6.25278*10^-13

The factorization took 0.0042 seconds. For comparision: SingularValueDecomposition[A] took 23.6235 seconds on my machine.