Normalgleichungen lösen

Das gelöste Problem

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)

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

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

Pseudoinverse nach Moore und Penrose

Die Moore-Penrose-Pseudoinverse ist für jede reelle Matrix definiert, auch für rechteckige Matrizen und Matrizen ohne vollen Rang. Schreibe ihre Singulärwertzerlegung als

A=U\Sigma V^T

Die dargestellte Faktorisierung ist eine Singulärwertzerlegung. Für die Pseudoinverse 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 Zerlegung mit QR

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

A=Q_1R_1,\qquad Q_1^TQ_1=I

Die Zielfunktion des Ausgleichsproblems zerfällt in

\begin{aligned} \|Ax-b\|_2^2 &=\|R_1x-Q_1^Tb\|_2^2 \\ &\quad{}+\|(I-Q_1Q_1^T)b\|_2^2 \end{aligned}

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

Die folgenden numerischen Faktoren und Koeffizienten sind $auf fünf Nachkommastellen gerundet$ . Die exakten Matrixgleichungen gelten vor dem Runden.

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 $\mathrm{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 $\mathrm{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$ :

\begin{aligned} &-1{,}73205x_3=3{,}48060 \\ &\Longrightarrow \\ &x_3=-2{,}00952 \end{aligned}

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

Die cyanfarbene Kurve ist das angepasste Polynom. Die orangefarbenen Punkte sind die sieben Messwerte aus

b

$x$: $[-3; 3]$
$y$: $[-4{,}35; 15{,}9]$

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 $O(mn^2)$ bilden, dann in etwa $\frac{1}{3}n^3$ Flops faktorisieren	schnell bei gut konditioniertem Problem und vollem Rang	$\kappa_2(A^TA)=\kappa_2(A)^2$
Householder- $\mathrm{QR}$	etwa $2mn^2-\frac{2}{3}n^3$ Flops	bildet $A^TA$ nicht und ist meist stabiler	bei abhängigen Spalten ist eine rangoffenbarende Variante nötig
$\mathrm{SVD}$	aufwendiger als $\mathrm{QR}$	gibt den Rang an und liefert die $Lösung kleinster Norm$	höchste Kosten der drei Verfahren

Verwende $\mathrm{QR}$ als übliche Wahl für dichte Probleme mit vollem Rang. Verwende die $\mathrm{SVD}$ , wenn die Matrix rangdefizient oder stark fehlkonditioniert ist oder eine Lösung kleinster Norm benötigt wird. Bilde die Normalmatrix nur, wenn der Geschwindigkeitsgewinn den Verlust an numerischer Stabilität rechtfertigt. Berechne in Produktionscode keine Inverse explizit.