Benutzer:TN/Iterative Refinement
Im Regelfall weicht die numerische Lösung eines linearen Gleichungssystems durch Rundungsfehler von der exakten Lösung ab. Diese Abweichung kann bei schlecht konditionierter Systemmatrix inakzeptabel groß sein. Mit Hilfe iterativer Verfahren, die meist mit dem englischen Begriff Iterative Refinement bezeichnet werden, lässt sich die numerische Lösung verbessern.
Grundverfahren
Zu lösen ist ein lineares Gleichungssystem
- Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle Ax = b}
mit regulärer, jedoch schlecht konditionierter Systemmatrix Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle A\in\mathbb{R}^{n\times n}} und rechter Seite Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle b\in\mathbb{R}^{n}} .
Die numerische Lösung ist fehlerbehaftet. Das kann so interpretiert werden, dass man ein Gleichungssystem
- Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \tilde A\tilde x=b}
mit gestörter Systemmatrix Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \tilde A} löst. Die Lösungen des originalen und des gestörten Gleichungssystems Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle x} bzw. Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \tilde x} weichen um einen Korrekturterm Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \Delta x} voneinander ab:
- Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle x=\tilde x + \Delta x}
Setzt man Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \tilde x + \Delta x} für Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle x} in die Originalgleichung ein erhält man
- Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle A(\tilde x + \Delta x) = b} .
Auflösen nach dem Korrekturterm Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \Delta x} liefert:
- Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle A\Delta x = b-A\tilde x}
Geht man davon aus, dass man die als Residuum bezeichnete rechte Seite Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle r:= b-A\tilde x} relativ genau berechnen kann, bei dem Lösen des linearen Gleichungssystems bzgl. Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \Delta x} jedoch nur mit der gestörten Systemmatrix Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \tilde A} rechnen kann, so erhält man das folgende iterative Verfahren:
- Startwert: Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \tilde x^{(0)}=0\in\mathbb{R}^n}
- Für Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle k=0,\ldots,N} :
- Berechne das Residuum Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle r^{(k)}=b-A\tilde x^{(k)}} in hoher Genauigkeit.
- Berechne die Korrektur Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \Delta x^{(k)}} als numerische Lösung des Gleichungssystems Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \tilde A\Delta x^{(k)} = r^{(k)}}
- Berechne die verbesserte numerische Lösung Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \tilde x^{(k+1)} = \tilde x^{(k)} + \Delta x^{(k)}} des Originalsystems.
- Abbruch für hinreichend kleinen Korrekturterm Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle |\Delta x^{(N)}|} und hinreichend kleines Residuum Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle |r^{(N)}|} .
Das folgende Scilab-Programm demonstriert die iterative Lösungsverbesserung (das ist keine numerisch effiziente Implementation, sondern dient ausschließlich der Demonstration des Verfahrens):
// Eingangsdaten:
N = 10; // Systemdimension
Cond = 1e7; // Kondition der Systemmatrix
noise = 0.01; // Störung durch LU-Zerlegung
maxIter = 10; // Begrenzung der Anzahl der Iterationsschritte
epsilon = 1e-8; // Genauigkeit mit der die Lösung berechnet werden soll.
// Generierung der Problemdaten:
// Schlecht konditionierte Systemmatrix generieren:
[U,S,VT] = svd(rand(N,N));
logCond = log(Cond);
S = exp(0:logCond/(N-1):logCond);
A = U * diag(S) * VT;
Anoisy = U * diag((1.0+noise*rand(1,N)).*S)*VT;
// Rechte Seite generieren:
b = rand(N,1);
// Eigentliches Verfahren "Iterative Refinement":
x = zeros(N,1); // Startwert für Iteration
for i=1:maxIter do
r = b - A*x; // Berechnung des Residuums
dx = Anoisy\r; // Lösung des linearen Gleichungssystems für den Korrekturterm dx
// (An dieser Stelle werden die Rundungsfehler bei der Lösung des
// linearen Gleichungssystems durch Rechnen mit der gestörten Matrix
// Anoisy simuliert.)
x = x + dx; // Verbesserung der numerischen Lösung des Originalproblems
norm_dx = norm(dx);
norm_r = norm(r);
print( %io(2), i, norm_dx, norm_r );
// Abbruchstest:
if norm_dx + norm_r < epsilon then
break;
end;
end;
Anwendung
In der LAPACK-Routine dsgesv wird das oben vorgestellte Grundverfahren für Iterative Refinement zur schnellen Lösung linearer Gleichungssysteme eingesetzt. Als eingebetteter direkter Lösungsalgorithmus kommt dabei eine LR-Zerlegung mit einfacher Genauigkeit zum Einsatz, während die nicht so aufwendige Residuenberechnung mit doppelter Genauigkeit durchgeführt wird. Die LR-Zerlegung wird nur einmal durchgeführt und innerhalb der Iterationsschleife werden die gewonnenen Dreiecksmatrizen zur Lösung des linearen Gleichungssystems genutzt.
Weitere Verfahren
Iterative Lösungsverfahren für lineare Gleichungssysteme (wie zum Beispiel das PCG-Verfahren) sind für Iterative Refinement einsetzbar. Die numerische Lösung des Gleichungssystems mittels LR-Zerlegung oder Cholesky-Verfahren dient dann als Präkonditionierer.
Literatur
- Germund Dahlquist, Ake Björck: Numerical Methods in Scientific Computing, Volume II. siam, Philadelphia 2007, [1].
- W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery: Nuerical Recipes in C. Second Edition, Cambridge University Press 1992, ISBN: 0 521 43108 5.
