10

I am new to Julia and primarily work in Mathematica, so I probably have a few elementary mistakes floating around. I attempted to time how long Julia took to compute the eigensystem of a random matrix, and found it was 5-6 times slower than in Mathematica.

In Julia:

D=1000*(rand(1000,1000)-0.5);
@time (E,F)=eig(D);

Out: elapsed time: 7.47950706 seconds (79638920 bytes allocated*)

In Mathematica:

First@Timing@Eigensystem[RandomReal[{-500, 500}, {1000, 1000}]]

Out: 1.310408

For 2000 x 2000 arrays it's similar, although the Julia result slowed down slightly less than the equivalent Mathematica call, but it's still slower; Julia takes 22 seconds, whereas Mathematica computes it in 8 seconds.

As far as I read in the Julia standard library for linear algebra, decompositions are implemented by calling LAPACK, which I thought was supposed to be very good, so I'm confused as to why the Julia code is running so much slower. Does anyone know why this is the case? Is it doing some kind of balancing or array-symmetry-detection that Mathematica doesn't do? Or is it actually slower?

Also, this is a syntax question and probably a silly error, but how do you change the balancing in Julia? I tried

@time (E,F)=eig(D[, balance=:nobalance]);

exactly as copied and pasted from the Julia manual, but it just gave a syntax error, so something's wrong.

I am using Windows 7 64-bit, with Julia version 0.2.0 64-bit, installed using the instructions at Steven Johnson's site, with Anaconda installed first to take care of prerequisites. I am using Mathematica student edition version 9.0.1.

EDIT 1:

Executing versioninfo() yielded

Julia Version 0.2.0
Commit 05c6461 (2013-11-16 23:44 UTC)
Platform Info:
System: Windows (x86_64-w64-mingw32)
WORD_SIZE: 64
BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY)
LAPACK: libopenblas
LIBM: libopenlibm

So it looks like I'm using the openBLAS for LAPACK and BLAS. Once I get the Mathematica implementation info I will add that as well.

EDIT 2:

It appears that Windows Mathematica probably uses Intel MKL BLAS.

Community
  • 1
  • 1
DumpsterDoofus
  • 1,132
  • 12
  • 30
  • 1
    Do you know what [BLAS](http://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms#Implementations) implementation your Julia is using? There can be big performance differences if it isn't one that is tuned for your architecture. – rhashimoto Feb 08 '14 at 05:48
  • One more evidence that the benchmarks on Julia website aren't fair. I experienced similar performance issue [here](http://stackoverflow.com/questions/19208014/is-the-julia-language-really-as-fast-as-it-claims). – juliohm Feb 08 '14 at 10:39
  • 1
    @julohm, the results are real and express genuine performance differences between languages. But the benchmarks don't test everything, and even in a fast language there can always be bugs or non-optimal algorithms. Instead of being dismissive, next time try filing an issue; you may be surprised at what's possible. – tholy Feb 08 '14 at 13:15
  • @rhashimoto: Based on the answers below I think that could be the issue. Is there a way to check what implementation of LAPACK/BLAS Julia is using? I searched my hard drive for "BLAS", and found a bunch of hits in C:\Users\Username\Anaconda\Lib\site-packages\scipy\lib (I installed Julia according to the MIT instructions by first installing Anaconda to get the prerequisites), and in C:\Users\Username\Downloads\julia-05c6461b55\share\julia\test\perf, which contained what appeared to be some performance tests, but no obvious mention of what implementation was used. – DumpsterDoofus Feb 08 '14 at 14:48
  • If you run `versioninfo()` at the julia prompt, it should show a bunch of information including what BLAS, LAPACK and LIBM version you're using. This kind of thing makes me wonder if we shouldn't just always use OpenBLAS unless explicitly asked to use something else. – StefanKarpinski Feb 08 '14 at 16:40
  • @StefanKarpinski: Thanks, I added the version info to my question. Once I get the info about Mathematica, I'll add that too. – DumpsterDoofus Feb 08 '14 at 17:47
  • @StefanKarpinski: So according to one member at MathematicaSE, Windows versions of Mathematica use the Intel Math Kernel Library BLAS. – DumpsterDoofus Feb 08 '14 at 19:31

2 Answers2

14

The eigen calculation in Julia is outsourced to LAPACK and BLAS, and I think it is also the case for Mathematica. Julia can use different versions of BLAS and LAPACK and you are therefore effectively comparing your choice of LAPACK and BLAS for Julia with Mathematica's LAPACK and BLAS (probably Intel MKL).

The default choice for Julia is OpenBLAS which is fast on most architectures and on my machine Julia is faster than Mathematica for the eigen calculation. If you are on Linux and have chosen BLAS and LAPACK from a repo, it is very likely that they are much slower than OpenBLAS.

The option for balancing has recently been added to Julia and mistakenly the option was not added to the function eig, which is only a MATLAB compatible interface to the eigfact function. Writing eigfact(A,balance=:nobalance) should work.

Edit 1: Further investigation has shown that the difference is due to a threading problem in OpenBLAS on Windows. If Julia's BLAS is restricted to one thread the timings are comparable to Mathematica, but if more threads are allowed the calculation slows down. This doesn't seem to be a problem on Mac or Linux, but as mentioned above, in general the performance of OpenBLAS depends on the architecture.

Edit 2: Recently, the balancing option has changed. Balancing can be switched off by writing eigfact(A,permute=false,scale=false).

Andreas Noack
  • 1,370
  • 6
  • 11
  • 1
    Interesting, I didn't know Mathematica outsources its dense array computations to BLAS/LAPACK, as the documentation for Eigensystem and similar operations doesn't mention this, although in the documentation's tutorial/SomeNotesOnInternalImplementation, it mentions "For dense arrays, LAPACK algorithms extended for arbitrary precision are used when appropriate" and "BLAS technology is used to optimize for particular machine architectures", but nothing more. I'll separately ask MathematicaSE how to find out what BLAS implementation I'm using, but do you know how to find out what Julia is using? – DumpsterDoofus Feb 08 '14 at 14:54
  • 1
    Software that solves eigen problems have many names, but behind the scene it is almost always LAPACK and BLAS doing the actual work. You can get some info about julia by running the function `versioninfo()`. If you are using OpenBLAS it should return something like `BLAS: libopenblas`. The slow reference BLAS could have the name `libgfortblas`. – Andreas Noack Feb 08 '14 at 18:07
6

There's definitely at least one thing wrong, and of couse the good news is that it's likely to be fixable. For this kind of thing, you're much better off filing an issue on GitHub. I've done that for you here.

Regarding the speed, you may want to check that the accuracy of the eigendecomposition is comparable. And of course, accuracy sometimes depends on the condition number of the matrix; it's possible that Julia is using a more careful algorithm, and your example may or may not reveal that if its condition number is not a problem. This is definitely something to discuss over at the issue.

Regarding :nobalance, in the documentation, [, something] is used to indicate that something is optional. You want to use @time (E,F)=eig(D, balance=:nobalance);. However, eig doesn't take keyword arguments, so the code and the documentation are not currently in sync.

tholy
  • 11,882
  • 1
  • 29
  • 42
  • Thank's for the GitHub thing, I'll check it out. I don't know anything about language internals, but does the fact that `eig` doesn't take keywords while `eigfact` does mean that `eig` is not simply a renamed version of `eigfact`? I guess it's just a minor documentation thing, so it's not too important. – DumpsterDoofus Feb 08 '14 at 15:00