4.1 Problemstellung

In der Praxis hat man häufig auch ein überbestimmtes lineares Gleichungssystem $Ax=b$ für eine Matrix $A\in \mathbb {R} ^{n\times k}$ mit $n>k$ und eine rechte Seite $b\in \mathbb {R} ^{n}$ zu „lösen“. Da ein solches System mehr Gleichungen als Unbekannte hat, ist im Allgemeinen $Ax-b\neq 0$ für alle $x\in \mathbb {R} ^{k}$ und besitzt es somit keine exakte Lösung. Es macht also Sinn, ein $x^{*}$ als „Lösung“ zu akzeptieren, für welches der Defekt $Ax-b$ hinsichtlich einer gewählten $l_{p}$ -Norm $\|\cdot \|_{p}$ auf dem $\mathbb {R} ^{k}$ minimal wird, also eine Lösung $x^{*}$ des Problems

\inf _{x\in \mathbb {R} ^{k}}\|b-Ax\|_{p}.

Im Fall der Verwendung der $l_{2}$ - bzw. Euklidischen Norm lautet dieses Problem

(4.1)

\inf _{x\in \mathbb {R} ^{k}}\|b-Ax\|_{2},

welches im Hinblick auf eine Lösung mit dem Problem

(4.2)

\inf _{x\in \mathbb {R} ^{k}}\|b-Ax\|_{2}^{2}

äquivalent ist, wobei die zu minimierende Funktion in (4.2) für alle $x\in \mathbb {R} ^{k}$ differenzierbar ist. (Die Funktionen in (4.1) und (4.2) haben offenbar dieselben Minimalpunkte $x$ , sofern solche existieren.) Bei Wahl der $l_{2}$ -Norm minimiert man also die Summe der Fehlerquadrate, und man spricht daher auch von Fehlerquadratmethode oder von diskreter $l_{2}$ -Approximation. Die Lösung des entsprechenden Problems für die $l_{1}$ - und die $l_{\infty }$ - bzw. Tschebyscheff-Norm ist ebenfalls von großem praktischen Interesse, führt aber auf nichtdifferenzierbare Funktionen, so dass diese Probleme schwieriger zu lösen sind. (Man kann letztere Probleme als lineare Optimierungsprobleme formulieren und beispielsweise mit dem sog. Simplexalgorithmus lösen.) Bevor wir nun das Problem in (4.2) weiter untersuchen, wollen wir zunächst zwei Aufgabenstellungen beschreiben und analysieren, die auf ein derartiges über-bestimmtes Gleichungssystem führen.

In Naturwissenschaft und Technik hat man oft das Problem, mit einer großen Zahl von Messwerten umgehen zu müssen. Ein anderes, sich häufig stellendes Problem ist es, eine in endlich vielen Punkten gegebene Funktion, welche durch eine rechenaufwändige Vorschrift bestimmt ist, durch eine einfacher zu berechnende zu ersetzen. Beide Probleme kann man gemeinsam angehen und als $l_{2}$ -Approximationsprobleme beschreiben.

Wir gehen dazu von $n$ Zahlenpaaren

(4.3)

(t_{j},y_{j}),\quad j=1,\ldots ,n

mit $t_{r}\neq t_{s}$ für $r\neq s$ aus, wobei üblicherweise $n$ groß ist. Beispielsweise können die $y_{j}\in \mathbb {R}$ irgendwelche zu unterschiedlichen Zeitpunkten $t_{j}\in \mathbb {R}$ gemessene Werte oder, im Hinblick auf die Approximation einer gegebenen Funktion $f$ , die Funktionswerte $y_{j}:=f(t_{j})$ zu gewissen Punkten $t_{j}$ des Definitionsbereichs von $f$ sein. Das Ziel der Ausgleichsrechnung ist es nun, durch geeignete Wahl eines Parametervektors $x:=(x_{1},\ldots ,x_{k})^{T}$ eine Funktion der Gestalt

(4.4)

z(x,t):=x_{1}v_{1}(t)+x_{2}v_{2}(t)+\ldots +x_{k}v_{k}(t),\quad t\in \mathbb {R}

mit $k$ gegebenen stetigen Ansatzfunktionen $v_{i}$ zu finden, so dass die Fehler

(4.5)

y_{j}-z(x,t_{j}),\quad j=1,\ldots ,n

möglichst klein ausfallen. Dabei sollte sinnvollerweise $k\leq n$ sein und ist zumeist $k\ll n$ . Hat man einen solchen Parametervektor $x$ bzw. eine solche Funktion $z(x,\cdot )$ gefunden und sind die Approximationsfehler in (4.5) nicht zu groß, so kann man statt mit den Daten (4.3) nur mit dieser Funktion arbeiten, für die im Fall $k\ll n$ erheblich weniger Information, nämlich nur der Vektor $x$ , abgespeichert werden muss. Weil $z(x,\cdot )$ eine stetige Funktion ist, erlaubt ein solches Vorgehen außerdem, Werte $y$ zu Werten bzw. „Zeiten“ $t$ zu berechnen, für die keine Messung vorliegt.

Die Ansatzfunktionen hat man so zu wählen, dass sie den gemessenen Prozess möglichst gut beschreiben. Häufig, insbesondere dann, wenn man wenig über den gegebenen Prozess weiß, verwendet man die Monome

(4.6)

v_{i}(t):=t^{i-1}\quad (i=1,\ldots ,k)

so dass

z(x,t):=x_{1}+x_{2}t+\ldots +x_{k}t^{k-1},\quad t\in \mathbb {R}

ein Polynom vom Grad $\leq k-1$ ist (Polynomapproximation). Wenn es sich um einen periodischen Prozess handelt, ist es aber beispielsweise günstiger, die $k:=2p+1$ Funktionen

(4.7)

1,\quad \sin(t),\quad \cos(t),\quad \sin(2t),\quad \cos(2t),\quad \ldots ,\quad \sin(pt),\quad cos(pt)

als Ansatzfunktionen zu wählen (trigonometrische Approximation), weil man dann im Allgemeinen bei gleicher Anzahl von Funktionen kleinere Fehler erhält. Andere Systeme von Ansatzfunktionen können ebenfalls vernünftig sein.

Nun ist es nicht sinnvoll, $x$ so zu wählen, dass die Summe aller Fehler

y_{j}-z(x,t_{j}),\quad j=1,\ldots ,n

möglichst klein wird, da diese Summe auch bei großen Einzelfehlern sehr klein werden kann, nämlich dann, wenn sich die positiven und negativen Fehler (nahezu) aufheben. Die Größe des Fehlervektors $(y_{j}-z(x,t_{j}))_{j=1,\ldots ,n}$ misst man daher mit einer $l_{p}$ -Norm im $\mathbb {R} ^{n}$ . Insbesondere führt dann die Verwendung der quadrierten $l_{2}$ - bzw. Euklidischen Norm (siehe den Kommentar zu (4.1) und (4.2)) auf die für alle $x\in \mathbb {R} ^{k}$ differenzierbare Funktion

{\tilde {F}}(x):=\sum _{j=1}^{n}(y_{j}-z(x,t_{j}))^{2}.

Mit den Setzungen

(4.8)

b:={\begin{pmatrix}y_{1}\\\vdots \\y_{j}\\\vdots \\y_{n}\end{pmatrix}},\quad A:={\begin{pmatrix}v_{1}(t_{1})&\ldots &v_{k}(t_{1})\\\vdots &&\vdots \\v_{1}(t_{j})&\ldots &v_{k}(t_{j})\\\vdots &&\vdots \\v_{1}(t_{n})&\ldots &v_{k}(t_{n})\end{pmatrix}}

kann diese geschrieben werden als

{\tilde {F}}(x):=\sum _{j=1}^{n}(y_{j}-z(x,t_{j}))^{2}=\sum _{j=1}^{n}(b_{j}-(Ax)_{j})^{2}=\|b-Ax\|_{2}^{2}.

Das beschriebene Problem der Ausgleichsrechnung ist also von der Form (4.2), wobei $A$ und $b$ durch (4.8) gegeben sind.

Wir betrachten nun allgemein das Problem (4.2). Die zu minimierende Funktion

F(x):=\|b-Ax\|_{2}^{2}

können wir wegen

F(x)=(b-Ax)^{T}(b-Ax)=b^{T}b-(Ax)^{T}b-b^{T}Ax+(Ax)^{T}Ax

als quadratische Funktion in $k$ Veränderlichen

(4.9)

F(x)={\frac {1}{2}}x^{T}(2A^{T}A)x-(2A^{T}b)^{T}x+b^{T}b

schreiben. Für die darin vorkommende Matrix $A^{T}A$ kann man aussagen:

Lemma 4.1

Sei $A\in \mathbb {R} ^{n\times k}$ mit $n\geq k$ und $\operatorname {Rang} (A)=k$ . Dann ist die Matrix $A^{T}A\in \mathbb {R} ^{k\times k}$ positiv definit.

Beweis.

Die Matrix $A^{T}A$ ist wegen $(A^{T}A)^{T}=A^{T}A$ symmetrisch. Weiter ist sie positiv semidefinit, d. h. es gilt

h^{T}(A^{T}A)h=(Ah)^{T}(Ah)=\|Ah\|_{2}^{2}\geq 0,\quad h\in \mathbb {R} ^{k}.

Wegen $\operatorname {Rang} (A)=k$ sind die $k$ Spalten $a^{1},\ldots ,a^{k}$ von $A$ linear unabhängig (Zeilenrang = Spaltenrang) und daher hat man

h^{T}(A^{T}A)h=0\Leftrightarrow \|Ah\|_{2}=0\Leftrightarrow Ah=0\Leftrightarrow h=0.

q.e.d.

Im Fall der Ausgleichsrechnung mit den Ansatzfunktionen (4.6) und (4.7) ist die Rangbedingung in Lemma 4.1 unter unserer Voraussetzung $n\geq k$ immer erfüllt, wie das folgende Lemma besagt.

Lemma 4.2

Für die mit den Ansatzfunktionen (4.6) oder (4.7) gebildete Matrix $A\in \mathbb {R} ^{n\times k}$ in (4.8) mit $n\geq k$ gilt $\operatorname {Rang} (A)=k$ .

Beweis.

Für $A$ in (4.8) und $h:=(h_{1},\ldots ,h_{k})^{T}\in \mathbb {R} ^{k}$ gilt:

Ah=0\Leftrightarrow \sum _{i=1}^{k}h_{i}v_{i}(t_{j})=0\quad (j=1,\ldots ,n).

Wird $A$ insbesondere durch (4.6) spezifiziert, so implizieren wegen $n\geq k$ letztere Gleichungen die Identität $h=0$ , da ein von Null verschiedenes Polynom vom Grad $\leq k-1$ höchstens $k-1$ Nullstellen besitzt. Für (4.7) schließt man analog mittels komplexer Darstellungen des Sinus und Kosinus (siehe z. B. Collatz/Krabs: Approximationstheorie, Teubner, Stuttgart, 1973).

q.e.d.

4.2 Lösung der Normalgleichungen

Wir betrachten nun die quadratische Funktion $F:\mathbb {R} ^{k}\to \mathbb {R}$ in (4.9), wobei wir $\operatorname {Rang} (A)=k$ voraussetzen. Sie hat bei $x$ den Gradienten und die Hessematrix

\nabla F(x)=2A^{T}Ax-2A^{T}b,\quad \nabla ^{2}F(x)=2A^{T}A.

Notwendige Bedingung dafür, dass $x^{*}$ Minimalpunkt von $F$ ist, ist die Bedingung $\nabla F(x^{*})=0$ bzw. äquivalent dazu, dass $x^{*}$ die sog. Normalgleichungen

A^{T}Ax=A^{T}b

erfüllt. Nach Lemma 4.1 ist dabei die (von $x$ unabhängige) Matrix $\nabla ^{2}F(x)$ positiv definit, so dass die eindeutige Lösung $x^{*}$ der Normalgleichungen auch der einzige (globale) Minimalpunkt von $F$ ist. Wir können somit festhalten:

Satz 4.3

Sei $A\in \mathbb {R} ^{n\times k}$ mit $n\geq k$ und $\operatorname {Rang} (A)=k$ . Dann besitzt das lineare Ausgleichsproblem

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

eine eindeutige Lösung $x^{*}\in \mathbb {R} ^{k}$ und diese ist eindeutige Lösung des linearen Gleichungssystems

(4.10)

A^{T}Ax=A^{T}b.

Wir betrachten dazu ein Beispiel der Ausgleichsrechnung.

Beispiel 4.4

Wir betrachten den Fall der sog. Ausgleichsgeraden. Wenn die $y_{j}$ $(j=1,\ldots ,n)$ mit $n\geq 2$ in (4.3) ungefähr auf einer Geraden liegen, macht es Sinn,

v_{1}(t):=1,\quad v_{2}(t):=t

d. h. die Ansatzfunktionen (4.6) mit $k=2$ zu wählen und somit

z(x,t):=x_{1}+x_{2}t,\quad t\in \mathbb {R}

zu setzen. Mit

e:=(1,1,\ldots ,1)^{T}\in \mathbb {R} ^{n},\quad d:=(t_{1},t_{2},\ldots ,t_{n})^{T}\in \mathbb {R} ^{n}

hat $A\in \mathbb {R} ^{n\times 2}$ in diesem Fall die Gestalt $A={\begin{pmatrix}e&d\end{pmatrix}}$ . Weiter ist dann

A^{T}A={\begin{pmatrix}e^{T}\\d^{T}\end{pmatrix}}{\begin{pmatrix}e&d\end{pmatrix}}={\begin{pmatrix}e^{T}e&e^{T}d\\e^{T}d&d^{T}d\end{pmatrix}},\quad A^{T}b={\begin{pmatrix}e^{T}b\\d^{T}b\end{pmatrix}}.

Für eine symmetrische invertierbare Matrix $B\in \mathbb {R} ^{2\times 2}$ kann man die Inverse explizit angeben:

B^{-1}={\begin{pmatrix}b_{11}&b_{12}\\b_{12}&b_{22}\end{pmatrix}}^{-1}={\frac {1}{b_{11}b_{22}-b_{12}^{2}}}{\begin{pmatrix}b_{22}&-b_{12}\\-b_{12}&b_{11}\end{pmatrix}}.

Somit lautet die Lösung der Normalgleichungen (4.10) in diesem Fall

x^{*}:=\left(A^{T}A\right)^{-1}A^{T}b={\frac {1}{(e^{T}e)(d^{T}d)-(e^{T}d)^{2}}}{\begin{pmatrix}d^{T}d&-e^{T}d\\-e^{T}d&e^{T}e\end{pmatrix}}{\begin{pmatrix}e^{T}b\\d^{T}b\end{pmatrix}}

und demzufolge

(4.11)

x^{*}={\frac {1}{(e^{T}e)(d^{T}d)-(e^{T}d)^{2}}}{\begin{pmatrix}\left(d^{T}d\right)\left(e^{T}b\right)-\left(d^{T}b\right)\left(e^{T}d\right)\\(e^{T}e)(d^{T}b)-(e^{T}d)(e^{T}b)\end{pmatrix}}.

Dabei hat man

(4.12)

e^{T}e=n,\quad e^{T}d=\sum _{j=1}^{n}t_{j},\quad e^{T}b=\sum _{j=1}^{n}y_{j},\quad d^{T}d=\sum _{j=1}^{n}t_{j}^{2},\quad d^{T}b=\sum _{j=1}^{n}t_{j}y_{j}.

Beispielsweise für die $n=8$ Wertepaare

{\begin{array}{|l||c|c|c|c|c|c|c|c|}\hline t_{j}&1&2&3&4&5&6&7&8\\\hline y_{j}&1.75&2.18&2.63&3.24&3.69&4.16&4.55&5.29\\\hline \end{array}}

errechnet man

\sum _{j=1}^{8}t_{j}=36,\quad \sum _{j=1}^{8}y_{j}=27.49,\quad \sum _{j=1}^{8}t_{j}^{2}=204,\quad \sum _{j=1}^{8}t_{j}y_{j}=144.54.

Mit (4.11) und (4.12) ergibt sich somit

x^{*}={\frac {1}{8\cdot 204-36^{2}}}{\begin{pmatrix}204\cdot 27.49-36\cdot 144.54\\8\cdot 144.54-36\cdot 27.49\end{pmatrix}}={\begin{pmatrix}1.203\ 929\\0.496\ 071\end{pmatrix}}.

Die Ausgleichsgerade zu den gegebenen Daten lautet folglich

z(x^{*},t):=1.203\ 929+0.496\ 071t,\quad t\in \mathbb {R} .

Der maximale relative Fehler der $z(x^{*},t_{j})$ bezüglich der $y_{j}$ beträgt

\max _{1\leq j\leq 8}{\frac {|y_{j}-z(x^{*},t_{j})|}{|z(x^{*},t_{j})|}}=0.016\ 243

bzw. ungefähr 1.6%.

Für $k>2$ könnte man die Normalgleichungen (4.10) mittels einer Cholesky-Zerlegung lösen. Diese selbst ist, wie man zeigen kann, numerisch stabil. Leider ist das Ausgleichproblem selbst aber häufig schlecht konditioniert. Man betrachte z. B. die Matrix $A$ , die sog. Vandermonde-Matrix, die man für $n=k$ im Fall der Wahl der Monome (4.6) als Ansatzfunktionen erhält:

A:={\begin{pmatrix}1&t_{1}&\ldots &t_{1}^{k-1}\\1&t_{2}&\ldots &t_{2}^{k-1}\\\vdots &\vdots &\ddots &\vdots \\1&t_{k}&\ldots &t_{k}^{k-1}\end{pmatrix}}.

Für $t\in [0,1]$ unterscheiden sich die Funktionen $t^{r-1}$ und $t^{r}$ bereits für nicht allzu großes $r$ kaum, so dass die $r$ -te und $(r+1)$ -te Spalten in $A$ für solche $r$ nahezu linear abhängig sind. Die oft große Kondition von $A$ geht außerdem noch im Fall $n=k$ bei der Lösung der Normalgleichungen quadratisch ein, denn es gilt:

Lemma 4.5

Für eine reguläre Matrix $A\in \mathbb {R} ^{k\times k}$ gilt

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

Beweis.

Nach Satz 3.19 hat eine positiv definite Matrix $B\in \mathbb {R} ^{k\times k}$ Eigenwerte $\lambda _{i}:=\lambda _{i}(B)>0$ $(i=1,\ldots ,k)$ . Weiter hat wegen

Bx^{i}=\lambda _{i}x^{i}\Leftrightarrow B^{-1}x^{i}=[\lambda _{i}]^{-1}x^{i}

für Eigenvektoren $x^{i}$ zu $\lambda _{i}$ die Inverse $B^{-1}$ , die somit auch positiv definit ist, die Eigenwerte $[\lambda _{i}]^{-1}$ $(i=1,\ldots ,k)$ . Es gilt folglich nach Satz 2.15

\operatorname {cond} _{2}(B)=\|B\|_{2}\left\|B^{-1}\right\|_{2}=\lambda _{\max }/\lambda _{\min },

wobei $\lambda _{\max }:=\lambda _{\max(}B)$ und $\lambda _{\min }:=\lambda _{\min(}B)$ einen größten und kleinsten Eigenwert von $B$ bezeichnen. Indem man $x$ mittels einer Orthonormalbasis von Eigenvektoren darstellt, kann man ferner die Abschätzungen

\lambda _{\min }\|x\|_{2}^{2}\leq x^{T}Bx\leq \lambda _{\max }\|x\|_{2}^{2},\quad x\in \mathbb {R} ^{k}

beweisen, wobei offenbar Gleichheit in der ersten bzw. zweiten Ungleichung für einen zu $\lambda _{\min }$ bzw. $\lambda _{\max }$ gehörenden Eigenvektor angenommen wird. Folglich schließt man

\lambda _{\min }=\min _{\|x\|_{2}=1}x^{T}Bx,\quad \lambda _{\max }=\max _{\|x\|_{2}=1}x^{T}Bx.

Wendet man diese Ergebnisse auf die nach Lemma 4.1 positiv definite Matrix $A^{T}A\in \mathbb {R} ^{k\times k}$ an, so erhält man mit Satz 2.19

\operatorname {cond} _{2}(A^{T}A)={\frac {\lambda _{\max(}A^{T}A)}{\lambda _{\min(}A^{T}A)}}={\frac {\max \limits _{\|x\|_{2}=1}x^{T}(A^{T}A)x}{\min \limits _{\|x\|_{2}=1}x^{T}(A^{T}A)x}}={\frac {\max \limits _{\|x\|_{2}=1}\|Ax\|_{2}^{2}}{\min \limits _{\|x\|_{2}=1}\|Ax\|_{2}^{2}}}=[\operatorname {cond} _{2}(A)]^{2}.

q.e.d.

Es ist daher große Vorsicht bei Anwendung der Cholesky-Zerlegung für die Lösung der Normalgleichungen geboten. Prinzipiell ist sie nur zu empfehlen, wenn große Residuen $b_{i}-(Ax)_{i}$ $(i=1,\ldots ,n)$ in der Lösung des Ausgleichsproblems zu erwarten sind (s. Deuflhard/Hohmann). Sicherer ist es aber, so vorzugehen, wie es im folgenden Abschnitt beschrieben ist.

4.3 Lösung mittels QS-Zerlegung

Für eine QS-Zerlegung stellen wir zunächst fest:

Lemma 4.6

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 (vgl. (3.39)):

(4.13)

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

Dann gilt:

(i) $\det(R)\neq 0$ ,

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

Beweis.

Da $Q$ eine orthogonale Matrix ist, gilt

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

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

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

also (i) richtig. Weiter hat man

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

Wendet man nun Lemma 4.5 auf die nach Lemma 4.1 positiv definite Matrix $R^{T}R$ an, so gelangt man zu (ii).

q.e.d.

Der folgende Satz gibt an, wie Fehlerquadratprobleme mittels einer QS-Zerlegung von $A$ gelöst werden können.

Satz 4.7

Es seien $A\in \mathbb {R} ^{n\times k}$ mit $n\geq k$ und $\operatorname {Rang} (A)=k$ , eine orthogonale Matrix $Q\in \mathbb {R} ^{n\times n}$ und $S\in \mathbb {R} ^{n\times k}$ gegeben wie in Lemma 4.6. 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},\quad y^{1}\in \mathbb {R} ^{k},\quad y^{2}\in \mathbb {R} ^{n-k}.

Dann ist der Vektor $x^{*}\in \mathbb {R} ^{k}$ genau dann die eindeutige Lösung des linearen Ausgleichsproblems

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

wenn $x^{*}$ eindeutige die Lösung des Systems

Rx=y^{1}

ist. Außerdem hat man

\|b-Ax^{*}\|_{2}=\left\|y^{2}\right\|_{2}.

Beweis.

Für einen beliebigen Vektor $x\in \mathbb {R} ^{k}$ erhalten wir unter Verwendung von (3.38) und Lemma 3.23

\|b-Ax^{*}\|_{2}^{2}=\left\|(QQ^{T})b-QSx\right\|_{2}^{2}=\left\|Q^{T}b-Sx\right\|_{2}^{2}=\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}.

Daraus folgt

\|b-Ax\|_{2}\geq \left\|y^{2}\right\|_{2}

sowie

\|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.

Nach Satz 4.3 und Lemma 4.6 besitzen das $l_{2}$ -Approximationsproblem sowie das System $Rx=y^{1}$ eine eindeutige Lösung. Wir geben zum letzten Satz ein Beispiel.

Beispiel 4.8

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}}.

Nach Beispiel 3.32 liefert eine QS-Zerlegung von $A$ mit gleichzeitiger Berechnung von $Q^{T}b$

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

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}-13/2\\-5/2\end{pmatrix}}

lautet

x^{*}=(x_{1}^{*},x_{2}^{*})^{T}:=(3/2,1/2)^{T}.

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

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

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 Satz 4.7 dargestellten Weg, vergleichen. 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 $k^{2}/2$ Multiplikationen und Divisionen erfordert, vernachlässigen wir hier diesen zusätzlichen Aufwand bei der 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 $b-Ax$ weitere $nk$ Multiplikationen, so dass die zuletzt genannten Teilaufgaben zusammen ungefähr

{\frac {1}{2}}nk^{2}+{\frac {1}{6}}k^{3}+nk

wesentliche Rechenoperationen erfordern, das sind

(4.14)

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

für

n\gg k

und

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

für

n\approx k

.

Für das sich aus Satz 4.7 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$ . Wie ein Vergleich von (3.55) und (4.14) liefert, müssen demnach für beide Wege im Fall $n\approx k$ etwa gleich viele wesentliche Rechenoperationen durchgeführt werden, während der Weg über die QS-Zerlegung für $n\gg k$ etwa doppelt so viele wesentliche Rechenoperationen erfordert wie der über die Normalgleichungen. 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 4.9

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}}.

Damit ergibt sich

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}}.

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 :=.5_{10}-5$ und damit $\varepsilon ^{2}=.25_{10}-10$ . Dann erhält man

gl(1+\varepsilon ^{2})=1,\quad gl(A^{T}A)=(a_{ij})

mit

a_{ij}:=1

(i,j=1,\ldots ,5)

.

Die Matrix $gl(A^{T}A)$ hat offenbar den Rang 1 und ist singulär. Man errechnet hier

\operatorname {cond} _{2}(ATA)={\frac {5+\varepsilon ^{2}}{\varepsilon ^{2}}}\approx 2_{10}+11,

\left[\operatorname {cond} _{2}(A^{T}A)\right]^{1/2}=\left[{\frac {5+\varepsilon ^{2}}{\varepsilon ^{2}}}\right]^{1/2}\approx 4.5_{10}+5.

This article is issued from Wikiversity. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.