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

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

iteratively using the following matrix iterative method


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.

  • 1
    $\begingroup$ 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/… $\endgroup$
    – Artes
    Commented Mar 11, 2013 at 19:49
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Mar 11, 2013 at 19:53
  • $\begingroup$ 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? $\endgroup$
    – Faz
    Commented Mar 11, 2013 at 20:26
  • $\begingroup$ "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. $\endgroup$ Commented Apr 25, 2013 at 17:29
  • $\begingroup$ 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. $\endgroup$
    – Faz
    Commented Apr 26, 2013 at 17:48

1 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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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