More efficient matrix-vector product

Question

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 \begin{equation} A = \sigma_1 u_1 v_1^\text{T} + \sigma_2 u_2 v_2^\text{T} + \ldots \end{equation} so \begin{equation} Ax = (\sigma_1 v_1^\text{T}x) u_1 + (\sigma_2 v_2^\text{T})u_2 + \ldots \end{equation} 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.

Answer 1

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.