Verfasst am: 21.04.2014, 13:03
Titel: Effizienz: Matrixberechnung mit Funktion abhängig vom Index
Hallo zusammen,
ich muss für meine Arbeit sehr große NxN Matrizen (N=10^4 bis 10^6) aufstellen und anschließend invertieren.
(Thematisch geht es darum, einen durch nicht-verschwindenden Anteil eines Zentralspins in einem Spinbad (Größenordnung N) zu berechnen, um wichtige Informationen für die Speicherung von Information bei der Realisierung von Quantencomputern zu gewinnen.)
Die einzelnen Elemente müssen, abhängig von ihrem Indize, durch unterschiedliche Funktionen berechnet werden. Die Aufteilung ist wie folgt:
1. erstes Element
2. erste Zeile
3. Diagonale
4. restliche obere Dreiecksmatrix
Da die Matrix hermitesch ist, benötige ich nur die obere Dreiecksmatrix und kann den Rest dann einfach kopieren.
Momentan löse ich das Problem mit vielen for-Schleifen, was aber nicht effizient sein kann. Effizienz ist jedoch genau das, was ich benötige, da bereits bei N=1000 etwa 9 Minuten für die Berechnung benötigt werden.
Ich suche also nach einer Möglichkeit, meine Matrix möglichst ohne for-Schleifen zu berechnen.
J_vector = 1:n; %Kopplungsvektor proportional zu k
M(1,1)=(n+1)/4; %erstes Element
for l=1:n %Diagonale
y=S_l(l,J_vector,n);
M(l+1,l+1)=(2*y*y + 3*(n-1)*Q_l(l,J_vector,n))/64;
end
for l=1:n %erste Zeile
M(1,l+1)=S_l(l,J_vector,n)/8;
M(l+1,1)=M(1,l+1);
end
for i=1:(n-1) %obere Hälfte (ohne erste Zeile) %Zeilennummer
for j=(i+1):n %Spaltennummer
x= J_j_l(J_vector(i),J_vector(j));
M(i+1,j+1)=(x * (S_l(j,J_vector,n) - S_l(i,J_vector,n)))/16 - (3*(n-3) * x*x)/64;
M(j+1,i+1)=M(i+1,j+1);
end end end
function[J] = J_j_l(x_l, x_j) %berechnet J_j^(l)
J=1/(-1/x_l + 1/x_j);
end
function[S] = S_l(l,J_vector,n) %berechnet S^(l)
S=-J_vector(l); %Epsilon_0 = 0 for j=1:n
if(j~=l) %da i =/= l sein muss
S=S+J_j_l(J_vector(l),J_vector(j));
end end end
function[Q] = Q_l(l,J_vector,n) %berechnet Q^(l)
Q =J_vector(l)*J_vector(l); %Epsilon_0 = 0 for j=1:n
if(j~=l) %da i =/= l sein muss
k = J_j_l(J_vector(l),J_vector(j));
Q=Q+k*k;
end end end
Mein neuer Ansatz zur Realisierung wäre bsxfun, jedoch weiß ich nicht wie ich die Methode hier genau einbauen soll, da die Funktion nicht immer die gleiche ist.
Eine Idee wäre jetzt die Funktion mit if-Bedingungen so zu gestalten, dass sie abfragt, welches Element gerade berechnet wird, und dann die passende Berechnung ausführt. Dann weiß ich jedoch nicht, wie ich ausnutzen kann, dass die Matrix hermitesch ist, ich also quasi nur die Hälfte der Matrix berechnen muss.
Über weitere Vorschläge zur Effizienzsteigerung würde ich mich auch freuen. Versuche, die Summenbildung in S_l und Q_l ohne for-Schleife durchzuführen, führten leider nur zu einer Verschlechterung der Effizienz. Die Methode über die for-Schleife scheint wohl die schnellste zu sein.
Verfasst am: 21.04.2014, 22:26
Titel: Re: Effizienz: Matrixberechnung mit Funktion abhängig vom I
Hallo Philsch,
Zitat:
ich muss für meine Arbeit sehr große NxN Matrizen (N=10^4 bis 10^6) aufstellen und anschließend invertieren.
Hier würde ich gleich mal einhaken:
Eine 1e6 x 1e6 Matrix benötigt 8TB, wenn Du DOUBLE als Typ verwendest. Zum Invertieren würde ich mindestens nochmal den gleichen Speicher. Die explizite Berechnung der Inversen -übrigens eine sehr sehr seltene Aufgabe in echten Anwendungen- per LU-Dekomposition DGETRF und danach DGETRI würde wohl noch mehr Speicher benötigen.
Hast Du eine > 20TB Ram Maschine zur Verfügung? Wenn nicht, liegt hier bereits ein ernstes Design-Problem vor.
Sind die benötigten Matrizen dünn besetzt? Dann wären Sparse-Matrizen unbedingt zu empfehlen!
Der Aufruf von Unterfunktionen benötigt viel Zeit. Hier kann man viel erreichen:
Code:
function S = S_l(l, J_vector,n) %berechnet S^(l)
S = -J_vector(l); %Epsilon_0 = 0 for j=1:n
if(j~=l) %da i =/= l sein muss
% S = S + J_j_l(J_vector(l), J_vector(j));
S = S + 1 / (1/J_vector(j) - 1/J_vector(l));
end end end
Ich habe nur den Aufruf von J_j_l durch die explizite Code-Zeile ersetzt und das Programm läuft 10 mal schneller, für 1000 komme ich auf meiner betagten Core2Duo-Maschine auf 52 Sekunden.
Dann muss man natürlich "1/J_vector(l)" nicht in jeder Iteration neu berechnen: Einmal vor der Schleife und dann in einer Variablen speichern. Das bringt weitere 20% mehr Geschwindigkeit.
Das gleiche nochmal in der Hauptfunktion:
Code:
... for i=1:(n-1) %obere Hälfte (ohne erste Zeile) %Zeilennummer
c1 = S_l(i,J_vector,n);
c2 = 1/J_vector(i);
for j=(i+1):n %Spaltennummer
% x= J_j_l(J_vector(i),J_vector(j));
x = 1 / (1/J_vector(j) - c2);
M(i+1,j+1) = (x * (S_l(j, J_vector, n) - c1))/16 - (3*(n-3) * x*x)/64;
M(j+1,i+1) = M(i+1, j+1);
end end
Es sieht so aus, als würde sich die Laufzeit um den Faktor 5 verlängern, wenn man den Input doppelt so groß wählt. Grob überschlagen käme ich dann auf 2700 Tage Laufzeit für die 1e6 x 1e6 Matrix, wenn man denn genug RAM hätte. Vielleicht lässt sich noch ein Faktor 10 rausholen, wenn man ganz schlau vorgeht.
Mit weniger RAM, wenn man die Festplatte als virtuellen Speicher verwenden muss, kommt wieder ein Faktor 100 oder 1000 oben drauf.
Erstmal riesigen Dank, das Programm läuft jetzt tatsächlich 30 mal schneller. n=1000 in 14s, so seh ich das gerne Da lernt man doch immer, dass man schön Funktionen nutzen soll, und dann sowas - hat mich gut überrascht.
Zum Glück muss ich nicht wirklich 10^6 Matrizen aufstellen, das hab ich leider unglücklich formuliert. Tatsächlich ist das zwar die Anzahl an Badspins, die eine Rolle spielt, aber das Ergebnis scheint ohnehin stark zu konvergieren (jedenfalls für diesen Kopplungsvektor, aber für den neuen, der proportional zu exp(-k) ist wohl auch).
Außerdem verwende ich später eine Methode, um auf N=inf zu schließen. Wenn ich jetzt also mit N=10^4 jetzt arbeiten kann, dann bin ich gut zufrieden, und ich denke mein Prof auch.
Die Matrix ist leider nicht dünn besetzt sondern ziemlich voll besetzt.
Einen eleganten Weg, um die ganzen for-Schleifen zu umgehen, siehst du aber auch nicht, oder?
Vielleicht kann man an meinem Skript (dort findet die Invertierung statt) aber auch noch etwas verbessern, das sieht wie folgt aus:
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
Da lernt man doch immer, dass man schön Funktionen nutzen soll, und dann sowas - hat mich gut überrascht.
RSS