Kaczmarz-Methode

Die auf der Arbeit des polnischen Mathematikers Stefan Kaczmarz basierte Kaczmarz-Methode dient der iterativen Lösung linearer Gleichungssysteme der Form <math>Ax = b</math>, wobei <math>A\in\R^{m\times n}</math> eine, evtl. nicht-quadratische, Matrix, b die gegebene rechte Seite und x der gesuchte Lösungsvektor ist. Es handelt sich dabei um einen iterativen Algorithmus, der unter anderem Einzug in die Computertomographie und digitale Signalverarbeitung gefunden hat.

Im Gegensatz zu den meisten anderen iterativen Lösungsverfahren, wie zum Beispiel dem Gauß-Seidel-Verfahren, benötigt der Kaczmarz-Algorithmus keine invertierbare Matrix A. Insbesondere im unter-bestimmten Fall mit <math>m<n</math> lässt sich damit die Lösung <math>x^+=A^+ b</math> kleinster Norm zur Pseudoinversen <math>A^+</math> berechnen. Aus diesem Grund kann das Verfahren für praktisch alle Anwendungen problemlos eingesetzt werden, obwohl Verfahren für spezielle Problemstellungen wesentlich schneller sein können. Allerdings zeigt eine randomisierte Variante des Kaczmarz-Verfahrens wesentlich bessere Konvergenz im Falle überbestimmter Gleichungssysteme als andere iterative Lösungsverfahren.

Der ursprüngliche Algorithmus

Gegeben ist eine reelle (oder auch komplexe) m × n Matrix <math>A\in\R^{m\times n}</math> von vollem Rang, sowie ein Vektor <math>b=(b_i)\in\R^m</math>. Die folgende Iterationsformel verbessert bei jedem Durchlauf die Annäherung <math>x_k</math> an eine Lösung x der Gleichung Ax = b:

<math>

 x_{k+1} = x_{k} + \omega_k \frac{b_{i} - a_{i}^T x_{k} }{\lVert a_{i} \rVert^2} a_{i}.

</math>, Bei der deterministischen Variante wird der Index i zyklisch gewählt,

<math>

i = (k \, \bmod \, m) + 1 </math>

und <math>\omega_k\in\R</math> ist ein positiver Relaxationsparameter, der in der ursprünglichen Version gleich eins war. Dabei ist mit <math>a_{i}\in\R^n</math> die transponierte i-te Zeile der Matrix A gemeint, <math>{\lVert a_{i} \rVert^2}</math> ist die quadrierte euklidische Norm von <math> a_{i} </math>, die Summe der quadrierten Einträge von <math>a_{i}</math> bzw. das Skalarprodukt <math>\langle a_{i},a_{i} \rangle=a_i^T a_i</math>. Das Verfahren ist besonders effizient bei dünn-besetzter Matrix, wenn die Vektoren <math>a_i</math> nur wenige nichttriviale Elemente besitzen.

Geometrisch projiziert der einzelne Iterationsschritt den aktuellen Näherungsvektor <math>x_k</math> orthogonal auf die i-te Hyperebene, es gilt <math>a_i^T x_{k+1}=b_i.</math>

Mit dem Startwert <math>x_0=0</math> konvergiert die Iteration im unterbestimmten Fall <math>m<n</math> gegen die Lösung <math>x^+=A^+ b</math> kleinster Norm zur Pseudoinversen <math>A^+</math>. Im überbestimmten, inkonsistenten Fall <math>m>n</math> springen die Iterierten am Ende allerdings zwischen den Hyperebenen hin und her, und man bekommt nur bei starker Unter-Relaxation mit <math>\omega_k\ll 1</math> eine brauchbare Näherung an die Kleinste-Quadrate-Lösung <math>x^+=A^+b</math>.

Querverbindung zu anderen Iterationsverfahren

Im unterbestimmten Fall <math>m<n</math> mit vollem Rang von A liegt die Lösung minimaler Norm im orthogonalen Komplement des Kerns/Nullraums <math>x^+\in N(A)^\perp=B(A^T)</math>, der mit dem Bildraum der transponierten Matrix übereinstimmt. Daher lässt sich <math>x^+</math> darstellen als <math>x^+=A^T y^+,\,y^+\in\R^m</math>. Somit ist <math>x^+</math> die eindeutige Lösung des Gleichungssystems

<math>

AA^T y=b,\quad x=A^T y,

</math>

mit positiv definiter Matrix <math>AA^T</math>. Wenn man die Iteration mit dem Startwert <math>x_0=0</math> beginnt, liegen offensichtlich auch alle Iterierten im Bildraum von <math>A^T</math>. Daher hat jede Iterierte eine eindeutige Darstellung <math>x_k=A^Ty_k,\,y_k\in\R^m</math>. Mit dieser wird der Iterationsschritt zu

<math>

 A^Ty_{k+1} = A^T\big(y_{k} + \omega_k \frac{b_{i} - a_{i}^T A^Ty_{k} }{\lVert a_{i} \rVert^2}e_{i}\big) ,

</math>

wobei <math>e_i</math> den i-ten Einheitsvektor bezeichnet. Da A vollen Rang besitzt, bedeutet das

<math>

 y_{k+1} = y_{k} + \omega_k \frac{b_{i} - a_{i}^T A^Ty_{k} }{\lVert a_{i} \rVert^2} e_{i}.

</math>

Dies ist aber genau der Teilschritt in einer einzelnen Komponente für das Einzelschritt-/Gauß-Seidel-Verfahren (<math>\omega_k=1</math>) beziehungsweise mit allgemeinem <math>\omega_k</math> das SOR-Verfahren zum System <math>AA^T y=b</math>, deren wohlbekannte Konvergenzaussagen sich also auf die Kaczmarz-Iteration übertragen. Z.B. konvergiert diese für festes <math>\omega_k\equiv\omega\in(0,2)</math>.

Randomisierter Algorithmus

Wie bereits weiter oben erwähnt, gibt es eine Variante des Verfahrens, die im Falle überbestimmter Gleichungssysteme wesentlich schneller konvergiert.<ref> Thomas Strohmer, Roman Vershynin: A randomized solver for linear systems with exponential convergence. (PDF; 155 kB)</ref> Dabei wird bessere Konvergenz nicht nur für das Lösen hochgradig überbestimmter Gleichungssysteme erreicht, sondern bspw. auch beim Lösen von Systemen mit 10 % mehr Gleichungen als Unbekannten. Hierzu muss bei der obigen Iterationsgleichung das i nicht abhängig von k, sondern zufällig (daher der Name der Methode) gewählt werden. Die Wahrscheinlichkeit für ein bestimmtes i ist hierbei proportional zu <math> \lVert a_{i} \rVert ^2 </math> für jede Iteration.

Ungleichungssysteme

Eine einfache Modifikation des Verfahrens kann zulässige Lösungen für Ungleichungssysteme <math>Ax\ge b,\,m\ge n</math>, berechnen. Beim Schritt

<math>

 x_{k+1} = x_{k} + \frac{\max\{0,b_{i} - a_{i}^T x_{k}\}}{\lVert a_{i} \rVert^2} a_{i},

</math>

wird eine Änderung nur durchgeführt, wenn <math>b_{i} - a_{i}^T x_{k}>0</math> ist, also die i-te Ungleichung <math>a_i^T x_k\ge b_i</math> verletzt ist. Dann wird <math>x_k</math> von außen auf den i-ten Halbraum projiziert.

Auf diese Weise lassen sich mit dem Verfahren sogar Lineare Programme lösen über die Äquivalenz von Optimierungs- und Zulässigkeitsproblemen.

Quellen

S. Kaczmarz: Przyblizone rozwiazywanie ukladów równan liniowych. – Angenäherte Auflösung von Systemen linearer Gleichungen. Bulletin International de l’Académie Polonaise des Sciences et des Lettres. Classe des Sciences Mathématiques et Naturelles. Série A, Sciences Mathématiques, S. 355–357, 1937. (Achtung: Dieser Artikel wird sehr verbreitet mit der falschen Seitenangabe 335–357 zitiert.)
M. Hanke-Bourgeois; W. Niethammer: On the acceleration of Kaczmarz's method for inconsistent linear systems; Linear Algebra Applcs 130 (1990), 83–98
A.Björck, T. Elfving: Accelerated projection methods for computing pseudoinverse solutions of systems of linear equations, BIT 19 (1979), 145–163
T. Strohmer, R. Vershynin: A randomized Kaczmarz algorithm with exponential convergence, J. Fourier Anal. Appl. 15(1), 262–278, 2009
W. Hackbusch: Iterative Lösung großer schwachbesetzter Gleichungssysteme. Teubner Studienbücher, 1993, S. 202–203

Einzelnachweis