Kurs:Numerik I/Orthonormalisierungsverfahren

Einleitung

Diese Seite kann als Wiki2Reveal Folien angezeigt werden. Einzelne Abschnitte werden als Folien betrachtet und Änderungen an den Folien wirken sich sofort auf den Inhalt der Folien aus.

Einführung - Orthogonalisierungsverfahren

In diesem Abschnitt soll für eine gegebene Matrix $A\in \mathbb {R} ^{n\times k}$ mit $1\leq k\leq n$ eine Zerlegung der Form

$A=QS$

bestimmt werden, wobei $Q\in \mathbb {R} ^{n\times n}$ eine orthogonale Matrix ist und $S\in \mathbb {R} ^{n\times n}$ eine um Nullzeilen ergänzte obere Dreiecksmatrix ist,

Orthogonale Matrix

Eine orthogonole eine reguläre (invertierbare) Matrix mit der Eigenschaft

$Q^{-1}=Q^{T}$

Wenn man nun das Matrixprodukt $Q^{-1}\cdot Q=Q^{T}\cdot Q=E_{n}$ berechnet, dann bilden die Spaltenvektoren $(q_{1},\ldots ,q_{n})$ Orthonormalbasis vom $\mathbb {R} ^{n}$ mit:

$Q^{T}\cdot Q={\begin{pmatrix}\langle q_{1},q_{1}\rangle &\langle q_{1},q_{2}\rangle &\cdots &\langle q_{1},q_{n}\rangle \\\langle q_{2},q_{1}\rangle &\langle q_{2},q_{2}\rangle &\cdots &\langle q_{2},q_{n}\rangle \\\vdots &\ddots &\cdots &\vdots \\\langle q_{n},q_{1}\rangle &\langle q_{n},q_{2}\rangle &\cdots &\langle q_{n},q_{n}\rangle \\\end{pmatrix}}$.

Orthonormal Vektoren und Einheitsmatrix

• Die Orthogonalität der Spaltenvektoren $(q_{1},\ldots ,q_{n})$ führt zu der Eigenschaft $\langle q_{j},q_{k}\rangle =0$ für $j\not =k$.
• Die Normalisierung der Basisvektoren $q_{k}\in \mathbb {R} ^{n}$ führt zu der Eigenschaft $\|q_{k}\|=1$ und damit auch $\|q_{k}\|^{2}=\langle q_{k},q_{k}\rangle =1$, was den Diagonalelementen der Matrix $Q^{T}\cdot Q=E_{n}$ entspricht.

Um Nullzeilen ergänzte obere Dreieckmatrix

In der Zerlegung $A=QS$ ist $S$ mit einer oberen Dreiecksmatrix $R$, die wie folgt mit Nullzeilen ergänzt wird:

$S={\begin{pmatrix}R\\\mathbf {0} \end{pmatrix}}\in \mathbb {R} ^{n\times k},\quad R:={\begin{pmatrix}\ast &...&*\\&\ddots &\vdots \\0&&*\end{pmatrix}}\in \mathbb {R} ^{k\times k},\quad \mathbf {0} :=(0)\in \mathbb {R} ^{(n-k)\times k}.$

Im Fall $n=k$ ist demnach insbesondere $S=R$.

Ziel der QR/QS Zerlegung

Wie wir zeigen wollen, ermöglicht eine solche QR- bzw. QS-Zerlegung von $A$ sowohl die stabile Lösung linearer Gleichungssysteme $Ax=b$ als auch die stabile Lösung von Ausgleichsproblemen

$\min _{x\in \mathbb {R} ^{k}}\|b-Ax\|_{2}.$

Dafür benötigen wir einige Eigenschaften orthogonaler Matrizen.

Elementare Eigenschaften orthogonaler Matrizen

In der Geometrie der Ebene kennt man die Kongruenzabbildungen, bei denen die Bilder von Strecken wieder Strecken der gleichen Länge abgebildet. Kongruenzabbildungen, wie Geradenspiegelungen sind dabei längentreu und winkeltreu, d.h. unter der Abbildung bleiben Winkel zwischen Strecken und die Länge von Strecken unverändert (invariant). Orthogonale Matrizen besitzen als Abbildung die Eigenschaft der Längeninvarianz.

Lemma - Längeninvarianz/Isometrie für orthogonale Matrizen

Sei $Q\in \mathbb {R} ^{n\times n}$ eine orthogonale Matrix. Dann ist auch $Q^{T}$ eine orthogonale Matrix und es gilt

$\|Qx\|_{2}=\|x\|_{2}=\left\|Q^{T}x\right\|_{2},\quad x\in \mathbb {R} ^{n}.$

Bemerkung

Eine orthogonal Matrix kann man als Abbildung von $\mathbb {R} ^{n}$ nach $\mathbb {R} ^{n}$ auffassen. Eine orthogonal Matrix verändert dabei die Länge des Vektors nicht in der Euklidischen Norm $\|\cdot \|_{2}$. Der Begriff der Isometrie ist dabei allgemein für metrische Räumen definiert. Ein zentrische Streckung $Z_{\lambda }$ mit einem Streckungsfaktor $\lambda \notin \{+1,-1\}$ ist z.B. keine Isometrie, da unter der Abbildung der Bildvektor $Z_{\lambda }(x)$ die $|\lambda |$-fache Länge von $x$ besitzt.

Beweis

Der Beweis erfolgt in zwei Schritten:

• Man zeigt, dass auch die transponierte Matrix eine orthogonale Matrix ist.
• Man zeigt die Längeninvarianz unter Darstellung des Skalarproduktes als Matrixmultiplikation (d.h. $\langle x,y\rangle =x^{T}\cdot y$, wobei $x,y\in \mathbb {R} ^{n\times 1}$ als Spaltenvektoren aufgefasst werden und $x^{T}\cdot y$ als Matrizenprodukt einer $1\times n$-Matrix $x^{T}$ und einer $n\times 1$-Matrix $y$ ist.

Beweis 1 - inverse Matrix - orthogonale Matrix

Wegen

$(Q^{T})^{-1}=(Q^{-1})^{-1}=Q=(Q^{T})^{T}$

ist auch $Q^{T}$ eine orthogonale Matrix.

Beweis 2 - Umformungen in einer Matrixmultiplikation

Des Weiteren hat man für $Q$

$\|Qx\|_{2}^{2}=(Qx)^{T}Qx=x^{T}Q^{T}Qx=x^{T}Q^{-1}Qx=x^{T}x=\|x\|_{2}^{2}$

und somit auch $\left\|Q^{T}x\right\|_{2}=\|x\|_{2}$ für die orthogonale Matrix $Q^{T}$.

q.e.d.

Bermerkung

Die Eigenschaft der Längentreue der durch $Q$ und $Q^{T}$ definierten linearen Abbildungen bezeichnet man auch als isometrisch bezüglich der Euklidischen Norm, bzw. als Isometrie.

Lemma - Euklidische Norm und Determinante orthogonaler Matrizen

Für die Operatornorm $\|\cdot \|_{2}:\mathbb {R} ^{n\times n}\to \mathbb {R} _{0}^{+}$ einer linearen Abbildung und für die Determinante orthogonaler Matrizen gilt:

• (i) $\|Q\|_{2}=\left\|Q^{T}\right\|_{2}=1$,
• (ii) $|\det(Q)|=\left|\det(Q^{T})\right|=1$.

Beweis - Euklidische Norm und Determinante orthogonaler Matrizen

Der Beweis gliedert sich in zwei Teilaussagen

• (i) zur Operatornorm für lineare Abbildung und
• (ii) zur Determinate von orthogonalen Matrizen.

Beweis (i) Längeninvarianz/Isometrie

Wenn man das Lemma Längeninvarianz bzw. Isometrie für orthogonale Matrizen auf $Q$ und $Q^{T}$ anwendet, erhält man die Aussage (i) über die Definition der Operatornorm mit:

${\begin{array}{rcl}\|Q\|_{2}&:=&\max _{\|x\|_{2}=1}\|Qx\|_{2}\\&{\stackrel {Isometrie\,Q}{=}}&\max _{\|x\|_{2}=1}\underbrace {\|x\|_{2}} _{=1}\\&{\stackrel {Isometrie\,Q^{T}}{=}}&\max _{\|x\|_{2}=1}\left\|Q^{T}x\right\|_{2}=1\\\end{array}}$

Beweis (ii) Längeninvarianz/Isometrie

Mit der Einheitsmatrix $E_{n}\in \mathbb {R} ^{n\times n}$ und der Multiplikativität der Determinanten $\det(A\cdot B)=\det(A)\cdot \det(B)$ folgt die Aussage (ii) über die Eigenschaft, dass die transponierte Matrix $Q^{T}$ multiplikativ invers zu der Matrix $Q$ ist.

${\begin{array}{rcl}1&=&\det(E_{n})=\det(Q\cdot Q^{T})=\det(Q)\cdot \det(Q^{T})\\&=&[\det(Q)]^{2}=\left[\det(Q^{T})\right]^{2}.\\\end{array}}$

Dabei geht in die Gleichungskette ein, dass die Determinante der transponierten Matrix und der Matrix selbst gleich sind. q.e.d.

Korollar - Konditionszahl orthogonaler Matrizen

Seien $A\in \mathbb {R} ^{n\times n}$ eine reguläre und $Q\in \mathbb {R} ^{n\times n}$ eine orthogonale Matrix. Dann gilt für die Konditionszahl

$\operatorname {cond} _{2}(QA)=\operatorname {cond} _{2}(A).$

Beweis

Der Beweis gliedert sich in die folgenden Teile:

• Verwendung der Längeninvarianz/Isometrie für orthogonale Abbildungen

Beweis 1 - Längeninvarianz orthogonaler Abbildungen

Nach Lemma zur Längeninvarianz/Isometrie für orthogonale Abbildungen gilt $\|Qy\|_{2}=\|Ay\|_{2}$ für alle $y\in \mathbb {R} ^{n}$. Mit $y:=Ax\in \mathbb {R} ^{n}$ erhält man

$\|(Q\cdot A)\cdot x\|_{2}=\|Q\cdot (A\cdot x)\|_{2}=\|Ax\|_{2}$

für $x\in \mathbb {R} ^{n}$ und somit

Beweis 2 - Längeninvarianz orthogonaler Abbildungen

$\|A\|_{2}=\max _{\|x\|_{2}=1}\|Ax\|_{2}=\max _{\|x\|_{2}=1}\|QAx\|2=\|QA\|_{2}.$

Weiter gilt mit Lemma 3.24 (i)

$\left\|A^{-1}\right\|_{2}=\left\|A^{-1}Q^{T}Q\right\|_{2}\leq \left\|A^{-1}Q^{T}\right\|_{2}\|Q\|_{2}=\left\|A^{-1}Q^{T}\right\|_{2}\leq \left\|A^{-1}\right\|_{2}\left\|Q^{T}\right\|_{2}=\left\|A^{-1}\right\|_{2}$

und folglich $\|A^{-1}\|_{2}=\left\|A^{-1}Q^{T}\right\|_{2}$. Damit ergibt sich

$\operatorname {cond} _{2}(QA)=\|QA\|_{2}\left\|A^{-1}Q^{-1}\right\|_{2}=\|QA\|_{2}\left\|A^{-1}Q^{T}\right\|_{2}=\|A\|_{2}\left\|A^{-1}\right\|_{2}=\operatorname {cond} _{2}(A).$

q.e.d.

Multiplikation einer quadratischen Matrix $A$ von links mit einer orthogonalen Matrix führt also auf eine Matrix mit derselben Kondition wie $A$. Weiter gilt:

Lemma - Produkt orthogonaler Matrizen

Seien $Q_{1},Q_{2}\in \mathbb {R} ^{n\times n}$ orthogonale Matrizen. Dann ist auch $Q_{1}Q_{2}$ eine orthogonale Matrix.

Beweis.

Die Aussage des Lemma folgt aus

$(Q_{1}Q_{2})^{-1}=Q_{2}^{-1}Q_{1}^{-1}=Q_{2}^{T}Q_{1}^{T}=(Q_{1}Q_{2})^{T}.$

q.e.d.

3.4.2 Lösung linearer Gleichungssysteme mittels QR-Zerlegung

Wir betrachten nun wieder die Lösung eines linearen Gleichungssystems $Ax=b$ für eine reguläre Matrix $A\in \mathbb {R} ^{n\times n}$ und beliebige rechte Seite $b\in \mathbb {R} ^{n}$ und gehen davon aus, dass eine Zerlegung der Form $A=QR$ von $A$ mit einer orthogonalen Matrix $Q\in \mathbb {R} ^{n\times n}$ und einer oberen Dreiecksmatrix $R\in \mathbb {R} ^{n\times n}$ gegeben ist (vgl. (3.37) - (3.39)). Wegen Lemma 3.24 gilt offenbar

$\det(A)\neq 0\Leftrightarrow \det(R)\neq 0.$

Weiter hat man die Äquivalenzen

$Ax=b\Leftrightarrow QRx=b\Leftrightarrow Rx=Q^{T}b$

und mit Korollar 3.25 die Beziehungen

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

D. h., das Gleichungssystem $Ax=b$ ist äquivalent zu dem gestaffelten Gleichungssystem $Rx=Q^{T}b$, wobei die Matrix $R$ bezüglich der Spektralnorm keine schlechtere Kondition als $A$ aufweist und gemäß 3.40 $\left\|Q^{T}b\right\|_{2}=\|b\|_{2}$ gilt.

3.4.3 QR-Zerlegung mittels Gram-Schmidt-Orthogonalisierung

Es sei nun $A\in \mathbb {R} ^{n\times n}$ eine quadratische Matrix mit $\det(A)\neq 0$ und es seien eine orthogonale Matrix $Q\in \mathbb {R} ^{n\times n}$ und eine obere Dreiecksmatrix $R\in \mathbb {R} ^{n\times n}$ gesucht, so dass

(3.41) $A=QR$

gilt. Schreiben wir $A,Q$ und $R$ mit Vektoren $a^{j},q^{j}\in \mathbb {R} ^{n}$ in der Form

$A={\begin{pmatrix}a^{1}&...&a^{n}\end{pmatrix}},\quad Q={\begin{pmatrix}q^{1}&...&q^{n}\end{pmatrix}},\quad R={\begin{pmatrix}r_{11}&...&r_{1n}\\&\ddots &\vdots \\&&r_{nn}\end{pmatrix}},$

so ist die Gleichung (3.41) offenbar äquivalent mit den Gleichungen

(3.42) $a^{j}=\sum _{i=1}^{j}r_{ij}q^{i},\quad j=1,...,n,$

wobei $r_{ij}=0$ für $i>j$ berücksichtigt wurde. Dabei sind die $a^{j}$ $(j=1,...,n)$ linear unabhängig und die $q^{j}$ $(j=1,...,n)$ wegen der Orthogonalität von $Q$ paarweise orthonormal und damit auch linear unabhängig. Die Gleichungen (3.42) implizieren somit für $1\leq s\leq n$

(3.43) $\operatorname {span} \{a^{1},...,a^{s}\}=\operatorname {span} \{q^{1},...,q^{s}\}.$

Folglich suchen wir eine orthogonale Basis $\{q^{1},...,q^{n}\}$ des $\mathbb {R} ^{n}$, für welche die Gleichungen (3.42) erfüllt sind. Wir wollen im Folgenden zeigen, dass man eine solche durch Orthogonalisierung der $a^{j}$ $(j=1,...,n)$ mittels des aus der Linearen Algebra bekannten Gram-Schmidt-Orthonormierungsverfahrens erhält.

Sind $a^{j}$ $(j=1,...,n)$ die gegebenen Spalten von $A$, welche nach Voraussetzung eines Basis des $\mathbb {R} ^{n}$ bilden, so lautet für diese das Gram-Schmidt-Orthonormierungsverfahren, wobei wir die durch das Verfahren erzeugten Vektoren gleich mit $q^{j}$ bezeichnen:

(3.44) ${\hat {q}}^{j}:=aj-\sum _{i=1}^{j-1}\left[(a^{j})^{T}q^{i}\right]q^{i},\quad q^{j}:={\frac {{\hat {q}}^{j}}{\|{\hat {q}}^{j}\|_{2}}},\quad j=1,2,...,n.$

Bekanntlich gilt für die so erzeugten Vektoren $q^{j}$ $(j=1,...,n)$ die Identität (3.43) und sind diese paarweise orthonormal.

Aus den Gleichungen (3.44) leitet man unmittelbar die folgenden Gleichungen her:

$a^{j}=\left\|{\hat {q}}^{j}\right\|_{2}q^{j}+\sum _{i=1}^{j-1}\left[(a^{j})^{T}q^{i}\right]q^{i},\quad j=1,2,...,n.$

Wie ein Vergleich mit den Gleichungen (3.42) zeigt, erhält man demnach die gewünschte QR-Faktorisierung von $A$ für die Matrix $Q$ mit den durch das Gram-Schmidt-Verfahren erzeugten Vektoren $q^{j}$ $(j=1,...,n)$ und die obere Dreiecksmatrix $R$, welche folgende Elemente hat:

$r_{ij}:=(a^{j})^{T}q^{i},\quad i=1,...,j-1,\quad r_{jj}:=\left\|{\hat {q}}^{j}\right\|_{2}\quad (j=1,...,n).$

Der hier beschriebene Orthogonalisierungsprozess ist aber unter Umständen nicht gutartig, etwa für $\|{\hat {q}}^{j}\|_{2}\ll 1$ und damit $r_{jj}\ll 1$ (s. Beispiel 3.4 in Band 1 von Quarteroni/Sacco/Saleri und S. 177 bei Stoer). Mit wachsendem $j$ kann die Orthogonalität immer stärker verloren gehen. Deshalb werden zumeist andere Methoden, wie die im folgenden Abschnitt dargestellte, zur QR-Faktorisierung von $A$ herangezogen.

3.4.4 QS-Zerlegung mittels Householder-Transformationen

Sei nun allgemeiner $A\in \mathbb {R} ^{n\times k}$ mit $1\leq k\leq n$ gegeben. Gegenstand dieses Abschnitts ist die Bestimmung einer Zerlegung der Form $A=QS$ mit einer orthogonalen Matrix $Q\in \mathbb {R} ^{n\times n}$ und einer im Fall $n>k$ nichtquadratischen Matrix $S\in \mathbb {R} ^{n\times k}$ wie in (3.39) mittels sog. Householder-Transformationen. In diesem Zusammenhang zeigen wir zunächst:

Lemma 3.27

Es sei $\omega \in \mathbb {R} ^{n}$ ein Vektor mit $\|\omega \|_{2}=1$ und ${\mathcal {H}}\in \mathbb {R} ^{n\times n}$ die Matrix

(3.45) ${\mathcal {H}}:=I-2\omega \omega ^{T}.$

Dann gilt

(3.46) ${\mathcal {H}}^{T}={\mathcal {H}}$ (${\mathcal {H}}$ ist symmetrisch),

(3.47) ${\mathcal {H}}^{2}=I$ (${\mathcal {H}}$ ist involutorisch),

(3.48) ${\mathcal {H}}^{T}{\mathcal {H}}=I$ (${\mathcal {H}}$ ist orthogonal).

Beweis.

Die Beziehungen (3.46) und (3.47) ergeben sich aus

${\mathcal {H}}^{T}=I-2(\omega \omega ^{T})^{T}=I-2\omega \omega ^{T}={\mathcal {H}},$

${\mathcal {H}}^{2}=(I-2\omega \omega ^{T})(I-2\omega \omega ^{T})=I-2\omega \omega ^{T}-2\omega \omega ^{T}+4\omega \underbrace {(\omega ^{T}\omega )} _{=1}\omega ^{T}=I.$

Die Identität (3.48) folgt aus (3.46) und (3.47).

q.e.d.

Eine Matrix vom Typ ${\mathcal {H}}$ (3.45) für ein $\omega \in \mathbb {R} ^{n}$ mit $\|\omega \|_{2}=1$ nennen wir Householder-Matrix und eine lineare Abbildung

$h:\mathbb {R} ^{n}\to \mathbb {R} ^{n},\quad h(x):={\mathcal {H}}x$

mit einer Householder-Matrix ${\mathcal {H}}\in \mathbb {R} ^{n\times n}$ bezeichnen wir als Householder-Transformation.

Householder-Matrizen kann man dazu verwenden, einen Block von Komponenten eines gegebenen Vektors $x$ (durch Multiplikation mit einer solchen Matrix von links) zu Null zu setzen. Will man insbesondere alle Komponenten von $x$ außer der $m$-ten Komponente Null setzen, so muss man den Vektor $\omega$ zur Konstruktion von ${\mathcal {H}}$ so wählen, wie er im folgenden Lemma angegeben wird. Dabei ist wieder mit $e^{m}\in \mathbb {R} ^{n}$ die $m$-te Spalte von $I\in \mathbb {R} ^{n\times n}$ gemeint.

Lemma 3.28

Gegeben sei $x\in \mathbb {R} ^{n}\setminus \{0\}$ mit $x\notin \operatorname {span} \{e^{m}\}$. Weiter sei

(3.49) $\omega :={\frac {x+\sigma e^{m}}{\|x+\sigma e^{m}\|_{2}}}$

mit

(3.50) $\sigma :=\pm \|x\|_{2}.$

Dann hat man $\|\omega \|_{2}=1$ und für ${\mathcal {H}}:=I-2\omega \omega ^{T}$

(3.51) ${\mathcal {H}}x=-\sigma e^{m}.$

Beweis.

Wegen $x\notin \operatorname {span} \{e^{m}\}$ ist der Nenner in (3.49) ungleich Null und damit $\omega$ in (3.49) wohldefiniert. Offenbar ist $\|\omega \|_{2}=1$. Ferner gilt

$\|x+\sigma e^{m}\|_{2}^{2}=\|x\|_{2}^{2}+2\sigma (e^{m})^{T}x+\sigma ^{2}=2\|x\|_{2}^{2}+2\sigma (e^{m})^{T}x=2(x+\sigma e^{m})^{T}x$

und damit

$2\omega ^{T}x={\frac {2(x+\sigma e^{m})^{T}x}{\|x+\sigma e^{m}\|_{2}}}=\|x+\sigma e^{m}\|_{2}.$

Zusammen mit (3.49) gelangt man zu der Darstellung

${\mathcal {H}}x=\left(I-2\omega \omega ^{T}\right)x=x-\|x+\sigma e^{m}\|_{2}\omega =x-(x+\sigma e^{m})=-\sigma e^{m}.$

Bemerkung 3.29

Zur Vermeidung von Stellenauslöschungen bei der Berechnung von (3.50) ist es offenbar sinnvoll, $\sigma :=\operatorname {sgn}(x_{m})\|x\|_{2}$ zu wählen, wobei $x_{m}$ die $m$-te Komponente von $x$ bezeichnet.

Wir geben ein Beispiel.

Beispiel 3.30

Sei $x:=(1,1,1,1)^{T}$ und $m:=3$. Dann errechnen wir $\|x\|_{2}=2$ und setzen wir somit $\sigma :=\operatorname {sgn}(x_{3})\|x\|_{2}=2$ sowie

$\omega :={\frac {x+\sigma e^{m}}{\|x+\sigma e^{m}\|_{2}}}={\frac {1}{2{\sqrt {3}}}}(1,1,3,1)^{T}.$

Es ergibt sich

$2\omega \omega ^{T}={\frac {1}{6}}{\begin{pmatrix}1\\1\\3\\1\end{pmatrix}}(1,1,3,1)={\frac {1}{6}}{\begin{pmatrix}1&1&3&1\\1&1&3&1\\3&3&9&3\\1&1&3&1\end{pmatrix}}$

und demnach

${\mathcal {H}}:=I-2\omega \omega ^{T}={\frac {1}{6}}{\begin{pmatrix}5&-1&-3&-1\\-1&5&-3&-1\\-3&-3&-3&-3\\-1&-1&-3&5\end{pmatrix}},\quad {\mathcal {H}}x={\begin{pmatrix}0\\0\\-2\\0\end{pmatrix}}.$

---

Will man für ein $k\geq 1$ die Komponenten $x_{i}$ $(i=1,...,k-1)$ von $x\in \mathbb {R} ^{n}$ unverändert lassen und gleichzeitig $x_{i}=0$ $(i=k+1,...,n)$ erreichen, bekommt man dies, indem man $x$ von links mit der Transformationsmatrix

$P_{k}:={\begin{pmatrix}I_{k-1}&0\\0&{\mathcal {H}}_{k}\end{pmatrix}}$

multipliziert. Dabei ist $I_{\ell }$ die $(\ell \times \ell )$-Einheitsmatrix und ${\mathcal {H}}_{k}$ die $(n-k+1)\times (n-k+1)$-Householder-Matrix der Form

${\mathcal {H}}_{k}:=I_{n-(k-1)}-2\omega ^{k}(\omega ^{k})^{T},\quad \omega ^{k}:={\frac {x^{n-k+1}+\sigma _{k}e^{n-k+1,1}}{\|x^{n-k+1}+\sigma _{k}e^{n-k+1,1}\|_{2}}},\quad \sigma _{k}:=\pm \left\|x^{n-k+1}\right\|_{2},$

welche mit dem aus den letzten $n-k+1$ Komponenten von $x$ bestehenden Vektor $x^{n-k+1}\in \mathbb {R} ^{n-k+1}$ und der ersten Spalte $e^{n-k+1,1}\in \mathbb {R} ^{n-k+1}$ der Einheitsmatrix $I_{n-k+1}$ gebildet wird. Denn ist $x^{k-1}\in \mathbb {R} ^{k-1}$ der aus den ersten $k-1$ Komponenten von $x$ bestehende Vektor, so hat man nach Lemma 3.28

$P_{k}x={\begin{pmatrix}I_{k}&0\\0&{\mathcal {H}}_{k}\end{pmatrix}}{\begin{pmatrix}x^{k-1}\\x^{n-k+1}\end{pmatrix}}={\begin{pmatrix}x^{k-1}\\{\mathcal {H}}_{k}x^{n-k+1}\end{pmatrix}}={\begin{pmatrix}x^{k-1}\\\sigma _{k}e^{n-k+1,1}\end{pmatrix}}.$

Nun ist klar, wie man eine Matrix $A\in \mathbb {R} ^{n\times k}$ in die Form $QS$ mit einer orthogonalen Matrix $Q\in R^{n\times n}$ und eine Matrix $S\in \mathbb {R} ^{n\times k}$ der Gestalt (3.39) zerlegen kann. Ausgehend von $A^{(1)}:=A$ werden sukzessive Matrizen der Form

(3.52) $A^{(j)}:={\begin{pmatrix}a_{11}^{(j)}&...&...&...&...&a_{1k}^{(j)}\\&\ddots &&&&\vdots \\&&a_{j-1,j-1}^{(j)}&...&...&a_{j-1,k}^{(j)}\\&&&a_{jj}^{(j)}&...&a_{jk}^{(j)}\\&&&\vdots &\ddots &\vdots \\&&&a_{nj}^{(j)}&...&a_{nk}^{(j)}\end{pmatrix}}\in \mathbb {R} ^{n\times k},\quad j=2,...,k+1$

bestimmt, so dass dann schließlich $S:=A^{(k+1)}$ die gewünschte Gestalt hat. Die Matrizen in (3.52) werden dabei sukzessive durch Transformationen der Form

(3.53) $A^{(j+1)}:={\mathcal {P}}^{(j)}A^{(j)},\quad {\mathcal {P}}^{(j)}:={\begin{pmatrix}I_{j-1}&0\\0&{\mathcal {H}}^{(j)}\end{pmatrix}}$

mit Householder-Matrizen

${\mathcal {H}}^{(j)}:=I_{n-(j-1)}-2\omega ^{(j)}(\omega ^{(j)})^{T}$

erzeugt, wobei $\omega ^{(j)}\in \mathbb {R} ^{n-(j-1)}$ mit

${\hat {a}}^{(j)}:=(a_{jj}^{(j)},...,a_{nj}^{(j)})^{T}$

so zu wählen ist, dass mit einem $\sigma _{j}\in \mathbb {R}$

${\mathcal {H}}^{(j)}{\hat {a}}(j)=-\sigma _{j}e^{n-(j-1),1}$

gilt. Im Fall ${\hat {a}}^{(j)}\notin \operatorname {span} \{e^{n-(j-1),1}\}$ erreicht man dies gemäß Lemma 3.28 mit

$\omega ^{(j)}:={\frac {{\hat {a}}^{(j)}+\sigma _{j}e^{n-(j-1),1}}{\|{\hat {a}}^{(j)}+\sigma _{j}e^{n-(j-1),1}\|_{2}}},\quad \sigma _{j}:=\operatorname {sgn}({\hat {a}}_{1}^{(j)})\left\|{\hat {a}}^{(j)}\right\|_{2}.$

Im selteneren Fall ${\hat {a}}^{(j)}\in \operatorname {span} \{e^{n-(j-1),1}\}$ kann man offenbar ${\mathcal {H}}^{(j)}:=I_{n-(j-1)}$ bzw. ${\mathcal {P}}^{(j)}:=I$ in (3.53) wählen. Nach Lemma 3.27 sind dann die Matrizen ${\mathcal {H}}^{(1)},...,{\mathcal {H}}^{(k)}$ und damit die Matrizen ${\mathcal {P}}^{(1)},...,{\mathcal {P}}^{(k)}$ orthogonal und symmetrisch, so dass man wegen

$A^{(1)}=A,$

$A^{(2)}={\mathcal {P}}^{(1)}A^{(1)}={\mathcal {P}}^{(1)}A,$

$A^{(3)}={\mathcal {P}}^{(2)}A^{(2)}={\mathcal {P}}^{(2)}{\mathcal {P}}^{(1)}A,$

$\vdots$

$A^{(k+1)}={\mathcal {P}}^{(k)}A^{(k)}={\mathcal {P}}^{(k)}{\mathcal {P}}^{(k-1)}\cdots {\mathcal {P}}^{(1)}A$

mit

$S:={\mathcal {P}}^{(k)}{\mathcal {P}}^{(k-1)}\cdots {\mathcal {P}}^{(1)}A,\quad Q:={\mathcal {P}}^{(1)}{\mathcal {P}}^{(2)}\cdots {\mathcal {P}}^{(k)}$

wegen (3.47) die gewünschte Zerlegung $A=QS$ erhält. Gemäß Lemma 3.26 ist $Q$ tatsächlich eine orthogonale Matrix.

Man beachte, dass man $Q$ für die Lösung des Gleichungssystems $Ax=b$ mit $n=k$ (sowie für die des Ausgleichsproblems in Abschnitt 4.3) gar nicht benötigt, sondern nur den Vektor $Q^{T}b$. Wegen $\left({\mathcal {P}}^{(j)}\right)^{T}={\mathcal {P}}^{(j)}$ ist dabei

$Q^{T}=\left({\mathcal {P}}^{(k)}\right)^{T}\left({\mathcal {P}}^{(k-1)}\right)^{T}\cdots \left({\mathcal {P}}^{(1)}\right)^{T}={\mathcal {P}}^{(k)}{\mathcal {P}}^{(k-1)}\cdots {\mathcal {P}}^{(1)},$

so dass man mit $b$ wie mit $A$ verfahren kann:

$b^{(1)}:=b,$

$b^{(2)}:={\mathcal {P}}^{(1)}b^{(1)}={\mathcal {P}}^{(1)}b,$

$b^{(3)}:={\mathcal {P}}^{(2)}b^{(2)}={\mathcal {P}}^{(2)}{\mathcal {P}}^{(1)}b,$

$\vdots$

$b^{(k+1)}:={\mathcal {P}}^{(k)}b^{(k)}{\mathcal {P}}^{(k)}{\mathcal {P}}^{(k-1)}\cdots {\mathcal {P}}^{(1)}b=Q^{T}b.$

Man beachte, dass mit

$B^{(j)}:={\begin{pmatrix}a_{jj}^{(j)}&...&a_{jk}^{(j)}\\\vdots &\ddots &\vdots \\a_{nj}^{(j)}&...&a_{nk}^{(j)}\end{pmatrix}}$

gemäß (3.52) und (3.53)

$A^{(j+1)}={\begin{pmatrix}I_{j-1}&0\\0&{\mathcal {H}}^{(j)}\end{pmatrix}}{\begin{pmatrix}\#&*\\0&B^{(j)}\end{pmatrix}}={\begin{pmatrix}\#&*\\0&{\mathcal {H}}^{(j)}B^{(j)}\end{pmatrix}}$

gilt, also wie beim Gauß-Algorithmus nur die verbleibende Restmatrix $B^{(j)}$ in $A^{(j)}$ und analog der verbleibende Anteil der rechten Seite zu transformieren ist. Weiter sei gesagt, dass Householder-Transformationsmatrizen niemals explizit durch Matrixmultiplikation gebildet werden sollten. Denn eine Matrixmultiplikation ${\mathcal {H}}B$ mit ${\mathcal {H}}\in \mathbb {R} ^{m\times m}$ und $B\in \mathbb {R} ^{m\times \ell }$ kostet ${\mathcal {O}}(m^{2}\ell )$ Multiplikationen. Die benötigten Matrixmultiplikationen für ${\mathcal {H}}:=I-2\omega \omega ^{T}$ berechne man daher wie folgt:

(3.54) ${\mathcal {H}}B=(I-2\omega \omega ^{T})B=B-\omega v^{T},\quad v^{T}:=2\omega ^{T}B.$

Zunächst ermittele man also den Vektor $v$ und dann "aufdatiere" man die Matrix $B$. Für ein solches Vorgehen benötigt man nur $2m\ell$ Multiplikationen.

Bemerkung 3.31

Will man die Vektoren $\omega ^{(j)}$ aufbewahren, weil man z. B. ein Gleichungssystem $Ax=b$ mit derselben Matrix und verschiedenen rechten Seiten zu lösen hat, so kann man für jedes $j\in \{1,...,k\}$ das Diagonalelement $a_{jj}^{(j+1)}$ zunächst gesondert speichern und anschließend den Vektor $\omega ^{(j)}$ in der $j$-ten Spalte der Matrix $A^{(j+1)}$ eintragen (s. Beispiel 3.32). Auf diese Weise spart man Speicherplatz, was bei sehr großen Matrizen relevant ist.

Beispiel 3.32

Man berechne eine QS-Zerlegung von $A$ und $Q^{T}b$ für

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

Wir setzen $A^{(1)}:=A$ und $b^{(1)}:=b$ und bekommen zunächst

${\hat {a}}^{(1)}:=(1,1,1,1)^{T},\quad \left\|{\hat {a}}^{(1)}\right\|_{2}=2,\quad \sigma _{1}=\operatorname {sgn}(1)\cdot 2=2.$

Damit folgt

${\tilde {\omega }}^{(1)}:={\hat {a}}^{(1)}+\sigma _{1}e^{4,1}=(1,1,1,1)^{T}+2(1,0,0,0)^{T}=(3,1,1,1)^{T},$

$\left\|{\tilde {\omega }}^{(1)}\right\|_{2}=2{\sqrt {3}},$

$\omega ^{(1)}:={\frac {{\tilde {\omega }}^{(1)}}{\left\|{\tilde {\omega }}^{(1)}\right\|_{2}}}={\frac {1}{2{\sqrt {3}}}}(3,1,1,1)^{T}.$

Für ${\mathcal {H}}^{(1)}:=I-2\omega ^{(1)}(\omega ^{(1)})^{T}$ benötigt man nun die Produkte ${\mathcal {H}}^{(1)}A^{(1)}$ und ${\mathcal {H}}^{(1)}b^{(1)}$, die man analog zu (3.54) berechnet. Es ist

$\left(v^{(1)}\right)^{T}:=2\left(\omega ^{(1)}\right)^{T}A^{(1)}={\frac {1}{\sqrt {3}}}{\begin{pmatrix}3&1&1&1\end{pmatrix}}{\begin{pmatrix}1&0\\1&3\\1&4\\1&7\end{pmatrix}}={\frac {1}{\sqrt {3}}}{\begin{pmatrix}6&14\end{pmatrix}},$

$\left(u^{(1)}\right)^{T}:=2\left(\omega ^{(1)}\right)^{T}b^{(1)}={\frac {1}{\sqrt {3}}}{\begin{pmatrix}3&1&1&1\end{pmatrix}}{\begin{pmatrix}1\\2\\6\\4\end{pmatrix}}={\frac {15}{\sqrt {3}}}$

und demnach

$\omega ^{(1)}\left(v^{(1)}\right)^{T}={\frac {1}{6}}{\begin{pmatrix}3\\1\\1\\1\end{pmatrix}}{\begin{pmatrix}6&14\end{pmatrix}}={\begin{pmatrix}3&7\\1&7/3\\1&7/3\\1&7/3\end{pmatrix}},$

$\omega ^{(1)}\left(u^{(1)}\right)^{T}={\frac {15}{6}}{\begin{pmatrix}3\\1\\1\\1\end{pmatrix}}={\begin{pmatrix}15/2\\5/2\\5/2\\5/2\end{pmatrix}}.$

Demzufolge ergibt sich

$A^{(2)}:={\mathcal {H}}^{(1)}A^{(1)}=A^{(1)}-\omega ^{(1)}\left(v^{(1)}\right)^{T}={\begin{pmatrix}-2&-7\\0&2/3\\0&5/3\\0&14/3\end{pmatrix}},$

$b^{(2)}:={\mathcal {H}}^{(1)}b^{(1)}=b^{(1)}-\omega ^{(1)}\left(u^{(1)}\right)^{T}={\begin{pmatrix}-13/2\\-1/2\\7/2\\3/2\end{pmatrix}}.$

Wollte man den Vektor $\omega ^{(1)}$ aufbewahren, um damit zu einem späteren Zeitpunkt die Zerlegung von $A$ wieder erzeugen zu können, so legt man einen Arbeitsvektor $d\in \mathbb {R} ^{2}$ an, setzt man $d_{1}:=a_{11}^{(2)}=-2$ und trägt man $\omega ^{(1)}$ in die erste Spalte von $A^{(2)}$ ein, so dass sich ergibt:

${\hat {A}}^{(2)}:={\begin{pmatrix}{\frac {1}{2}}{\sqrt {3}}&-7\\{\frac {1}{6}}{\sqrt {3}}&2/3\\{\frac {1}{6}}{\sqrt {3}}&5/3\\{\frac {1}{6}}{\sqrt {3}}&14/3\end{pmatrix}}$

In der Praxis kann man $A^{(2)}$ und auch ${\hat {A}}^{(2)}$ wieder $A$ nennen, da man $A$ selbst nicht mehr benötigt.

Nun verfährt man analog mit der Restmatrix bzw. dem Restvektor

${\tilde {A}}^{(2)}:={\begin{pmatrix}2/3\\5/3\\14/3\end{pmatrix}},\quad {\tilde {b}}^{(2)}:={\begin{pmatrix}-1/2\\7/2\\3/2\end{pmatrix}}.$

Man bekommt

${\hat {a}}^{(2)}:=\left({\frac {2}{3}},{\frac {5}{3}},{\frac {14}{3}}\right)^{T},\quad \left\|{\hat {a}}^{(2)}\right\|_{2}=5,\quad \sigma _{2}=\operatorname {sgn} \left({\frac {2}{3}}\right)\cdot 5=5,$

${\tilde {\omega }}^{(2)}:={\hat {a}}^{(2)}+\sigma _{2}e^{3,1}=\left({\frac {2}{3}},{\frac {5}{3}},{\frac {14}{3}}\right)^{T}+5(1,0,0)^{T}=\left({\frac {17}{3}},{\frac {5}{3}},{\frac {14}{3}}\right)^{T},$

$\left\|{\tilde {\omega }}^{(2)}\right\|_{2}={\frac {\sqrt {510}}{3}},$

$\omega ^{(2)}:={\frac {{\tilde {\omega }}^{(2)}}{\left\|{\tilde {\omega }}^{(2)}\right\|_{2}}}={\frac {1}{\sqrt {510}}}(17,5,14)^{T},$

$\left(v^{(2)}\right)^{T}:=2\left(\omega ^{(2)}\right)^{T}{\tilde {A}}^{(2)}={\frac {2}{\sqrt {510}}}{\begin{pmatrix}17&5&14\end{pmatrix}}{\begin{pmatrix}2/3\\5/3\\14/3\end{pmatrix}}={\frac {510}{3{\sqrt {510}}}},$

$\left(u^{(2)}\right)^{T}:=2\left(\omega ^{(2)}\right)^{T}{\tilde {b}}^{(1)}={\frac {1}{\sqrt {510}}}{\begin{pmatrix}17&5&14\end{pmatrix}}{\begin{pmatrix}-1\\7\\3\end{pmatrix}}={\frac {60}{\sqrt {510}}},$

${\mathcal {H}}^{(2)}{\tilde {A}}^{(2)}={\tilde {A}}^{(2)}-\omega ^{(2)}\left(v^{(2)}\right)^{T}={\begin{pmatrix}2/3\\5/3\\14/3\end{pmatrix}}-{\frac {1}{3}}{\begin{pmatrix}17\\5\\14\end{pmatrix}}={\begin{pmatrix}-5\\0\\0\end{pmatrix}},$

${\mathcal {H}}^{(2)}{\tilde {b}}^{(2)}={\tilde {b}}^{(2)}-\omega ^{(2)}\left(u^{(2)}\right)^{T}={\begin{pmatrix}-1/2\\7/2\\3/2\end{pmatrix}}-{\frac {60}{510}}{\begin{pmatrix}17\\5\\14\end{pmatrix}}={\begin{pmatrix}-5/2\\99/34\\-5/34\end{pmatrix}}.$

Damit ergibt sich nun

$S:=A^{(3)}={\begin{pmatrix}-2&-7\\0&-5\\0&0\\0&0\end{pmatrix}},\quad R:={\begin{pmatrix}-2&-7\\0&-5\end{pmatrix}},\quad Q^{T}b=b^{(3)}:={\begin{pmatrix}-13/2\\-5/2\\99/34\\-5/34\end{pmatrix}}.$

Will man die $\omega ^{(j)}$ wiederverwenden, so bilde man mit ${\hat {A}}^{(2)}$ und dem oben definierten $d_{1}$

${\hat {A}}^{(3)}={\begin{pmatrix}{\frac {1}{2}}{\sqrt {3}}&-7\\{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{\sqrt {510}}}17\\{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{\sqrt {510}}}5\\{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{\sqrt {510}}}14\end{pmatrix}},\quad d:={\begin{pmatrix}-2\\-5\end{pmatrix}}.$

Aus ${\hat {A}}^{(3)}$ und $d$ könnte man $R$ in offensichtlicher Weise gewinnen.

---

Man kann sich überlegen, dass eine QS-Zerlegung von $A\in \mathbb {R} ^{n\times k}$ mittels Householder-Transformationen etwa

$2\left[{\frac {k^{3}}{3}}+(n-k){\frac {k^{2}}{2}}\right]$

Multiplikationen benötigt, also

(3.55) $\sim nk^{2}$ für $n\gg k$ und $\sim {\frac {2}{3}}n^{3}$ für $n\approx k$

und damit bei quadratischen Matrizen etwa doppelt so viele wesentliche Rechenoperationen wie der Gauß-Algorithmus. Neben dem Einsatz für die bereits besprochene stabile Lösung von linearen Gleichungssystemen werden wir im nächsten Abschnitt eine weitere wichtige Anwendung einer solchen QS-Zerlegung geben.

Siehe auch

• normalisierte Gleitkomma-Darstellung
• Lösung linearer Gleichungssysteme - (Foliensatz)
• Zerlegung PA = LR - (Foliensatz)

Seiteninformation

Diese Lernresource können Sie als Wiki2Reveal-Foliensatz darstellen.

Wiki2Reveal

Dieser Wiki2Reveal Foliensatz wurde für den Lerneinheit Kurs:Numerik I' erstellt der Link für die Wiki2Reveal-Folien wurde mit dem Wiki2Reveal-Linkgenerator erstellt.

• Die Seite wurde als Dokumententyp PanDocElectron-SLIDE erstellt.
• Link zur Quelle in Wikiversity: https://de.wikiversity.org/wiki/Kurs:Numerik%20I/Orthonormalisierungsverfahren
• siehe auch weitere Informationen zu Wiki2Reveal und unter Wiki2Reveal-Linkgenerator.