2017-07-25 3 views
4

Ich versuche derzeit, riesige Matrizen der Größenordnung 1 Million um 1 Million zu invertieren, und ich dachte mir, dass der Backslash-Operator dabei hilfreich sein wird. Irgendeine Idee, wie es umgesetzt wird ?. Ich habe keine konkreten Beispiele gefunden, daher wird jede Hilfe sehr geschätzt.Wie verwende ich den Backslash-Operator in Julia?

Antwort

8

Irgendeine Idee, wie es umgesetzt wird?

Es ist ein Multialgorithmus. Dies zeigt, wie es zu benutzen:

julia> A = rand(10,10) 
10×10 Array{Float64,2}: 
0.330453 0.294142 0.682869 0.991427 … 0.533443 0.876566 0.157157 
0.666233 0.47974 0.172657 0.427015  0.501511 0.0978822 0.634164 
0.829653 0.38.589555 0.480963  0.606704 0.642441 0.159564 
0.709197 0.570496 0.484826 0.17325  0.699379 0.0281233 0.66744 
0.478663 0.87298 0.488389 0.188844  0.38193 0.641309 0.448757 
0.471705 0.804767 0.420039 0.0528729 … 0.658368 0.911007 0.705696 
0.679734 0.542958 0.22658 0.977581  0.197043 0.717683 0.21933 
0.771544 0.326557 0.863982 0.641557  0.969889 0.382148 0.508773 
0.932684 0.531116 0.838293 0.031451  0.242338 0.663352 0.784813 
0.283031 0.754613 0.938358 0.0408097  0.609105 0.325545 0.671151 

julia> b = rand(10) 
10-element Array{Float64,1}: 
0.0795157 
0.219318 
0.965155 
0.896807 
0.701626 
0.741823 
0.954437 
0.573683 
0.493615 
0.0821557 


julia> A\b 
10-element Array{Float64,1}: 
    1.47909 
    2.39816 
-0.15789 
    0.144003 
-1.10083 
-0.273698 
-0.775122 
    0.590762 
-0.0266894 
-2.36216 

können Sie @which zu sehen, wie es definiert ist:

julia> @which A\b 
\(A::AbstractArray{T,2} where T, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) in Base.LinAlg at linalg\generic.jl:805 

Was uns hier führt: https://github.com/JuliaLang/julia/blob/master/base/linalg/generic.jl#L827 (Zeilennummern leicht wegen Versionsunterschiede ändern). Wie Sie sehen können, führt es einige schnelle Funktionsaufrufe durch, um festzustellen, um welchen Matrixtyp es sich handelt. istril findet aus seinem unteren Dreieck: https://github.com/JuliaLang/julia/blob/master/base/linalg/generic.jl#L987, usw. Sobald es den Matrix-Typ bestimmt, spezialisiert es die Matrix so viel wie möglich, so dass es effizient sein kann, und ruft dann \. Diese spezialisierten Matrixtypen führen entweder eine Faktorisierung durch, die dann die Rückwärtssubstitution durchführt (was eine gute Möglichkeit ist, \ auf Ihrem eigenen BTW zu verwenden, um die Faktorisierung erneut zu verwenden), oder sie "weiß" direkt die Antwort, wie für Dreiecks- oder Diagonalmatrizen .

Konkreter als die Quelle kann nicht erhalten werden.

Beachten Sie, dass \ etwas anders als nur invertieren ist. Normalerweise möchten Sie keine Matrix invertieren, geschweige denn eine große Matrix. Diese Faktorisierungen sind viel numerischer stabil. inv wird jedoch eine Inversion durchführen, die einer LU-Faktorisierung sehr ähnlich ist (was in Julia lufact ist). Sie können auch in pinv für die Psudo-Inverse in einigen Fällen suchen, wo die Matrix Singular oder nahe Singular ist, aber Sie sollten dies wirklich vermeiden und stattdessen das System faktorisieren + lösen anstatt das Inverse.

Für sehr große dünn besetzte Matrizen sollten Sie iterative Solver verwenden. Sie finden viele Implementierungen in IterativeSolvers.jl

+5

Eine 10^6 x 10^6 dichten 'Float64'-Matrix benötigt 7,3 Terabyte RAM. Da das Invertieren einer dünnen Matrix eine dichte Matrix erzeugt, möchten Sie Ihre Matrix nicht wirklich invertieren. Dies ist eine gute Faustregel: Nicht invertieren, lösen. – StefanKarpinski

+0

Danke für die Hilfe. Ich arbeite jetzt mit dem Backslash-Operator. Ich probiere auch die inv() -Aufruf, aber es gibt mir Stackoverflow-Fehler. Ich suchte nach einem Mechanismus, um das System-RAM, das von Julia benutzt werden soll, etwas zu begrenzen (Mein System hat einen RAM von 1TB) Ich frage mich, ob wir solche Sprachbeschränkungen in Julia haben. –

Verwandte Themen