Eigenkompositionen sind in Julia fünfmal langsamer als in Mathematica?

Question

Ich bin neu bei Julia und arbeite hauptsächlich in Mathematica, also habe ich wahrscheinlich einige elementare Fehler im Umlauf. Ich versuchte herauszufinden, wie lange Julia brauchte, um das Eigensystem einer zufälligen Matrix zu berechnen, und fand es 5-6 mal langsamer als 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 

Für 2000 x 2000 Arrays es ähnlich ist, obwohl die Julia Ergebnis etwas weniger verlangsamt als das Äquivalent Mathematica Ruf, aber es ist immer noch langsamer ; Julia braucht 22 Sekunden, während Mathematica es in 8 Sekunden berechnet.

Soweit ich in der Julia standard library for linear algebra gelesen habe, werden Dekompositionen durch den Aufruf von LAPACK, die ich dachte, dass es sehr gut sein sollte, implementiert, so dass ich verwirrt bin, warum der Julia-Code so viel langsamer läuft. Weiß jemand, warum das so ist? Führt es eine Art Symmetrie- oder Array-Symmetrie-Erkennung durch, die Mathematica nicht tut? Oder ist es eigentlich langsamer?

Auch dies ist eine Syntaxfrage und wahrscheinlich ein dummer Fehler, aber wie ändert man den Ausgleich in Julia? Ich habe versucht,

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

genau wie kopiert und aus der Julia manuell eingefügt, aber es gab nur einen Syntaxfehler, so etwas nicht stimmt.

Ich benutze Windows 7 64-Bit, mit Julia Version 0.2.0 64-Bit, installiert mit der instructions at Steven Johnson's site, mit Anaconda installiert zuerst, um die Voraussetzungen zu kümmern. Ich benutze Mathematica Student Edition Version 9.0.1.

EDIT 1:

Executing versioninfo() ergab

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 sieht es aus wie ich die openBLAS für LAPACK und BLAS bin mit. Sobald ich die Mathematica-Implementierungsinformation bekomme, werde ich das auch hinzufügen.

EDIT 2:

Es scheint, dass Windows Mathematica probably uses Intel MKL BLAS.

Answer 1

Die Eigenberechnung in Julia ist in LAPACK und BLAS ausgelagert, und ich denke, das ist auch bei Mathematica der Fall. Julia kann verschiedene Versionen von BLAS und LAPACK verwenden und vergleicht daher effektiv Ihre Wahl von LAPACK und BLAS für Julia mit Mathematicas LAPACK und BLAS (wahrscheinlich Intel MKL).

Die Standardauswahl für Julia ist OpenBLAS, die auf den meisten Architekturen schnell ist und auf meiner Maschine Julia ist schneller als Mathematica für die Eigenberechnung. Wenn Sie unter Linux sind und BLAS und LAPACK aus einem Repo ausgewählt haben, ist es sehr wahrscheinlich, dass sie viel langsamer als OpenBLAS sind.

Die Option zum Ausgleich wurde vor kurzem zu Julia hinzugefügt, und fälschlicherweise wurde die Option nicht zur Funktion eig hinzugefügt, die nur eine MATLAB-kompatible Schnittstelle zur eigfact-Funktion ist. Schreiben eigfact(A,balance=:nobalance) sollte funktionieren.

Bearbeiten 1: Weitere Untersuchungen haben gezeigt, dass der Unterschied auf ein Threading-Problem in OpenBLAS unter Windows zurückzuführen ist. Wenn Julias BLAS auf einen Thread beschränkt ist, sind die Timings vergleichbar mit Mathematica, aber wenn mehr Threads erlaubt sind, wird die Berechnung langsamer. Dies scheint kein Problem auf Mac oder Linux zu sein, aber wie oben erwähnt, hängt die Leistung von OpenBLAS generell von der Architektur ab.

Bearbeiten 2: In letzter Zeit hat sich die Auswuchtoption geändert. Das Auswuchten kann durch Schreiben von eigfact(A,permute=false,scale=false) ausgeschaltet werden.

Answer 2

Es ist definitiv mindestens eine Sache falsch, und die gute Nachricht ist, dass es wahrscheinlich reparierbar ist. Für solche Dinge ist es viel besser, ein Problem auf GitHub zu stellen. Ich habe das für Sie getan here.

In Bezug auf die Geschwindigkeit möchten Sie vielleicht überprüfen, dass die Genauigkeit der Eigenkomposition vergleichbar ist. Und natürlich hängt die Genauigkeit manchmal von der Zustandsnummer der Matrix ab; Es ist möglich, dass Julia einen sorgfältigeren Algorithmus verwendet, und Ihr Beispiel zeigt dies möglicherweise nicht oder zeigt es nicht an, wenn seine Bedingungsnummer kein Problem ist. Dies ist definitiv etwas, worüber in der Sache diskutiert werden kann.

In Bezug auf :nobalance wird in der Dokumentation [, something] verwendet, um anzuzeigen, dass something optional ist. Sie möchten @time (E,F)=eig(D, balance=:nobalance); verwenden. eig akzeptiert jedoch keine Schlüsselwortargumente, sodass der Code und die Dokumentation derzeit nicht synchronisiert sind.