CPU timing for matrix rank calculation difference between 10.3 and 10.4 and 11.0

Question

I run a small rank CPU test when new Mathematica version comes out. I noticed when version 10.1 came out that Mathematica rank calculation suddenly became extremely fast on a matrix when the size became over 2500 by 2500. (link at the very end). Here is from the link:

Mathematica 10.1 was surprisingly much faster on this test than 10.0.2. It seems Mathematica 10.1 is using different algorithm to compute the rank now to account for this drastic difference in speed improvement.

The speed boost was observed to occur at certain matrix size. At matrix size of 2500 or less, the same speed was obtained as with version 10.0.2. At matrix size over 2500, even by just one, a dramatic speed increase was seen. For n = 2500 Mathematica CPU was around 4.6 seconds which is the same as in 10.0.2, but by increasing the matrix size to n = 2501 , CPU time went down to about 1.4 seconds. This is 3 times as fast for essentially the same matrix size. This result was reproducible. This seems to indicate that Mathematica internally uses the same algorithm as previous version for smaller size matrices, and then switches to different algorithm for larger matrices.

The above was the case on windows up to 10.3.

Suddenly, in 10.4 and 11, this sudden speed up in rank calculation went away, and the CPU time went back to what it was before version 10.

Would some one be able to explain what happened to the speed improvement? and why now Mathematica is slower on this test than on 10.3? This is all on same PC, all on 64 bit windows 7.

This is the source code for the test

Remove["Global`*"];
$HistoryLength = 0;
Share[];
n = 2500;
m = RandomReal[{}, {n, n}];
AbsoluteTiming[MatrixRank[m];]

Here is the result. Notice the speed up when the size becomes 2501

version 10.3

version 11

Similar result for 10.4 as above. So something changed between 10.3 and 10.4 which carried over to 11.

The full detailed report can be found here

---

Added data:

size = Range[500, 8000, 500];
ver10 = {0.03, 0.14, 0.65, 2.15, 4.66, 2.25, 3.42, 4.98, 6.97, 8.8, 
   9.66, 12.81, 14.7, 17.82, 22.49, 27.04};
ver11 = {0.024, 0.134, 0.616, 2.033, 4.505, 8.393, 13.865, 22.088, 
   30.846, 43.079, 58.216, 75.655, 96.048, 120.505, 148.593, 180.311};
data = Transpose[{size, ver10, ver11}]
Grid[Join[{{"N", "version 10.1 CPU", "version 11 CPU"}}, data], 
 Frame -> All]

p1 = ListPlot[Transpose@{size, ver10}, Mesh -> All, Joined -> True, 
   PlotStyle -> Green, Frame -> True, 
   FrameLabel -> {{"CPU time (sec)", None}, {"Matrix size", 
      "Comparing CPU time for rank calculation"}}];
p2 = ListPlot[Transpose@{size, ver11}, Mesh -> All, Joined -> True, 
   PlotStyle -> Red];
Show[p1, p2, PlotRange -> All]

UPDATE

---

I think I have found something very important and I have a theory about what has happened. I run the rank test now on Matlab 2016a and also on Maple 2016, and discovered that the speed boost now is present in both Maple and Matlab. This tells me that the speed boost in rank calculation came from external library, which would be intel MKL math library which all three M's link to. So the speed boost happend before in Mathematica, simply because 10.3 was shipped linked to that library.

Now the newer versions of Maple and Matlab are also linked to it. This is why the speed boost did not show up in the earlier versions of Maple and Matlab I used when I run the test before.

The question then, why did the speed boost go away in 10.4 and 11.0? The only thing that could explain this, is that Mathematica is not linked to the current intel MKL math library which had this speed in rank calculations? What else could explain this?

Here is Maple result and Matlab result with the new versions. Notice the sudden speed in rank calculation now, the same as used to be with 10.3 and at the same exact size: All on same PC, all 64 bit.

Maple 2016.1:

Matlab 2016a:

>> clear all; n=2500; A=rand(n,n); tic();rank(A);toc()
Elapsed time is 4.549381 seconds.

>> clear all; n=2501; A=rand(n,n); tic();rank(A);toc()
Elapsed time is 1.436199 seconds.
>> 

Question is: Is Mathematica 11 linked to the correct/latest intel MKL library on windows? What does the test show on mac with version 11? I need now to try to find out what intel MKL library Mathematica 11 in linked to and compare that to 10.3

Update

---

Here is more information. It seems the MKL version 11.2 had the speed boost in it based on checking the MKL DLL file version on the CAS I have. And that MKL 11.3 do not have this speed boost for some reason. This is assuming the issue is to due to MKL

From the above one can see that when MKL 11.2 is used, the three M's had fast rank calculations for Matrix size over 2500. Mathematica 10.4 and 11 changed from MKL 11.2 to 11.3 and I am guessing this is the reason why rank got slower. This is only a theory ofcourse.

Answer 1

Edit 2019-09-28: OK, now I got a full picture:

It is the changes of SVD algorithm in MKL that affect speed of MatrixRank[] in Mathematica.

Particularly, the SVD routine dgesdd got speed boost only in MKL version 11.2, hence also the Mathematica.

A picture (Hum...a table) worth a thousand words:

This table shows the implementations of the functions in left most column, in different time period marked by MKL versions.

The red dgerdb is the core of fast SVD computation. These cipher-like dge??? names are LAPACK functions, see original post below to understand them.

Clearly, MKL v11.2 speed boost all its downstream softwares (the "3M"). But after that, only who use dgesvd got benefits. Mathematica always use dgesdd for SVD. But for an unclear reason, Matlab switch to the slow dgesdd at some point. I did a preliminary test for the fast code, and found no precision/accuracy lost compared to the slow one.

This table is obtained through techniques mention below, particularly the LD_PRELOAD trick. Here is some code I used. And thanks the university HPC for Matlab/Mathematica.

∎

---

(Original post till the end)

It is probably attribute to the change of SVD algorithm in MKL.

I don't have multiple versions of Mathematica or Matlab, and I run only Linux. But some implementation details are dig out:

• Mathematica and Matlab use MKL for rank computation.
• The rank is computed through SVD.
• In Mathematica8 the LAPACK routine DGESDD is called for SVD.
• In Matlab2017b the LAPACK routine DGESVD is called for SVD.
• In reference LAPACK implementation version 3.8, both DGESDD and DGESVD eventually call DBDSQR (implicit zero-shift QR algorithm) for actual singular value computation, as long as no singular vector is required. This explains why Mathematica and Matlab behavior the same even using different SVD routines.

MKL updated SVD code to "improve" the speed during that period:

What's New in Intel MKL 11.2

Improved performance of (S/D)GE(SVD/SDD) when M>=N and singular vectors are not needed

What's New in Intel MKL 11.3

Improved performance of SVD for cases where singular vectors are computed on Intel AVX or Intel AVX2 architectures, and when M>=N and singular vectors are not needed.

The DLL or SO file of corresponding MKL versions might be needed to further find out what happen to them.

Details

Here is a detailed steps to find out how Mathematica computes MatrixRank.

• Start Mathematica, and run the MatrixRank when necessary.

• top to find the Process ID that consumes CPU during computation.

  In my case, the PID is 8820 named MathKernel.

• perf top -p 8820 to find the time consuming function and library. 1

  The first line (most time consuming) is "34.34% libmkl_avx.so [.] LA16Y8_Loop_M16gas_1". So we have a clear evidence that Mathematica is using MKL. But the function "LA16Y8_Loop_M16gas_1" is not documented, so at this point it is not so clear that how the rank is computed.

• ltrace -T -p 8820 [2]. Here we trace how Mathematica calls its library functions.

  Among the recorded calls, we can find the related function dgesdd. The dgesdd consumes 6.486 sec which occupies majority of the computation time (8.850 sec).

  "dgesdd" is the name of standard LAPACK subroutine to compute singular value decomposition (SVD), see LAPACK dgesdd doc.

  Also, look around the trace records, there is a call pattern dgesdd, malloc then dgesdd which is a typical usage of dgesdd, the first two calls are for workspace query and memory allocation.

  As we know, SVD is the standard means to get matrix rank. Now we have a very strong evidence that Mathematica is using SVD (particully dgesdd) to compute matrix rank.

For Matlab, similar steps above also works, except that ltrace is failed to discover all calls to DGESVD. To fix that, the LD_PRELOAD trick is used to further confirm that DGESVD is the workhorse.

1 Linux perf tool, specially the top command.

[2] The actual command I used is ltrace -T -o ltr.txt -e -MLUserData-pthread_mutex_lock-gettimeofday-pthread_mutex_unlock-pthread_cond_timedwait -p 8820. For matlab it is ltrace -o ltrM.txt --demangle -e '@*lapack*' -e '@*mkl*' -e '@*blas*' -f -T /usr/local/MATLAB/R2017b/bin/glnxa64/MATLAB -nojvm.

Discussion

Now the question is: why DGESVD or DGESDD are changing speed, assuming Mathematica/Matlab did not change its MatrixRank implementation which is suggested by the version table in the asked question, then the only possibility is: MKL changed. MKL's release note does suggest this hypothesis (see the quotes in beginning), but we still don't know what happen to DGESVD/DGESDD in terms of algorithm or implementation.

The (reference) LAPACK project is sponsored in part by MathWorks and Intel, the reference implementation and MKL implementation might shared something in common. But after looking into (reference) LAPACK update history, I find nothing but bug fix about SVD.

Edit 2019-09-26: The below part (now quoted in yellow) is proven wrong, see the update after it instead.

I would like to add a wild guess here, in MKL 11.2, probably someone decided to use the "divide and conquer method" (wiki) for SVD even when only singular value is needed. See also the suggestion in comment of DBDSDC (which will be called by DGESDD).

The code currently calls DLASDQ if singular values only are desired. However, it can be slightly modified to compute singular values using the divide and conquer method.

Then one day, the MKL team found that divide-and-conquer is not accurate in many cases, so went back to the slightly slower but more accurate QR algorithm.

Edit 2019-09-26 continued.

I actually have the "fast" version MKL, so investigate further here.

The new findings are:

• Matlab2017b use MKL 11.3.1, but it is "fast".

  Computing rank of matrix size n=2501 (0.989 sec) is much fast than n=2500 (4.98 sec).

• In MKL 11.3.1, the dgesvd behave differently for n=2500 and n=2501.

  For n=2500, dgesvd will call mkl_lapack_dgebrd, dgebrd is a standard LAPACK routine to reduces a general real matrix A to upper or lower bidiagonal form B by an orthogonal transformation Q**T * A * P = B.

  For n=2501, dgesvd will call mkl_lapack_dgerdb for the same purpose (presumably). mkl_lapack_dgerdb and dgerdb are undocumented, but dsyrdb is documented, it can be found in MKL Reference. With help of that manual and following LAPACK naming convention, dgerdb means double type - general matrix - successive bandwidth reduction.

  The bidiagonalization here is a standard way to compute SVD. The LAPACK routine dbdsqr is specifically for SVD of bidiagonal matrix, and dbdsqr is been called during SVD computation in my Matlab/Mathematica. This helps us to further confirm the meaning and function of dgerdb.

  The tool SystemTap is used to discover the detail of how these functions are called.

• Bidiagonalization is the most time consuming part in SVD.

  In my case, for n=2500 bidiagonalization cost 4.845 sec, 97% of SVD time (4.973 sec), or 87% for n=2501.

  Hence my previous theory is wrong, 'divide-and-conquer' cannot help much since it need bidiagonalization first.

• Mathematica is slow for MKL 11.3, because it is using dgesdd.

  dgesdd always call dgebrd for bidiagonalization, thus no speed boost for n>=2501. This can be confirmed by replace Matlab call from dgesvd to dgesdd (again using the LD_PRELOAD trick), which make Matlab slows down like Mathematica.

• dgerdb (the fast one) do not exist in old MKL (e.g. Mathematica 8).

  Confirmed by look into the symbol table of MKL library of Mathematica 8.

As a summary of this edit, MKL is observed improving SVD algorithm (dgesvd), but Mathematica probably missed these benefits by using a different SVD routine (dgesdd).