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

$Ax=b$

mit regulärer, jedoch schlecht konditionierter Systemmatrix $A\in \mathbb {R} ^{n\times n}$ und rechter Seite $b\in \mathbb {R} ^{n}$.

Die numerische Lösung ist fehlerbehaftet. Das kann so interpretiert werden, dass man ein Gleichungssystem

${\tilde {A}}{\tilde {x}}=b$

mit gestörter Systemmatrix ${\tilde {A}}$ löst. Die Lösungen des originalen und des gestörten Gleichungssystems $x$ bzw. ${\tilde {x}}$ weichen um einen Korrekturterm $\Delta x$ voneinander ab:

$x={\tilde {x}}+\Delta x$

Setzt man ${\tilde {x}}+\Delta x$ für $x$ in die Originalgleichung ein erhält man

$A({\tilde {x}}+\Delta x)=b$.

Auflösen nach dem Korrekturterm $\Delta x$ liefert:

$A\Delta x=b-A{\tilde {x}}$

Geht man davon aus, dass man die als Residuum bezeichnete rechte Seite $r:=b-A{\tilde {x}}$ relativ genau berechnen kann, bei dem Lösen des linearen Gleichungssystems bzgl. $\Delta x$ jedoch nur mit der gestörten Systemmatrix ${\tilde {A}}$ rechnen kann, so erhält man das folgende iterative Verfahren:

• Startwert: ${\tilde {x}}^{(0)}=0\in \mathbb {R} ^{n}$
• Für $k=0,\ldots ,N$:

1. Berechne das Residuum $r^{(k)}=b-A{\tilde {x}}^{(k)}$ in hoher Genauigkeit.
2. Berechne die Korrektur $\Delta x^{(k)}$ als numerische Lösung des Gleichungssystems ${\tilde {A}}\Delta x^{(k)}=r^{(k)}$
3. Berechne die verbesserte numerische Lösung ${\tilde {x}}^{(k+1)}={\tilde {x}}^{(k)}+\Delta x^{(k)}$ des Originalsystems.
4. Abbruch für hinreichend kleinen Korrekturterm $|\Delta x^{(N)}|$ und hinreichend kleines Residuum $|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.

Kategorie:Numerische lineare Algebra