Kurs:Numerik I/Lösung mit QS-Zerlegung

Einleitung

Diese Lernressoure zum Thema QS-Zerlegung kann als Wiki2Reveal Folien angezeigt werden. Einzelne Abschnitte werden als Folien betrachtet und Folien können in der Wiki2Reveal-Darstellung handschriftlich im Browser beschriftet werden.

Lösung mittels QS-Zerlegung

Für eine QS-Zerlegung stellen wir zunächst fest, dass sich eine Matrix mit maximalen Spaltenrang in ein Produkt von

einer orthogonalen quadratischen Matrix $Q\in \mathbb {R} ^{n\times n}$ und
erweiterten oberen Dreiecksmatrix $S\in \mathbb {R} ^{n\times k}$

Lemma - QS-Zerlegung

Sei $A\in \mathbb {R} ^{n\times k}$ mit $n\geq k$ und $\operatorname {Rang} (A)=k$ gegeben und $A=QS$ eine Zerlegung von $A$ , wobei $Q\in \mathbb {R} ^{n\times n}$ eine orthogonale und $S\in \mathbb {R} ^{n\times k}$ eine Matrix der folgenden Gestalt ist mit $\mathbf {0} :=(0)\in \mathbb {R} ^{(n-k)\times k}$ :

S={\begin{pmatrix}R\\\mathbf {0} \end{pmatrix}}\in \mathbb {R} ^{n\times k},R:={\begin{pmatrix}*&\ldots &*\\\vdots &\ddots &\vdots \\0&\ldots &*\end{pmatrix}}\in \mathbb {R} ^{k\times k},

Dann gilt für die obere Dreiecksmatrix $R$ :

(i)

\det(R)\neq 0

,

(ii)

\operatorname {cond} _{2}(A^{T}A)={\big (}\operatorname {cond} _{2}(R){\big )}^{2}

.

Beweis - QS-Zerlegungslemma

Der Beweis gliedert sich in folgende Schritte:

Beweisschritt 1 - Einsetzung

Da $Q$ eine orthogonale Matrix und $A=QS$ ist, gilt mit $(QS)^{T}=S^{T}Q^{T}$ :

A^{T}A=\underbrace {S^{T}Q^{T}} _{=A^{T}}\,\underbrace {QS} _{=A}=S^{T}S={\begin{pmatrix}R^{T}&\mathbf {0} \end{pmatrix}}{\begin{pmatrix}R\\\mathbf {0} \end{pmatrix}}=R^{T}R.

Dabei wurde u.a. verwendet, dass für orthogonale Matrizen $Q$ und $s\in \mathbb {R} ^{n}$ gilt:

(Qs)^{T}Qs=\langle Qs,Qs\rangle =\langle s,s\rangle =s^{T}s

Beweisschritt 2 - Betrachtung des Ranges

Wegen $\operatorname {Rang} (A)=k$ ist nach Lemma 4.1 insbesondere

0\neq \det(A^{T}A)=\det(R^{T}R)={\big (}\det(R){\big )}^{2}

.

Damit gilt $det(A^{T}A)>0$ und $det(R)\not =0$ und (i) ist korrekt.

Beweisschritt 3 - Kondition der Matrix

Ferner erhält man durch die Gleichheit von $A^{T}A=R^{T}R$ aus Beweisschritt 1 die Gleichheit der Kondition

\operatorname {cond} _{2}(A^{T}A)=\operatorname {cond} _{2}(R^{T}R).

Damit kann man die Konditionszahl von $A^{T}A$ über die Konditionszahl der invertierbaren oberen Dreiecksmatrix $R$ berechnen.

Beweisschritt 4 - Kondition der Matrix

Wendet man nun Lemma für die Konditionszahl auf die reguläre obere Dreiecksmatrix $R$ an, erhält man:

\operatorname {cond} _{2}(R^{T}R)=(\operatorname {cond} _{2}(R))^{2}.

Beweisschritt 5 - Kondition der Matrix

Nach dem Lemma über Eigenwerte positiv definiter Matrizen und mit der Invertierbarkeit von $R$ sind die Eigenwerte der Matrix $R^{T}R$ positiv. Insgesamt erhält man mit Schritt 4 die Eigenschaft (ii) aus dem Lemma.

\operatorname {cond} _{2}(A^{T}A)=\operatorname {cond} _{2}(R^{T}R)={\big (}\operatorname {cond} _{2}(R){\big )}^{2}

q.e.d.

Bemerkung - Fehlerquadratprobleme

In dem obige Satz für die QS-Zerlegung wurde der maximale Spaltenrang $k$ der Matrix $A\in \mathbb {R} ^{n\times k}$ mit $k\leq n$ vorausgesetzt. Der folgende Satz gibt an, wie Fehlerquadratprobleme mittels einer QS-Zerlegung von $A$ gelöst werden können.

Satz - QS-Zerlegung

Voraussetzungen Es seien $A\in \mathbb {R} ^{n\times k}$ mit $n\geq k$ , $\operatorname {Rang} (A)=k$ , $R\in \mathbb {R} ^{k\times k}$ eine obere Dreiecksmatrix und $Q\in \mathbb {R} ^{n\times n}$ eine orthogonale Matrix und $S\in \mathbb {R} ^{n\times k}$ gegeben wie in Lemma - QS-Zerlegung. Weiter sei der Vektor $Q^{T}b\in \mathbb {R} ^{n}$ dargestellt in der Form

Q^{T}b={\begin{pmatrix}y^{(1)}\\y^{(2)}\end{pmatrix}}\in \mathbb {R} ^{n},y^{(1)}\in \mathbb {R} ^{k},y^{(2)}\in \mathbb {R} ^{n-k}.

Folgerung QS-Zerlegung

Dann ist der Vektor $x^{*}\in \mathbb {R} ^{k}$ genau dann die eindeutige Lösung des linearen Ausgleichsproblems $\min _{x\in \mathbb {R} ^{k}}\|Ax-b\|_{2},$ , wenn $x^{*}$ eindeutige die Lösung des Systems $Rx=y^{(1)}$ ist. Außerdem gilt für den Fehler $\|Ax^{*}-b\|_{2}=\left\|y^{(2)}\right\|_{2}.$

Bemerkung QS-Zerlegung

Ziel des Lösungsverfahrens für ein gesuchtest $x$ ist es, sich möglichst gut bzgl. der euklidschen Norm mit $Ax$ an den Vektor $b$ anzunähern. Die eindeutige Lösung des linearen Ausgleichsproblems liefert mit $y^{(1)}$ und $y^{(2)}$ :

$Rx=y^{(1)}$ ein Gleichungssystem mit einer einfach zu lösenden oberen Dreiecksmatrix $R$ und
mit $y^{(2)}$ bzw. $\|y^{(2)}\|_{2}$ ein Maß für die minimale Abweichung von $Ax$ und $b$ über $\|Ax^{*}-b\|_{2}=\|y^{(2)}\|_{2}$ bzgl. der Lösung $x^{*}$ .

Beweis

Für einen beliebigen Vektor $v\in \mathbb {R} ^{k}$ erhalten wir für eine orthogonale Matrix $Q$ eine längetreue Abbildung über:

\|Qv\|_{2}^{2}=\langle Qv,Qv\rangle =\langle v,v\rangle =\|v\|_{2}^{2}\quad (\ast )

Ferner gilt für orthogonale Matrizen $Q^{T}Q=I$ , wobei $I$ die Einheitsmatrix als neutrales Element der Matrixmultiplikation ist.

Beweisschritt 1 - Anwendung QS-Zerlegung

unter Verwendung des QS-Zerlegungssatzes und der Ersetzung von $A=QS$ erhält man:

{\begin{array}{rcl}\|b-Ax^{*}\|_{2}^{2}&=&\left\|(QQ^{T})b-QSx\right\|_{2}^{2}={\bigg \|}Q{\big (}\underbrace {Q^{T}b-Sx} _{=v}{\big )}{\bigg \|}_{2}^{2}\\&=&{\big \|}\underbrace {Q^{T}b-Sx} _{=v}{\big \|}_{2}^{2}{\mbox{ mit QS-Zerlegungssatz}}\\&=&\left\|{\begin{pmatrix}y^{1}\\y^{2}\end{pmatrix}}-{\begin{pmatrix}R\\\mathbf {0} \end{pmatrix}}x\right\|_{2}^{2}=\left\|y^{(1)}-Rx\right\|_{2}^{2}+\left\|y^{(2)}\right\|_{2}^{2}\\\end{array}}

Beweisschritt 2 - Anwendung QS-Zerlegung

Daraus folgt für einen beliebigen Vektor ${\widehat {x}}\in \mathbb {R} ^{k}$

\|A{\widehat {x}}-b\|_{2}\geq \min _{x\in \mathbb {R} ^{k}}\|Ax-b\|_{2}=\left\|y^{(2)}\right\|_{2}

Beweisschritt 3 - Eindeutige Lösung QS-Zerlegung

Ferner gilt

\|b-Ax\|_{2}=\left\|y^{(2)}\right\|_{2}\Leftrightarrow \left\|y^{(1)}-Rx\right\|_{2}=0\Leftrightarrow Rx=y^{(1)}.

Damit ist alles gezeigt.

q.e.d.

Bemerkung

Nach dem Satz und dem Lemma zur QS-Zerlegung besitzen das $l_{2}$ -Approximationsproblem, wobei man versucht $Ax$ in der $l_{2}$ -Norm möglichst gut anzunähern, eine eindeutige Lösung. Dies ist äquivalent zu der eindeutigen Lösbarkeit des Systems $Rx=y^{1}$ . Abschließend geben wir zu dem letzten Satz noch ein Beispiel.

Beispiel - QS-Zerlegung

Wir betrachten das Fehlerquadratproblem mit den Daten

A:={\begin{pmatrix}1&0\\1&3\\1&4\\1&7\end{pmatrix}},\quad b:={\begin{pmatrix}1\\2\\6\\4\end{pmatrix}}.

Beispielschritt 1 - gesuchter Vektor

Wir suchen nun eine Vektor $x^{*}\in \mathbb {R} ^{2}$ , der den Abstand $\|Ax-b\|_{2}$ zwischen $Ax$ und $b$ minimiert.

Beispielschritt 2 - QS-Zerlegung

Der Spaltenrang von $A$ ist 2, da die beiden Spalten von $A$ nicht linear abhängig sind. Nun liefert die QS-Zerlegung von $A$ mit gleichzeitiger Berechnung von $Q^{T}b$ die folgende obere Dreiecksmatrix $R$ und die folgenden Vektoren $y^{(1)}$ und $y^{(2)}$ :

R:={\begin{pmatrix}-2&-7\\0&-5\end{pmatrix}},\quad y^{(1)}:={\begin{pmatrix}-{\frac {13}{2}}\\-{\frac {5}{2}}\end{pmatrix}},\quad y^{(2)}:={\begin{pmatrix}{\frac {99}{34}}\\-{\frac {5}{34}}\end{pmatrix}}.

Beispielschritt 3 - Lösung für das Ausgleichsproblem

Die Lösung von $Rx=y^{(1)}$ bzw.

{\begin{pmatrix}-2&-7\\0&-5\end{pmatrix}}{\begin{pmatrix}x_{1}\\x_{2}\end{pmatrix}}={\begin{pmatrix}-{\frac {13}{2}}\\-{\frac {5}{2}}\end{pmatrix}}

lautet

x^{*}=(x_{1}^{*},x_{2}^{*})^{T}:=\left({\frac {3}{2}},{\frac {1}{2}}\right)^{T}.

Beispielschritt 4 - Approximationsfehler des Ausgleichsproblems

Der zugehörige Approximationsfehler in der Euklidischen Norm ist gegeben durch

\|b-Ax^{*}\|_{2}=\left\|y^{2}\right\|_{2}=\left\|{\begin{pmatrix}{\frac {99}{34}}\\-{\frac {5}{34}}\end{pmatrix}}\right\|_{2}={\frac {\sqrt {9826}}{34}}={\frac {17{\sqrt {34}}}{34}}={\frac {1}{2}}{\sqrt {34}}\approx 2.915.

Vergleich der Lösungswege bzgl. Zerlegung

Wir wollen nun die beiden beschriebenen Lösungswege für das $l_{2}$ -Problem, die Lösung der Normalgleichungen mittels Cholesky-Zerlegung und die Lösung über den in QS-Zerlegungssatz dargestellten Weg, vergleichen.

Gemeinsamkeiten

In beiden Fällen muss die rechte Seite $b$ des Systems $Ax=b$ von links mit einer Matrix multipliziert werden. Weiter sind bei der Cholesky-Zerlegung zwei und bei dem Weg über die QS-Zerlegung ein gestaffeltes lineares Gleichungssystem zu lösen. Da die Lösung eines solchen Systems nur etwa ${\frac {k^{2}}{2}}$ Multiplikationen und Divisionen erfordert, vernachlässigen wir hier diesen zusätzlichen Aufwand bei der Cholesky-Zerlegung.

Daten für das Gleichungssystem

Der Lösungsvektor $x\in \mathbb {R} ^{k}$ hat in der Regel eine feste Dimension $k\leq n$ während $n$ die Anzahl der Gleichungen für Daten angibt, die bzgl. $\|Ax-b\|_{2}$ in der euklidischen Norm minimiert werden soll.

Berechnungsaufwand - Cholesky-Zerlegung

Im Fall der Lösung der Normalgleichungen benötigt man ferner zur Berechnung der symmetrischen Matrix $A^{T}A$ etwa ${\frac {1}{2}}nk^{2}$ und für die Cholesky-Zerlegung etwa ${\frac {1}{6}}k^{3}$ wesentliche Rechenoperationen. Daneben erfordert die Berechnung des Residuenvektors $Ax-b$ weitere $n\cdot k$ Multiplikationen, so dass die zuletzt genannten Teilaufgaben zusammen ungefähr

{\frac {1}{2}}n\cdot k^{2}+{\frac {1}{6}}\cdot k^{3}+n\cdot k

wesentliche Rechenoperationen erfordern, das sind

\sim {\frac {1}{2}}nk^{2}

für

n\gg k

und

\sim {\frac {2}{3}}n^{3}

für

n\approx k

.

Berechnungsaufwand - QS-Zerlegung

Für das sich aus dem QS-Zerlegungssatz ergebende Vorgehen benötigt man über die Matrix-Vektor-Multiplikation und die Lösung eines Dreieckssystems hinaus nur die Berechnung der QS-Zerlegung von $A$ .

Vergleich n nahe bei k

Wenn $n$ nahe bei $k$ liegt ( $n\approx k$ ) benötigen im Vergleich der beiden Verfahren demnach für beide Wege etwa gleich viele wesentliche Rechenoperationen.

Vergleich n groß im Vergleich zu k

Ist $n$ sehr groß im Vergleich zu $k$ müssen über den Weg der QS-Zerlegung etwa doppelt so viele wesentliche Rechenoperationen durchgeführt werden, wie der Lösung über die Normalgleichungen (Cholesky).

Berechnungsaufwand und Konditionszahl

Letzterem Nachteil der QS-Zerlegung steht jedoch entgegen, dass bei ihr nach Lemma 4.6 die Konditionszahl $\kappa :=\left[\operatorname {cond} _{2}\left(A^{T}A\right)\right]^{1/2}$ bzw. für quadratisches $A$ (vgl. Satz 4.5) die Zahl $\kappa :=\operatorname {cond} _{2}(A)$ den Rundungsfehlereinfluss bei der Lösung des Dreieckssystems bestimmt, während dies bei dem Weg über die Normalgleichungen die Zahl $\kappa ^{2}$ ist. Wir geben dazu ein (allerdings etwas extremes) Beispiel:

Beispiel - Vergleich Konditionszahl

Es sei $n:=6,k:=5$ und $A\in \mathbb {R} ^{6\times 5}$ mit einem $\varepsilon >0$ gegeben durch

A:={\begin{pmatrix}1&1&1&1&1\\\varepsilon &&&&\\&\varepsilon &&&\\&&\varepsilon &&\\&&&\varepsilon &\\&&&&\varepsilon \end{pmatrix}}.

Berechnung von A^TA

Damit ergibt sich für die Berechnung von $A^{T}A$ :

A^{T}A={\begin{pmatrix}1+\varepsilon ^{2}&1&1&1&1\\1&1+\varepsilon ^{2}&1&1&1\\1&1&1+\varepsilon ^{2}&1&1\\1&1&1&1+\varepsilon ^{2}&1\\1&1&1&1&1+\varepsilon ^{2}\end{pmatrix}}.

Maschinengenauigkeit 1

Es seien nun $b:=10$ und $r:=10$ die Basis und Mantissenlänge des verwendeten Computers, so dass $eps=5\cdot 10^{-10}$ die zugehörige Maschinengenauigkeit ist. Weiter sei $\varepsilon :=0.5_{10}-5$ und damit $\varepsilon ^{2}=.25_{10}-10$ . Damit liegt $\varepsilon ^{2}$ unterhalb der zugehörigen Maschinengenauigkeit und wird als 0 gespeichert.