How to use Compile for accelerating matrix multiplications?

Question

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

 m = n = 500;
 SeedRandom[123456];
 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][[1]],
    smin = SingularValueList[A, -1][[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.

Answer 1

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.