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.

  • 2
    $\begingroup$ 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. $\endgroup$
    – ssch
    Jan 15 '13 at 13:44
  • 1
    $\begingroup$ @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. $\endgroup$
    – Yves Klett
    Jan 15 '13 at 14:02
  • 2
    $\begingroup$ 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 LinearAlgebra`BLAS`GEMV. $\endgroup$ Jan 15 '13 at 14:34
  • 3
    $\begingroup$ 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. $\endgroup$
    – rcollyer
    Jan 15 '13 at 15:30
  • 2
    $\begingroup$ Related question: Efficiently Constructing Rank One Approximations for a Matrix using SVD $\endgroup$
    – Jens
    Jan 15 '13 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]]




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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.