3 Lineare Ausgleichsrechnung
Anstatt eine extakte Lösung zu finden, können wir auch eine Näherungslösung suchen, die die Summe der Abweichungen minimiert. Dies wird als lineare Ausgleichsrechnung bezeichnet.
3.1 Normalengleichung
\[ A^\top Ax = A^\top b \]
ist äquivalent zur Minimierung von \(\|Ax - b\|_2\)
- Wenn \(A^\top A\) invertierbar ist, ist die Lösung eindeutig.
- Wenn \(\operatorname{Kern} \{v \in \mathbb{R}^N : Av=0\} \neq \{0\}\), dann ist die Lösung nicht eindeutig. In diesem Fall ist das Problem nicht sachgemäß gestellt!
Lösungsverfahren:
Falls \(A^\top A\) invertierbar und gut konditioniert:
- Berechne \(QR\)-Zerlegung von \(A=QR\)
- \(R[1:N, 1:N]x = (Q^T b)[1:N]\) lösen
- Berechne \(A^\top A\) mittels Cholesky-Zerlegung
3.2 Singularitätszerlegung
(en. Singular Value Decomposition, SVD)
Falls \(A\) singulär (nicht invertierbar) oder schlecht konditioniert ist, kann die Singularitätszerlegung verwendet werden um die Normalengleichung zu lösen.
Sei \(A \in \mathbb{R}^{K \times N}\) mit \(R = \text{rang} A\). Dann existieren Singulärwerte \(\sigma_1, \dots, \sigma_R > 0\), Matrizen \(V \in \mathbb{R}^{K \times R}, U \in \mathbb{R}^{N \times R}\) mit \(V^\top V = I_R = U^\top U\) und eine Zerlegung
\[ A = U \Sigma V^\top \]
mit \(\Sigma = \text{diag}(\sigma_1, \dots, \sigma_R)= \text{diag}(\sqrt{\lambda_1}, \dots, \sqrt{\lambda_R}) \in \mathbb{R}^{R \times R}\).
Wir führen eine Drehung mit \(V^\top\) durch, skalieren mit den Singularitätswerten und passen die Dimension von \(N\) auf \(K\) mit \(\Sigma\) und drehen dann in \(\mathbb{R}^K\) mittels \(U\).
- \(\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_R > 0\) sind die geordneten Eigenwerte von \(A^\top A\) und \(A A^\top\).
- \(U\) enthält die Eigenvektoren von \(A A^\top\) als Spaltenvektoren und \(V\) die Eigenvektoren von \(A^\top A\) als Zeilenvektoren (transponiert).
- Es gilt \(A = \sum_{k=1}^R \sigma_k v_k u_k^\top\).
- Die Spalten \(v_1,...,v_R\) von V bilden eine ONB von \(\operatorname{Bild} A\).
- Die Spalten \(u_1,...,u_R\) von U lassen sich mit \(u_{R+1},...,u_N\) zu einer ONB von \(\mathbb{R}^N\) ergänzen. Dann ist \(u_{R+1},...,u_N\) eine ONB von \(\operatorname{Kern} A\).
- \(A^+ := \sum_{k=1}^R \frac{1}{\sigma_k} u_k v_k^T = U \Sigma^{-1} V^T\) heißt Pseudo-Inverse.
- Ist \(A\) regulär, so gilt \(A^{-1} = A^+\). Die Normalengleichung und damit \(\|Ax-b\|_2 = \min!\) wird durch \(x = A^+b\) gelöst.
- \(\|A\|_2 = \sigma_1\) und \(\kappa_2(A) = \frac{\sigma_1}{\sigma_R}\). Berechnung aufwendig, erfordert Eigenwertberechnung von \(S = A^T A\). Stabile Berechnung mit \(O(N^3)\) Operationen!
- Die Berechnung von \(A^+b\) ist für \(\sigma_R \ll 1\) numerisch nicht stabil.
3.3 Tikhonov-Regularisierung
Stabilisiert die Singularitätszerlegung, indem ein Regularisierungsterm hinzugefügt wird.
Zu \(A \in \mathbb{R}^{M \times N}, b \in \mathbb{R}^M, \alpha > 0\) existiert genau ein \(x^\alpha \in \mathbb{R}^N\), das die Tikhonov-Regularisierung \[ F_\alpha(x) := \frac{1}{2}|Ax-b|_2^2 + \frac{\alpha}{2}|x|_2^2 \] minimiert. Es gilt \[ x^\alpha = (A^\top A + \alpha I_N)^{-1} A^\top b. \]