Normalgleichungen lösen

Welches Problem wird gelöst?

Für $A \in \mathbb{R}^{m \times n}$ und $b \in \mathbb{R}^m$ sucht ein lineares Ausgleichsproblem den Parametervektor, der das Residuum möglichst klein macht:

\min_x \|Ax-b\|_2^2

Jeder Minimierer erfüllt die Normalgleichung

A^T A\hat{x}=A^T b

Diese Gleichung ist immer lösbar. Der Minimierer ist genau dann eindeutig, wenn $\operatorname{rank}(A)=n$ gilt. Das setzt $m\ge n$ voraus. Bei vollem Spaltenrang gilt

\hat{x}=(A^T A)^{-1}A^T b

Hat $A$ keinen vollen Spaltenrang, bleibt der angepasste Vektor $A\hat{x}$ eindeutig, der Parametervektor jedoch nicht. Ist $\hat{x}$ eine Lösung, dann liegen alle Lösungen in

\hat{x}+\ker(A)

Warum gibt es immer eine Lösung?

Der Spaltenraum von $A$ ist ein endlichdimensionaler abgeschlossener Unterraum von $\mathbb{R}^m$ . Deshalb besitzt $b$ eine eindeutige orthogonale Projektion auf diesen Raum. Jedes $\hat{x}$ mit $A\hat{x}=\operatorname{proj}_{\operatorname{col}(A)}b$ ist eine Ausgleichslösung.

Die Normalgleichung beschreibt dieselbe Geometrie, denn

A^T(A\hat{x}-b)=0

bedeutet, dass das Residuum $A\hat{x}-b$ zu jeder Spalte von $A$ orthogonal ist.

Moore-Penrose-Pseudoinverse

Die Moore-Penrose-Pseudoinverse ist für jede reelle Matrix definiert, auch für rechteckige Matrizen und Matrizen ohne vollen Rang. Ist

A=U\Sigma V^T

eine Singulärwertzerlegung, wird jeder von null verschiedene Singulärwert $\sigma_i$ durch $1/\sigma_i$ ersetzt. Singulärwerte gleich null bleiben null. Damit erhält man

\Sigma^{\dagger}_{ii}=\begin{cases}1/\sigma_i,&\sigma_i>0,\\0,&\sigma_i=0,\end{cases}

Der Vektor $x^{\dagger}$ ist die eindeutige Ausgleichslösung mit kleinster euklidischer Norm. Hat $A$ vollen Spaltenrang, vereinfacht sich die allgemeine Definition zu

A^{\dagger}=(A^T A)^{-1}A^T

Die Pseudoinverse wird durch die vier Penrose-Identitäten eindeutig bestimmt:

AA^{\dagger}A=A

Direkte Lösung mit der QR-Zerlegung

Für eine Matrix mit vollem Spaltenrang und $m\ge n$ verwendet man die dünne QR-Zerlegung

A=Q_1R_1,\qquad Q_1^TQ_1=I

Die Zielfunktion des Ausgleichsproblems zerfällt in

\|Ax-b\|_2^2=\|R_1x-Q_1^Tb\|_2^2+\|(I-Q_1Q_1^T)b\|_2^2

Der zweite Term hängt nicht von $x$ ab. Daher erhält man den Minimierer durch Lösen des oberen Dreieckssystems

R_1\hat{x}=Q_1^Tb

Dieser Weg löst das ursprüngliche Ausgleichsproblem direkt. Man bildet nicht zuerst $A^TA$ , um anschließend die Normalgleichung zu lösen.

Beispiel einer Polynomanpassung

Wir passen $y=at^2+bt+c$ an sieben Messwerte an. Die Designmatrix und der Beobachtungsvektor lauten

A=\begin{pmatrix}9&-3&1\\4&-2&1\\1&-1&1\\0&0&1\\1&1&1\\4&2&1\\9&3&1\end{pmatrix}

Jede Zeile von $A$ hat die Form $[t_i^2,t_i,1]$ . Eine gültige dünne QR-Zerlegung ist

Q_1=\begin{pmatrix}-0{,}64286&-0{,}56695&0{,}16496\\-0{,}28571&-0{,}37796&-0{,}24744\\-0{,}07143&-0{,}18898&-0{,}49487\\0&0&-0{,}57735\\-0{,}07143&0{,}18898&-0{,}49487\\-0{,}28571&0{,}37796&-0{,}24744\\-0{,}64286&0{,}56695&0{,}16496\end{pmatrix}

Die Vorzeichen zusammengehöriger Spalten von $Q_1$ und Zeilen von $R_1$ können sich zwischen gültigen QR-Implementierungen unterscheiden. Das Produkt $Q_1R_1$ und die endgültige Anpassung bleiben gleich.

Zuerst berechnen wir

Q_1^Tb=\begin{pmatrix}-9{,}71429\\16{,}02569\\3{,}48060\end{pmatrix}

Danach folgt die Rückwärtssubstitution:

-1{,}73205x_3=3{,}48060\quad\Longrightarrow\quad x_3=-2{,}00952

Damit gilt

\hat{x}=\begin{pmatrix}0{,}98095\\3{,}02857\\-2{,}00952\end{pmatrix}

Ausgleichskurve

\hat{y}(t)=0{,}98095t^2+3{,}02857t-2{,}00952

Wahl des numerischen Verfahrens

Verfahren	Ungefähre Kosten für dichte Matrizen mit $m\ge n$	Hauptvorteil	Wichtigste Einschränkung
Normalgleichungen und Cholesky	$A^TA$ in bilden, dann in etwa Flops faktorisieren

QR ist die übliche Wahl für dichte Probleme mit vollem Rang. Verwende die SVD, wenn Rangmangel, starke Fehlkonditionierung oder die Lösung kleinster Norm wichtig ist. Bilde die Normalmatrix nur, wenn ihr Geschwindigkeitsvorteil den Verlust an numerischer Robustheit rechtfertigt, und berechne in Produktionscode niemals explizit die Inverse.