# How to use Compile for accelerating matrix multiplications?

I am trying to find the inverse of the matrix $A$

 m = n = 500;
SeedRandom;
A = RandomReal[10, {m, n}];


iteratively using the following matrix iterative method

$$X_{k+1}=X_k(2I-AX_k)$$

when $X_0=\frac{2}{\sigma_1^2+\sigma_r^2}A^*$, with $\sigma_1$ and $\sigma_r$ as the largest and the smallest singular values of $A$, and $I$ is the identity matrix. This iterative methods can simply be coded in what follows:

 Schulz[X_] := With[{XX = A.X,Id = SparseArray[{{i_, i_} -> 1.}, {m, m}]}, X.(2 Id - XX)];


The whole of this iteration to reach the prescribed tolerance can then be given by the following two-argument function

inverse[A_, tolerance_] := Module[{smax = SingularValueList[A, 1][],
smin = SingularValueList[A, -1][]},
X0 = (2./(smax^2 + smin^2))*ConjugateTranspose[A];
Schulz[X_] := With[{XX = A.X,
Id = SparseArray[{{i_, i_} -> 1.}, {m, m}]}, X.(2 Id - XX)];
FixedPoint[(Schulz[#] &), X0,
SameTest -> (Norm[#1 - #2, Infinity] <= tolerance &)]];


This function could converge to the inverse $A^{-1}$. For example,

tolerance = 10^-6;
B = inverse[A, tolerance]; // AbsoluteTiming
B.A // Chop // MatrixPlot


My main question is here: the running time of finding the inverse using the above piece of code is too much in contrast to the built-in function Inverse[A]. In fact, the above implementation takes 6.2 seconds in my machine while Inverse[A] takes around 0.2 seconds. So, is there anyway to accelerate the above code by using Compile on the matrix iterations?

I think the most time consuming parts are the two matrix-matrix multiplications per step, but is there anyway to use RunTimeAttribute or Parallelization -> True, to become it faster than it is?

I will be thankful if anyone could revise the above implementation.

• I recommend to use your GPU resources, e.g. with CUDALink you can use e.g. CUDADot, for a benchmark take a look at this answer : stackoverflow.com/questions/8638905/… – Artes Mar 11 '13 at 19:49
• Even if you could make the iterations faster, how will you improve speed for finding that smallest singular value? Any power iteration based method will need to solve equations of the form A.x==b and that involves at least a third or so of the cost of matrix inversion. In any case I doubt you'll improve on Dot since that goes to BLAS library code for machine real matrices. – Daniel Lichtblau Mar 11 '13 at 19:53
• As far as I know, for smin and smax, Mathematica 8 uses the Arnoldi algorithm with at most 1000 MaxIterations. However, this part is quite fast. Do you think, is there any way to use Compile for speeding this code up? – Fazlollah Mar 11 '13 at 20:26
• "Do you think, is there any way to use Compile for speeding this code up?" -since it's already using ARPACK in that case you speak of, compilation isn't very useful to do. – J. M.'s discontentment Apr 25 '13 at 17:29
• So, if you think there is any other way to speed up the process, I am happy to be informed. Furthermore, I cannot use CUDADot in my MMA. It fails to be fully downloaded and installed after 30 minutes. Any tips to accelerate the process is fully appreciated. – Fazlollah Apr 26 '13 at 17:48

The answer to your question is 'No'. Compile will not make this code any faster.

If you are running on Windows, launch the Task Manager and watch what happens as you run your code. On my machine, I see that I get 100% CPU Utilisation across all 4 of my 4 processor cores. This is because your program spends the majority of its time running Matrix operations such as Dot which, in turn, call highly optimised, parallelised BLAS functions.

I believe that Mathematica uses the Intel MKL as its BLAS and you won't get much faster than that on Intel CPUs.

If you reach a point where the vast majority of your code's time is spent calling BLAS routines, it is highly unlikley that you'll be able to parallelise it further.

As one comment suggested, you may find a decent speed-up if you switch to using CUDA based routines such as CUDADot if you have the required hardware.

• So, is there any other way such as parallelization? – Fazlollah Mar 11 '13 at 20:55
• the Intel MKL BLAS functions are parallelized. – s0rce Apr 25 '13 at 17:49
• Please write a code solution based on your tip. I am not familiar with your suggestion. – Fazlollah Apr 26 '13 at 17:50