Kaczmarz-Methode

Die auf der Arbeit des polnischen Mathematikers Stefan Kaczmarz basierte Kaczmarz-Methode dient der iterativen Lösung linearer Gleichungssysteme der Form ${\displaystyle Ax=b}$, wobei ${\displaystyle A\in \mathbb {R} ^{m\times n}}$ 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 ${\displaystyle m<n}$ lässt sich damit die Lösung ${\displaystyle x^{+}=A^{+}b}$ kleinster Norm zur Pseudoinversen ${\displaystyle A^{+}}$ 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 ${\displaystyle A\in \mathbb {R} ^{m\times n}}$ von vollem Rang, sowie ein Vektor ${\displaystyle b=(b_{i})\in \mathbb {R} ^{m}}$. Die folgende Iterationsformel verbessert bei jedem Durchlauf die Annäherung ${\displaystyle x_{k}}$ an eine Lösung x der Gleichung Ax = b:

${\displaystyle x_{k+1}=x_{k}+\omega _{k}{\frac {b_{i}-a_{i}^{T}x_{k}}{\lVert a_{i}\rVert ^{2}}}a_{i}.}$,

Bei der deterministischen Variante wird der Index i zyklisch gewählt,

${\displaystyle i=(k\,{\bmod {\,}}m)+1}$

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

Geometrisch projiziert der einzelne Iterationsschritt den aktuellen Näherungsvektor ${\displaystyle x_{k}}$ orthogonal auf die i-te Hyperebene, es gilt ${\displaystyle a_{i}^{T}x_{k+1}=b_{i}.}$

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

Querverbindung zu anderen Iterationsverfahren

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

${\displaystyle AA^{T}y=b,\quad x=A^{T}y,}$

mit positiv definiter Matrix ${\displaystyle AA^{T}}$. Wenn man die Iteration mit dem Startwert ${\displaystyle x_{0}=0}$ beginnt, liegen offensichtlich auch alle Iterierten im Bildraum von ${\displaystyle A^{T}}$. Daher hat jede Iterierte eine eindeutige Darstellung ${\displaystyle x_{k}=A^{T}y_{k},\,y_{k}\in \mathbb {R} ^{m}}$. Mit dieser wird der Iterationsschritt zu

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

wobei ${\displaystyle e_{i}}$ den i-ten Einheitsvektor bezeichnet. Da A vollen Rang besitzt, bedeutet das

${\displaystyle y_{k+1}=y_{k}+\omega _{k}{\frac {b_{i}-a_{i}^{T}A^{T}y_{k}}{\lVert a_{i}\rVert ^{2}}}e_{i}.}$

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

Randomisierter Algorithmus

Wie bereits weiter oben erwähnt, gibt es eine Variante des Verfahrens, die im Falle überbestimmter Gleichungssysteme wesentlich schneller konvergiert. 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 ${\displaystyle \lVert a_{i}\rVert ^{2}}$ für jede Iteration.

Ungleichungssysteme

Eine einfache Modifikation des Verfahrens kann zulässige Lösungen für Ungleichungssysteme ${\displaystyle Ax\geq b,\,m\geq n}$, berechnen. Beim Schritt

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

wird eine Änderung nur durchgeführt, wenn ${\displaystyle b_{i}-a_{i}^{T}x_{k}>0}$ ist, also die i-te Ungleichung ${\displaystyle a_{i}^{T}x_{k}\geq b_{i}}$ verletzt ist. Dann wird ${\displaystyle x_{k}}$ 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