[WS] Lineare Ausgleichprobleme

Neue Frage »

26.07.2007, 02:13

tigerbine

Auf diesen Beitrag antworten »

[WS] Lineare Ausgleichprobleme

Gliederung

Wahl der Norm
Normalengleichungen

2a. Existenz einer Lösung

2b. Eindeutigkeit für Rang(A)=n
Lösen mit Cholesky-Zerlegung
QR-Zerlegungen

4a. Eindeutigkeit der reduzierten QR-Zerlegung

4b. Lösen der Normalengleichung - red. QR

4c. Lösen des Ausgleichsproblems - volle QR
Gram-Schmidt-Orthogonalisierung

5a. Berechnung der red. QR-Zerlegung

5b. Existenz der red. QR-Zerlegung

5c. Implementierung - Modifiziertes Gram-Schmidt Verfahren
Householder Spiegelungen

6a. Eigenschaften

6b. Konstruktion von Q

6c. Implementierung

6d. Geometrische Deutung
Givens-Rotationen

7a. Eigenschaften

7b. Konstruktion von Q

7c. Implementierung

7d. Geometrische Deutung
Gauß'sches Ausgleichsverfahren (Methode der kleinsten Quadrate)

8a. Lineare Modellfunktion

8b. Allgemeiner Linearer Fall

26.07.2007, 02:13

tigerbine

Auf diesen Beitrag antworten »

1. Wahl der Norm
Im Gegensatz zur Lösung von regulären Systemen $\begin{eqnarray*} Ax=b \end{eqnarray*}$ werden hier Systeme der Art $\begin{eqnarray*} Ax=b \end{eqnarray*}$ mit $\begin{eqnarray*} A \in \mathbb R^{m \times n},~m \geq n \end{eqnarray*}$ betrachtet. Im Allgemeinen wird es nun kein x* geben mit:

$\begin{eqnarray*} Ax^*=b \end{eqnarray*}$

Man könnte nun aber nach einer approximativen Lösung $\begin{eqnarray*} Ax \approx b \end{eqnarray*}$ suchen. Natürlich nicht nach irgendeiner, sondern nach derjenigen, die den Ausdruck

$\begin{eqnarray*} ||Ax-b|| \end{eqnarray*}$

minimiert. Es werden hier häufig die folgenden Normen verwendet:

$\begin{eqnarray*} ||.||_1 \end{eqnarray*}$
$\begin{eqnarray*} ||.||_2 \end{eqnarray*}$
$\begin{eqnarray*} ||.||_\infty \end{eqnarray*}$

Dabei hängt die Lösung eines Ausgleichsproblems von der verwendeten Norm ab. Ein Beispiel (Vorgriff):

Gegeben seien die sieben Punktepaare (Messreihe):

$\begin{eqnarray*} (0,0),~(\frac{1}{6},0),~(\frac{2}{6},0),(\frac{3}{6},0),(\frac{4}{6},0),(\frac{5}{6},0),(\frac{6}{6},1) \end{eqnarray*}$

Gesucht in nun eine konstante Funktion g, welche die Paare auf dem Intervall [0,1] am besten approximiert. Mit der Bezeichnung g(t):=c, sucht man also c derart, dass der Ausdruck:

$\begin{eqnarray*} \begin{Vmatrix}\begin{pmatrix} g(t_1)-y_1 \\ \vdots \\ g(t_m)-y_m \end{pmatrix} \end{Vmatrix} = \begin{Vmatrix} \begin{pmatrix}c-y_1 \\ \vdots \\ c-y_m \end{pmatrix}\end{Vmatrix} \end{eqnarray*}$

möglichst klein wird. Das soll nun für obige 3 Normen geschehen.

$\begin{eqnarray*} ||.||_1 \end{eqnarray*}$

$\begin{eqnarray*} \begin{Vmatrix} \begin{pmatrix}c-y_1 \\ \vdots \\ c-y_m \end{pmatrix}\end{Vmatrix}_1 = \sum_{i=1}^7~|c-y_i| = \sum_{i=1}^6~|c| + |c-1| = 6|c| + |c-1| \end{eqnarray*}$

Offensichtlicht lautet die Lösung dann $\begin{eqnarray*} c=0 \end{eqnarray*}$
$\begin{eqnarray*} ||.||_2 \end{eqnarray*}$

$\begin{eqnarray*} \begin{Vmatrix} \begin{pmatrix}c-y_1 \\ \vdots \\ c-y_m \end{pmatrix}\end{Vmatrix}_2 &&= \Bigl(\sum_{i=1}^m~(c-y_i)^2 \Bigr)^{\frac{1}{2}} =\Bigl(\sum_{i=1}^{m-1}~c^2 + (c-1)^2 \Bigr)^{\frac{1}{2}} \\&& = \Bigl(6c^2+ (c-1)^2 \Bigr)^{\frac{1}{2}} \end{eqnarray*}$

Dieser Ausdruck wird minimal, wenn der Term unter der Wurzel minimal wird. Mit der notwendigen und da quadratische Funktion hinreichenden Bedingung $\begin{eqnarray*} f'(c)=0 \end{eqnarray*}$ folgt dann $\begin{eqnarray*} c = \frac{1}{7} \end{eqnarray*}$
$\begin{eqnarray*} ||.||_\infty \end{eqnarray*}$

$\begin{eqnarray*} \begin{Vmatrix} \begin{pmatrix}c-y_1 \\ \vdots \\ c-y_m \end{pmatrix}\end{Vmatrix}_\infty &&= \max_{i=1,...,m} |c-y_i| = \max\{|c|,|1-c|\} \end{eqnarray*}$

Somit lautet hier die Lösung $\begin{eqnarray*} c=\frac{1}{2} \end{eqnarray*}$

Interpretation

Der letzte Wert ist wohl ein Ausreißer der Messreihe. Bei der Ausgleichsrechnung mittels $\begin{eqnarray*} \mathit l_1 \end{eqnarray*}$ Norm spielt dieser wohl keiner Rolle, bei der $\begin{eqnarray*} \mathit l_2 \end{eqnarray*}$ Norm fließt er nur geringfügig ein. Den größten Einfluss hat er wohl bei der $\begin{eqnarray*} \mathit l_\infty \end{eqnarray*}$ Norm.

Lehrer

Im weiteren soll die $\begin{eqnarray*} \mathit l_2 \end{eqnarray*}$ Norm verwendet werden, denn in der nun äquivalenten Formulierung:

$\begin{eqnarray*} \min!~||Ax-b||_2^2 \end{eqnarray*}$

liegt ein stetig differenzierbares Optimierungsproblem vor. Die beiden anderen Normen würden nicht differenzierbare Probleme liefern.

26.07.2007, 02:14

tigerbine

Auf diesen Beitrag antworten »

2. Normalengleichungen
Hier soll ein zum Minimierungsproblem äquivalentes aufgezeigt werden.

$\begin{eqnarray*} x^* \text{ löst }~\min!~||Ax-b||_2 ~\Leftrightarrow x^* \text{ löst } A^TAx = A^Tb \end{eqnarray*}$

Beweis:

$\begin{eqnarray*} \Leftarrow \end{eqnarray*}$

Es gilt nun $\begin{eqnarray*} A^TAx^*=A^Tb \end{eqnarray*}$ . Damit folgt auch die Gültigkeit von $\begin{eqnarray*} A^T(Ax^*-b)=0~(*) \end{eqnarray*}$ . Es ist nun:

$\begin{eqnarray*} ||Ax-b||_2^2 && = ||A(x-x^*) + (Ax^*-b) ||_2^2 \\ && = [A(x-x^*) + (Ax^*-b)]^T[A(x-x^*) + (Ax^*-b)] \\&& = ||A(x-x^*)||_2^2 + 2(x-x^*)^T\underbrace{A^T(Ax^*-b)}_{(*)} + ||Ax^*-b||_2^2\\ &&=||A(x-x^*)||_2^2 + ||Ax^*-b||_2^2 \\ && \geq ||Ax^*-b||_2^2 \\&&= ||Ax^*-b||_2^2 + 2td^T(-d) +t^2\cdot ||Ad||_2^2 \end{eqnarray*}$

$\begin{eqnarray*} \Rightarrow \end{eqnarray*}$

Sei nun x* die Lösung des lin. Ausgleichsproblems. Angenommen x^* genüge nicht der den Normalengleichungen, dann gilt:

$\begin{eqnarray*} A^T(Ax^*-b) \neq 0 \Rightarrow d:=-A^T(Ax^*-b) \neq 0 \end{eqnarray*}$

Mit der Defintion $\begin{eqnarray*} x:=x^*+t\cdot d \end{eqnarray*}$ folgt dann $\begin{eqnarray*} (x-x^*)=t \cdot d \end{eqnarray*}$ :

$\begin{eqnarray*} ||Ax-b||_2^2&&= ||(Ax^*-b)+A(x-x^*)||_2^2 \\ &&= ||Ax^*-b||_2^2+2(x-x^*)^TA^T(Ax^*-b)+||A(x-x^*)||_2^2 \\&& = ||Ax^*-b||_2^2+2t\cdot d^T(-d)+||A(x-x^*)||_2^2 \\ && = ||Ax^*-b||_2^2-2t\cdot ||d||_2^2+t^2\cdot ||Ad||_2^2 \end{eqnarray*}$

Ist nun t hinreichend klein gewählt, so gilt:

$\begin{eqnarray*} t\cdot ||Ad||_2^2 \leq t \cdot ||A||_2^2\cdot ||d||_2^2 < 2||d||_2^2 \end{eqnarray*}$

Damit folgt dann:

$\begin{eqnarray*} ||Ax-b||_2^2&& < ||Ax^*-b||_2^2 \end{eqnarray*}$

Dann ist aber x^* nicht die Lösung des lin. Ausgleichsproblems.

26.07.2007, 02:15

tigerbine

Auf diesen Beitrag antworten »

2a. Existenz einer Lösung
Die Normalengleichung kann man nun dazu verwenden, die folgende Aussage zu beweisen:

$\begin{eqnarray*} \fbox{$\displaystyle \text{Das Lineare Ausgleichsproblem besitzt stets eine Lösung} $} \end{eqnarray*}$

Beweis:

Die Behauptung ist äquivlatent zu $\begin{eqnarray*} \exists~x^* \text{ mit } A^T(Ax^*-b)=0,~\quad \forall ~ A \in \mathbb R^{m \times n},~b \in \mathbb R^m \end{eqnarray*}$ . Es ergibt sich:

$\begin{eqnarray*} &&A^TAx = A^Tb \\ && \Rightarrow A^Tb \in \text{Im}(A^T) \end{eqnarray*}$

Es gilt nun $\begin{eqnarray*} \text{Im}(A^T) = \text{Im}(A^TA) \end{eqnarray*}$ . Beweis folgt im Anschluss an diesen. Somit gilt:

$\begin{eqnarray*} A^Tb \in \text{Im}(A^TA) \end{eqnarray*}$

Und schließlich die entscheidende Folgergung:

$\begin{eqnarray*} \exists~x^* \in \mathbb R^n \text {mit }A^TAx^* = A^Tb \end{eqnarray*}$

Beweis (Bildräume):

Hierzu muss man sich mit orthogonalen Räumen beschäftigen und es sind weitere Aussagen nötig. Zunächst einmal die Logikkette:

$\begin{eqnarray*} \text{Im}(A^T) = \text{Kern}(A)^\perp = \text{Kern}(A^TA)^\perp = \text{Im}(A^TA) \end{eqnarray*}$

Gleichheit 1 folgt durch Elementenvergleich.

$\begin{eqnarray*} v \in \text{Kern}(A) && \Rightarrow Av = 0 \\ && \Rightarrow (Av)^Tx = 0 ~\forall x \in \mathbb R^m \\ && \Rightarrow v^TA^Tx = 0 ~\forall x \in \mathbb R^m \\&&\Rightarrow v \in \text{Im}(A^T)^\perp \end{eqnarray*}$

$\begin{eqnarray*} \text{Im}(A^T)^\perp && \Rightarrow v^TA^Tx = 0 ~\forall~x \in \mathbb R^m \\ && \Rightarrow x^TAv = 0~\forall~x \in \mathbb R^m \\ && \Rightarrow (Av)^T(Av) = 0 \\&&\Rightarrow Av = 0 \\&& \Rightarrow v \in \text{Kern}(A) \end{eqnarray*}$

Somit gilt $\begin{eqnarray*} \text{Kern}(A) = \text{Im}(A^T)^\perp \end{eqnarray*}$ - Damit folgt auch sofort $\begin{eqnarray*} \text{Kern}(A)^\perp = \text{Im}(A^T) \end{eqnarray*}$

Nun zur zweiten Gleichheit. Es gilt:

$\begin{eqnarray*} v \in \text{Kern}(A) \Rightarrow Av = 0 \Rightarrow A^TAv = 0 \Rightarrow v \in \text{Kern}(A^TA) \end{eqnarray*}$

$\begin{eqnarray*} v \in \text{Kern}(A^TA) \Rightarrow A^TAv = 0 \Rightarrow v^TA^TAv = 0 \Rightarrow Av=0 \Rightarrow v \in \text{Kern}(A) \end{eqnarray*}$

Somit folgt dann auch die gesuchte Behauptung.

26.07.2007, 02:16

tigerbine

Auf diesen Beitrag antworten »

2b. Eindeutigkeit für Rang(A)=n
$\begin{eqnarray*} \text{Rang}(A)=n \Rightarrow A^TA \text{ ist positiv definit }\Rightarrow A^TA \text{ ist regulär } \Rightarrow \text{NGL hat eind. Lösung} \end{eqnarray*}$

(Hierauf werden die theoretischen Aussagen angewandt)

$\begin{eqnarray*} \text{Rang}(A)<n \Rightarrow A^TA \text{ ist singulär }\Rightarrow \text{NGL hat unendl. viele Lösungen} \end{eqnarray*}$

(Dieser Fall wird im WS - Beispiele betrachtet)

26.07.2007, 02:16

tigerbine

Auf diesen Beitrag antworten »

3. Lösen mit Cholesky-Zerlegung
Im Falle Rang(A)=n könnte man nun die Normalengleichung $\begin{eqnarray*} A^TAx=A^Tb \end{eqnarray*}$ aufgrund der SPD-Eigenschaft von $\begin{eqnarray*} A^TA \end{eqnarray*}$ z.B. mit dem Cholesky-Verfahren lösen.

I.A. wird man dies aufgrund der schlechten Kondition jedoch nicht tun.

26.07.2007, 02:16

tigerbine

Auf diesen Beitrag antworten »

4. QR-Zerlegungen
Betrachtet wird nun der Fall $\begin{eqnarray*} \text{Rang}(A)=n \leq m \end{eqnarray*}$ , nach (2b) existiert dann also eine eindeutige Lösung des Linearen Ausgleichsproblems.

reduzierte QR-Zerlegung

$\begin{eqnarray*} A = \hat Q \cdot \hat R,~\quad \hat Q \in \mathbb R^{m \times n},~\hat R \in \mathbb R^{n \times n} \end{eqnarray*}$

volle QR-Zerlegung

$\begin{eqnarray*} A = Q \cdot R,\quad Q = [\hat Q \tilde Q] \in \mathbb R^{m \times n},~R = \begin{bmatrix} \hat R \\0 \end{bmatrix} \in \mathbb R^{m \times n} \end{eqnarray*}$

$\begin{eqnarray*} \Rightarrow A = QR = [\hat Q \tilde Q]\begin{bmatrix} \hat R \\0 \end{bmatrix} = \hat Q \hat R \end{eqnarray*}$

Dabei sind die Spalten von $\begin{eqnarray*} \hat Q \end{eqnarray*}$ orthonormiert und $\begin{eqnarray*} \hat R \end{eqnarray*}$ eine obere Dreiecksmatrix.

26.07.2007, 02:18

tigerbine

Auf diesen Beitrag antworten »

4a. Eindeutigkeit der reduzierten QR-Zerlegung
Sollte nun eine reduzierte QR-Zerlegung exisitieren, so ist sie im Wesentlichen eindeutig. Betrachtet man also 2 reduzierte QR-Zerlegungen, so existiert eine Diagonalmatrix D (mit Werten +/- 1), so dass gilt:

$\begin{eqnarray*} A = \hat Q_1 \hat R_1,~A = \hat Q_2 \hat R_2 \end{eqnarray*}$

$\begin{eqnarray*} \hat Q_1 = \hat Q_2D,~\hat R_2 = D \hat R_1 \end{eqnarray*}$

Beweis:

Es nun nach Voraussetzung $\begin{eqnarray*} \hat Q_1 \hat R_1 = \hat Q_2 \hat R_2 \end{eqnarray*}$ . Da nun die Spalten der Qs orthonormiert sind, gilt:

$\begin{eqnarray*} \hat Q_1^T \hat Q_1 = I_n,\quad \hat Q_2^T \hat Q_2 = I_n \end{eqnarray*}$

Damit erhält man dann aus der ersten Gleichheit:

$\begin{eqnarray*} \hat Q_2^T\hat Q_1 \hat R_1 &&= \hat Q_2^T\hat Q_2 \hat R_2 \\ \hat Q_2^T\hat Q_1 \hat R_1 &&= \hat R_2 \end{eqnarray*}$

und

$\begin{eqnarray*} \hat Q_1^T\hat Q_2 \hat R_2 &&= \hat R_1 \end{eqnarray*}$

Es sind nun die Rs regulär, sonst wäre Rang(A) <n, gilt:

$\begin{eqnarray*} \hat Q_2^T\hat Q_1 &&= \hat R_2 \hat R_1^{-1} \end{eqnarray*}$

$\begin{eqnarray*} \hat Q_1^T\hat Q_2 &&= \hat R_1 \hat R_2^{-1} \end{eqnarray*}$

Nun definiert man $\begin{eqnarray*} D:=\hat R_2 \hat R_1^{-1} \end{eqnarray*}$ und weist nach, dass es sich um eine orthogonale Diagonalmatrix handelt.

Mit dem Wissen über Dreiecksmatrizen folgt, dass auch die Inversen obere Dreiecksmatrizen sind. Es ist also auch D eine obere Dreiecksmatrix und ihre Transponierte somit eine untere Dreiecksmatrix.Es gilt weiter:

$\begin{eqnarray*} D^T = (\hat R_2 \hat R_1^{-1})^T = (\hat Q_2^T\hat Q_1)^T = \hat Q_1^T \hat Q_2 = \hat R_1 \hat R_2^{-1} \end{eqnarray*}$

Somit muss D eine Diagonalmatrix sein. Nun folgt weiter

$\begin{eqnarray*} D^T = \hat R_1 \hat R_2^{-1} = (\hat R_2 \hat R_1^{-1})^{-1} = D^{-1} \end{eqnarray*}$

Womit nun auch die Orthogonalität von d gezeigt wurde, d.h. die Diagonalwerte sind +/- 1.

Mit $\begin{eqnarray*} \hat Q_1 = \hat Q_2\hat R_2 \hat R_1^{-1} = \hat Q_2D \end{eqnarray*}$ folgt dann die Behauptung.

Im folgenden wird die Existenz der red. QR-Zerlegung erst einmal angenommen. Der Beweis folgt hier.

26.07.2007, 02:18

tigerbine

Auf diesen Beitrag antworten »

4b. Lösen der Normalengleichung - red. QR
Sei also Rang (A) = n und es existiere eine red. QR-Zerlegung. Dann folgt:

$\begin{eqnarray*} A^TAx &&= A^Tb \\ (\hat Q \hat R)^T(\hat Q \hat R) &&= (\hat Q \hat R)^Tb \\ \hat R^T\hat Q^T \hat Q \hat R x &&= (\hat Q \hat R)^Tb \\ \hat R^T \hat R x &&= \hat R^T \hat Q^T b\\ \hat Rx && = \hat Q^Tb \end{eqnarray*}$

Wie kann man nun also die eindeutige Lösung x* bestimmen?

Bestimme die red. QR-Zerlegung
$\begin{eqnarray*} \hat c:=\hat Q^T b \end{eqnarray*}$
Löse $\begin{eqnarray*} \hat Rx=\hat c \end{eqnarray*}$ durch Rückwärtssubstitution

26.07.2007, 02:18

tigerbine

Auf diesen Beitrag antworten »

4c. Lösen des Ausgleichsproblems - volle QR
$\begin{eqnarray*} A = Q \cdot R,\quad Q = [\hat Q \tilde Q] \in \mathbb R^{m \times n},~R = \begin{bmatrix} \hat R \\0 \end{bmatrix} \in \mathbb R^{m \times n},~c:=Q^Tb = \begin{pmatrix} \hat c \\ \tilde c\end{pmatrix} \end{eqnarray*}$

$\begin{eqnarray*} ||Ax-b||_2 && = ||QRx - b||_2 \\&& = ||QRx - QQ^Tb||_2 \\&&= ||Q(Rx-Q^Tb)||_2 \\ && = \underbrace{||Q||_2}_{=1} \cdot ||Rx-Q^Tb||_2 \\&& = ||Rx-c||_2 \\ &&= || \begin{bmatrix} \hat R \\0 \end{bmatrix} x - \begin{pmatrix} \hat c \\ \tilde c\end{pmatrix}||_2 \end{eqnarray*}$

$\begin{eqnarray*} \Rightarrow \min || \begin{bmatrix} \hat R \\0 \end{bmatrix} x - \begin{pmatrix} \hat c \\ \tilde c\end{pmatrix}||_2 \end{eqnarray*}$

$\begin{eqnarray*} \Rightarrow \min || \begin{bmatrix} \hat R \\0 \end{bmatrix} x - \begin{pmatrix} \hat c \\ \tilde c\end{pmatrix}||_2^2 &&= ||\begin{pmatrix} \hat Rx - \hat c \\ 0 - \tilde c \end{pmatrix} ||_2^2 \\&& = ||\hat R x- \hat c||_2^2 + ||\tilde c||_2^2 \end{eqnarray*}$

Da der zweite Summand eine Konstante ist folgt:

$\begin{eqnarray*} \Rightarrow \min ||\hat R x- \hat c||_2^2 \end{eqnarray*}$

Dies wird offensichtlich für $\begin{eqnarray*} \hat Rx = \hat c \end{eqnarray*}$ minimal und wir sind wieder bei der Lösung mittels reduzierten QR-Zerlegung aus (4b).

Lehrer

Für m = n können QR-Zerlegungen auch zur Lösung von Linearen Gleichungssystemen verwendet werden.

Wie man nun eine QR-Zerlegung bestimmt wird in den Abschnitten (5),(6),(7) gezeigt.

26.07.2007, 02:19

tigerbine

Auf diesen Beitrag antworten »

5. Gram-Schmidt-Orthogonalisierung
Gesucht ist nun die red. QR-Zerlegung. Zunächst notiert man:

$\begin{eqnarray*} A = \hat Q \hat R \Leftrightarrow \begin{pmatrix} \vdots & & \vdots \\ a^1& \hdots & a^n \\ \vdots&&\vdots \end{pmatrix} =\begin{pmatrix} \vdots & & \vdots \\ q^1& \hdots & q^n \\ \vdots&&\vdots \end{pmatrix}\begin{pmatrix} r_{11} & &r_{1n} \\ & \ddots & \vdots \\ &&r_{nn} \end{pmatrix} \end{eqnarray*}$

Lehrer

Die Potenzen Kennzeichen hier die Spaltenvektoren

ausgeschrieben bedeutet das:

$\begin{eqnarray*} a^1=r_{11}\cdot q^1 \end{eqnarray*}$

$\begin{eqnarray*} a^2=r_{12}\cdot q^1 + r_{22} \cdot q^2 \end{eqnarray*}$

$\begin{eqnarray*} a^3=r_{13}\cdot q^1 + r_{23} \cdot q^2 + r_{33} \cdot q^3 \end{eqnarray*}$

$\begin{eqnarray*} \vdots \end{eqnarray*}$

$\begin{eqnarray*} a^n=r_{1n}\cdot q^1 + r_{2n}\cdot q^2 + \hdots + r_{nn} \cdot q^n \end{eqnarray*}$

allgemein:

$\begin{eqnarray*} a^j = \sum_{i=1}^j~r_{ij} \cdot q^i \quad j=1,...,n \qquad (*) \end{eqnarray*}$

Damit erhält man nun eine Berechnungsvorschrift für die orthonormalen Vektoren "q" und die zugehörigen Koeffizienten "r". Für j=1 ergibt sich:

$\begin{eqnarray*} ||a^1||_2 = |r_{11}| \cdot ||q^1||_2 \end{eqnarray*}$

Da Q orthogonal ist, gilt

$\begin{eqnarray*} ||q^{j}||_2 = 1 ,~j=1,...,n \end{eqnarray*}$

$\begin{eqnarray*} \Rightarrow r_{11} = \pm ||a^1||_2,\quad q^1 = \frac{1}{r_{11}}\cdot a^1 \end{eqnarray*}$

Um nun eine eindeutige Lösung zu erhalten, sollen die Diagonalelemente von R positiv sein. Weiter erhält man durch Umstellen zunächst:

$\begin{eqnarray*} r_{22} = ||a^2-r_{12}\cdot q^1||_2 \quad \text{Normierung von } q^2 \end{eqnarray*}$

$\begin{eqnarray*} q^2 = \frac{1}{r_{22}}\cdot (a^2-r_{12}\cdot q^1)\quad \text{Umstellen von (*)} \end{eqnarray*}$

Nun die Orthogonalitätsbedingung:

$\begin{eqnarray*} 0 \stackrel{!}=(q^1)^Tq^2 \end{eqnarray*}$

$\begin{eqnarray*} 0 \stackrel{!}=(q^1)^T\frac{1}{r_{22}} \cdot (a^2-r_{12} \cdot q^1) \end{eqnarray*}$

$\begin{eqnarray*} 0 \stackrel{!}=\frac{1}{r_{22}} \cdot [(q^1)^Ta^2-r_{12} \cdot (q^1)^Tq^1] \end{eqnarray*}$

$\begin{eqnarray*} 0 \stackrel{!}=\frac{1}{r_{22}} \cdot [(q^1)^Ta^2-r_{12} ] \end{eqnarray*}$

$\begin{eqnarray*} r_{12} = (q^1)^Ta^2 \end{eqnarray*}$

Allgemein ergibt sich dann:

$\begin{eqnarray*} r_{jj} = ||a^j-\sum_{i=1}^{j-1}~r_{ij} \cdot q^{i}||_2 \end{eqnarray*}$

$\begin{eqnarray*} q^j:=\frac{1}{r_{jj}} \cdot \Bigl( a^j - \sum_{k=1}^{j-1}~r_{kj} \cdot q^k \Bigr) \end{eqnarray*}$

$\begin{eqnarray*} r_{ij} = (q^{i})^Ta^j~\quad \forall ~ i=1,...,j-1 \end{eqnarray*}$

26.07.2007, 02:20

tigerbine

Auf diesen Beitrag antworten »

5a. Berechnung der red. QR-Zerlegung
Fasst man nun die obige Konstruktion in einem Algorithmus zusammen, so erhält man das

Gram-Schmidtsche Orthogonalisierungsverfahren

26.07.2007, 02:21

tigerbine

Auf diesen Beitrag antworten »

5b. Existenz der red. QR-Zerlegung
Mit dieser Vorschrift (Gram-Schmidt-Verfahren) ist die Existenz einer red. QR-Zerlegung konstruktiv (fast) bewiesen worden. Es ist noch zu zeigen, dass das Verfahren wohldefiniert ist. Dies ist dann der Fall, wenn für die Nenner gilt $\begin{eqnarray*} r_{jj} \neq 0 \end{eqnarray*}$ .

Wäre nun $\begin{eqnarray*} r_{jj} = 0 \end{eqnarray*}$ für ein kleinstes $\begin{eqnarray*} j \in \{1,...,n \} \end{eqnarray*}$ , so folgt aus (*):

$\begin{eqnarray*} a^j \in \text{span} \{q^1,...,q^{j-1}\} = \text{span} \{a^1,...,a^{j-1}\} \end{eqnarray*}$

Dies steht aber im Widerspruch zu Rang(A) = n.

26.07.2007, 02:21

tigerbine

Auf diesen Beitrag antworten »

5c. Implementierung - Modifiziertes Gram-Schmidt Verfahren
Das Verfahren ist jedoch numerisch instabil, da aufgrund der Rechengenauigkeit die Orthogonalität der Vektoren schnell verloren geht. Ein Rechenbeispiel findet sich dazu im WS - Beispiele.

Zunächst einmal führt man einen Hilfsvektor $\begin{eqnarray*} p^j \end{eqnarray*}$ ein, um das Überschreiben der Vektoren $\begin{eqnarray*} q^j \end{eqnarray*}$ im Algorithmus zu vermeiden (siehe auch Rechenbeispiel)

Ausgeschrieben lautet die Berechnung der nicht normierten Vektoren

$\begin{eqnarray*} p^j && =a^j - \sum_{i=1}^{j-1}~r_{ij}\cdot q^i = a^j-\sum_{i=1}^{j-1}~[(q^i)^Ta^j]q^i = a^j-\sum_{i=1}^{j-1}~[q^i(q^i)^Ta^j] \\&&= \Bigl(I_n - \sum_{i=1}^{j-1} q^i(q^i)^T \Bigr) a^j \end{eqnarray*}$

Wegen der Orthogonalität der Vektoren q folgt (Beweis am Ende des Posts):

$\begin{eqnarray*} \Bigl(I_n - \sum_{i=1}^{j-1} q^i(q^i)^T \Bigr) a^j &&= \prod_{i=j-1}^1~\Bigl(I_n - \sum_{i=1}^{j-1} q^i(q^i)^T \Bigr)~a^j \\&&=\Bigl(I_n - q^{j-1}(q^{j-1})^T \Bigr) ... \Bigl(I_n -q^1(q^1)^T \Bigr)~a^j \end{eqnarray*}$

Daher kann man die Vektoren p auch wie folgt berechnen:

Es gilt dann:

$\begin{eqnarray*} p^{j,1}=a^j, \end{eqnarray*}$

$\begin{eqnarray*} p^{j,2}=\Bigl(I_n -q^1(q^1)^T \Bigr) ~a^j, \end{eqnarray*}$

$\begin{eqnarray*} \vdots \end{eqnarray*}$

$\begin{eqnarray*} p^{j,j}=\Bigl(I_n -q^{j-1}(q^{j-1})^T \Bigr)...\Bigl(I_n -q^1(q^1)^T \Bigr) ~ ~a^j \end{eqnarray*}$

Speichert man die Hilfsvektoren $\begin{eqnarray*} p^{j,i} \end{eqnarray*}$ nicht explizit ab und überschriebt stattdessen einen Vektor $\begin{eqnarray*} q^{i} \end{eqnarray*}$ mit den entsprechenden Einträgen, so erhält man das modifizierte Gram-Schmidt-Verfahren.

Beweis:

Ist durch Induktion zu führen. Wir verifizieren hier nur den Induktionsanfang. Für n=2 gilt, aufgrund der Orthogonalität der Vektoren

$\begin{eqnarray*} <q^1,q^2> = q_1^1q_1^2 + q_2^1q_2^2 = 0 \Leftrightarrow -q_1^1q_1^2 = q_2^1q_2^2 \end{eqnarray*}$

$\begin{eqnarray*} &&\begin{pmatrix}1-q_1^1q_1^1 & -q_1^1q_2^1 \\ -q_1^1q_2^1 & 1 - q_2^1q_2^1 \end{pmatrix}\begin{pmatrix} 1-q_1^2q_1^2 & -q_1^2q_2^2 \\ -q_1^2q_2^2 & 1 - q_2^2q_2^2 \end{pmatrix} =\\ \\&& \begin{pmatrix} 1-(q_1^2)^2 - (q_1^2)^2 + (q_1^1)^2(q_1^2)^2 + q_1^1q_1^2q_2^1q_2^2 &-q_1^1q_2^1+q_1^1q_2^1q_2^2q_2^2 - q_1^1q_2^1+q_1^1q_1^1q_1^2q_2^2 \\-q_1^1q_2^1 + q_1^1q_2^1q_1^2q_1^2 - q_1^2q_2^2+q_2^1q_2^1q_1^2q_2^2 & q_1^1q_2^1q_1^2q_2^2+1-q_2^1q_2^1-q_2^2q_2^2+q_2^1q_2^1q_2^2q_2^2 \end{pmatrix} =\\ \\ && \begin{pmatrix} 1-(q_1^2)^2 - (q_1^2)^2 &-q_1^1q_2^1-q_1^2q_2^2 \\ -q_1^1q_2^1-q_1^2q_2^2& 1-q_2^1q_2^1 - q_2^2q_2^2\end{pmatrix} \end{eqnarray*}$

Somit ist der Anfang gemacht. Viel Spass beim Rest der Indize-Schlacht Augenzwinkern

26.07.2007, 02:21

tigerbine

Auf diesen Beitrag antworten »

6. Householder Spiegelungen
Nun soll eine volle QR-Zerlegung bestimmt werden. Wieder hat A den Rang n. Es ist dann

$\begin{eqnarray*} A = QR \Leftrightarrow Q^{-1}A = R \Leftrightarrow Q^{T}A = R \end{eqnarray*}$

Ähnlich dem Prinzip beim Gaußalgorithmus ist es nun das Ziel die Matrix $\begin{eqnarray*} Q^T \end{eqnarray*}$ zu konstruieren, so dass $\begin{eqnarray*} Q^TA \end{eqnarray*}$ eine obere Dreiecksmatrix ist. Man schreibt:

$\begin{eqnarray*} Q^T = Q_l \cdot .... \cdot Q_1,\quad l:=\min\{m-1,n\} \end{eqnarray*}$

Dabei annulliert $\begin{eqnarray*} Q_j \end{eqnarray*}$ die Eintrage der j-ten Spalte unterhalb der Diagonalen.

26.07.2007, 02:23

tigerbine

Auf diesen Beitrag antworten »

6a. Eigenschaften
Householder-Spiegelung (Definition)

$\begin{eqnarray*} H:=I_n - 2\cdot \underbrace{uu^T}_{Dyade}\quad \text{mit } u \in \mathbb R^n,~||u||_2=1 \end{eqnarray*}$

Konstruktion von u

x sei nun der Teil des Spaltenvektors von "A", den es zu annullieren gilt, zzgl. des Diagonalelements.

$\begin{eqnarray*} u:=x + \text{sign}(x_1)\cdot ||x||_2 \cdot e_1 \quad \text{(Bedingungen für die zu erzeugenden Nullen)} \end{eqnarray*}$

$\begin{eqnarray*} \beta:=\frac{2}{u^Tu} \quad \text{(Normierung)} \end{eqnarray*}$

$\begin{eqnarray*} H=I-\beta\cdot uu^T \end{eqnarray*}$

Eigenschaften der Matrix H

H ist symmetrisch, dies folgt direkt aus der Symmetrie von I und der Dyade.
H ist orthogonal, denn
$\begin{eqnarray*} H^TH = (I_n-2uu^T)^T(I_n-2uu^T) = I - 4uu^T + 4u\underbrace{(u^Tu)}_{1}u^T = I_n \end{eqnarray*}$
u ist ein Eigenvektor von H, denn
$\begin{eqnarray*} Hu = (I_n-2uu^T) u = u - 2u(u^Tu) = (-1) \cdot u \end{eqnarray*}$
Der Eigenraum von (-1) ist eindimensional.
H hat sonst nur noch den Eigenwert 1. Die Dimension seines Eigenraums ist (n-1). Da H symmetrisch und regulär ist, gibt es eine ONB aus Eigenvektoren von H, so dass H bzgl. dieser Diagonalagestalt hat .Sei v ein zu u orthogonaler Vektor. Dann gilt:
$\begin{eqnarray*} Hv = v - 2u\underbrace{u^Tv}_{0} = v \end{eqnarray*}$
H ist eine Spiegelung, denn
$\begin{eqnarray*} \det(H) = -1 \end{eqnarray*}$
Hx ist ein Vielfaches des ersten Einheitsvektors mit dem Faktor $\begin{eqnarray*} \alpha = - \text{sign}(x_1)||x||_2 \end{eqnarray*}$

Beweis:

Die Vektoren x und u unterscheiden sich nur im ersten Eintrag, d.h. es gilt $\begin{eqnarray*} u_j=x_j,~j=2,...,n \end{eqnarray*}$ . So kann man Dyade wie folgt schreiben:

$\begin{eqnarray*} U = \begin{pmatrix}u_1u_1 &u_1x_2 & \hdots & u_1x_n \\ x_2u_1 &x_2x_2 & \hdots & x_2x_n \\ \vdots & \vdots & \ddots & \vdots \\ x_nu_1 & x_nx_2 & \hdots & x_nx_n \end{pmatrix} \end{eqnarray*}$

Nun ist

$\begin{eqnarray*} Hx = (I-\beta uu^T)x = x - \frac{2}{u^Tu}\cdot Ux \end{eqnarray*}$ .

Für den Nenner des Faktors ß gilt:

$\begin{eqnarray*} u^Tu &&= \sum_{k=1}^n~u_k^2 = u_1^2 + \sum_{k=2}^n~x_k^2 \\ &&=\Bigl(x_1+\text{sign}(x_1)||x||_2\Bigr)^2 + \sum_{k=2}^n~x_k^2 \\&&= x_1^2 + 2x_1\cdot \text{sign}(x_1)||x||_2 +(\text{sign}(x_1)||x||_2)^2 + \sum_{k=2}^n~x_k^2 \\&& = 2\cdot \Bigl( ||x||_2^2+x_1 \cdot \text{sign}(x_1)||x||_2 \Bigr) \end{eqnarray*}$

Für den j-ten Eintrag des Vektors Ux gilt:

$\begin{eqnarray*} (Ux)_j = \sum_{k=1}^{n}~(u_ju_k)x_k \end{eqnarray*}$

Für j=1 gilt:

$\begin{eqnarray*} (Ux)_1&& = \sum_{k=1}^n~(u_1u_k)x_k = u_1 \cdot \Bigl(u_1x_1 + \sum_{k=2}^n~x_kx_k \Bigr)\\ &&=u_1 \cdot \Bigl((x_1 + \text{sign}(x_1)\cdot ||x||_2)x_1 +\sum_{k=2}^n~x_kx_k \Bigr) \\&& =u_1\cdot \Bigl( ||x||_2^2 + x_1 \cdot \text{sign}(x_1)||x||_2\Bigr) \end{eqnarray*}$

$\begin{eqnarray*} (Hx)_1&&=[(I-\beta\cdot uu^T)x]_1 = x_1-\beta\cdot (Ux)_1 \\&& = x_1 - \beta\cdot u_1\cdot \Bigl( ||x||_2^2 + x_1 \cdot \text{sign}(x_1)||x||_2\Bigr) \\&& = x_1 - u_1\\&& = -\text{sign}(x_1)||x||_2 \end{eqnarray*}$

Für j =2,...,n gilt:

$\begin{eqnarray*} (Ux)_1&& = \sum_{k=1}^n~(x_ju_k)x_k = \\&&=x_ju_1x_1 + \sum_{k=2}^n~(x_jx_k)x_k \\&&=x_ju_1x_1 + x_j\cdot \sum_{k=2}^n~(x_k)^2 = x_j\cdot (u_1x_1 + \sum_{k=2}^n~(x_k)^2) \\&&= x_j\cdot (x_1^2+x_1\cdot \text{sign}(x_1)||x||_2 + \sum_{k=2}^n~(x_k)^2) \\&& = x_j\cdot (||x||_2^2 + x_1\cdot \text{sign}(x_1)||x||_2) \end{eqnarray*}$

$\begin{eqnarray*} (Hx)_j=x_j -x_j = 0 \end{eqnarray*}$

Somit ist Hx ein Vielfaches des ersten Einheitsvektors mit obigem Faktor.

26.07.2007, 02:23

tigerbine

Auf diesen Beitrag antworten »

6b. Konstruktion von Q
Wie hängen nun die eingangs erwähnten Matrizen $\begin{eqnarray*} Q_i \end{eqnarray*}$ und die Matrizen H zusammen? Mit entsprechenden Dimensionen gilt:

$\begin{eqnarray*} Q_i = \begin{pmatrix}\underbrace{I}_{(i-1) \times (i-1)} &\underbrace{0}_{(i-1) \times (n+1-i)}\\ \underbrace{0}_{(n+1-i) \times (i-1)}&\underbrace{H_i}_{(n+1-i)\times (n+1-i)} \end{pmatrix} \end{eqnarray*}$

Auch die Matrizen $\begin{eqnarray*} Q_i \end{eqnarray*}$ sind symmetrisch. Es gilt:

$\begin{eqnarray*} Q_i^T&&= \begin{pmatrix}\underbrace{I}_{(i-1) \times (i-1)} &\underbrace{0}_{(i-1) \times (n+1-i)}\\ \underbrace{0}_{(n+1-i) \times (i-1)}&\underbrace{H_i^T}_{(n+1-i)\times (n+1-i)} \end{pmatrix} \\&&=\begin{pmatrix}\underbrace{I}_{(i-1) \times (i-1)} &\underbrace{0}_{(i-1) \times (n+1-i)}\\ \underbrace{0}_{(n+1-i) \times (i-1)}&\underbrace{H_i}_{(n+1-i)\times (n+1-i)} \end{pmatrix} \\&&=Q_i \end{eqnarray*}$

Ebenso sind sie orthogonal:

$\begin{eqnarray*} Q_i^TQ_i &&= \begin{pmatrix}\underbrace{I^TI}_{(i-1) \times (i-1)} &\underbrace{0}_{(i-1) \times (n+1-i)}\\ \underbrace{0}_{(n+1-i) \times (i-1)}&\underbrace{H_i^TH_i}_{(n+1-i)\times (n+1-i)} \end{pmatrix} \\&&= I \end{eqnarray*}$

Und mit der Rechenregel für Blockdiagonalmatrizen folgt auch sofort:

$\begin{eqnarray*} \det(Q_i)=-1 \end{eqnarray*}$

Handelt es sich auch um eine Householder-Matrix? Dazu muss man nur einen Vektor $\begin{eqnarray*} \hat u_i \end{eqnarray*}$ finden mit:

$\begin{eqnarray*} Q_i = I - \frac{2}{\hat u_i^T \hat u_i} \cdot \hat u_i \hat u_i^T \end{eqnarray*}$

Wählt man nun:

$\begin{eqnarray*} (\hat u_i)^T = ( \underbrace{0}_{(i-1)},~\underbrace{u_i}_{(n+1-i)})^T \end{eqnarray*}$

Dann ist:

$\begin{eqnarray*} \frac{2}{\hat u_i^T \hat u_i} = \frac{2}{u_i^T u_i} \end{eqnarray*}$

$\begin{eqnarray*} \hat U_i = \begin{pmatrix} \underbrace{0}_{(i-1) \times (i-1)} & \underbrace{0}_{(i-1) \times (n+1-i)}\\ \underbrace{0}_{(n+1-i) \times (i-1)} &\underbrace{U_i}_{(n+1-i) \times (n+1-i)} \end{pmatrix} \end{eqnarray*}$

Damit ergibt:

$\begin{eqnarray*} I_n - \beta \cdot \hat U_i &&= \begin{pmatrix} \underbrace{I_{i-1}}_{(i-1) \times (i-1)} & \underbrace{0}_{(i-1) \times (n+1-i)}\\ \underbrace{0}_{(n+1-i) \times (i-1)} &\underbrace{I_{n+1-i}-\beta \cdot U_i}_{(n+1-i) \times (n+1-i)} \end{pmatrix} =\begin{pmatrix} \underbrace{I_{i-1}}_{(i-1) \times (i-1)} & \underbrace{0}_{(i-1) \times (n+1-i)}\\ \underbrace{0}_{(n+1-i) \times (i-1)} &\underbrace{H_i}_{(n+1-i) \times (n+1-i)} \end{pmatrix} \\ && = Q_i \end{eqnarray*}$

Somit sind die $\begin{eqnarray*} Q_i \end{eqnarray*}$ auch Householder-Matrizen.

26.07.2007, 02:23

tigerbine

Auf diesen Beitrag antworten »

6c. Implementierung
Wie im Beispiel Workshop angemerkt, muss hier nicht die Matrix Q explizit berechnet werden. Desweiteren reicht es sogar aus, nur die zur Konstruktion benötigten Vektoren u abzuspeichern. Im folgenden wird die Matrix A mit diesen Überschrieben. Im oberen Dreieck wird die Matrix R gespeichert. Ihre Diagonale speichert man im Vektor p.

26.07.2007, 02:23

tigerbine

Auf diesen Beitrag antworten »

6d. Geometrische Deutung
Wir hatten festgestellt, dass gilt:

Hx ist ein Vielfaches des ersten Einheitsvektors mit dem Faktor

Zitat:

$\begin{eqnarray*} \alpha = - \text{sign}(x_1)||x||_2 \end{eqnarray*}$

Die Konstruktion von u mittels Signum hat numersische Gründe. Es soll hier

Zitat:

$\begin{eqnarray*} u:=x + \text{sign}(x_1)\cdot ||x||_2 \cdot e_1 \quad \text{(Bedingungen für die zu erzeugenden Nullen)} \end{eqnarray*}$

eine Auslöschung der ersten Komponente von u vermieden werden. Theoretisch hätte man jedoch 2 Möglichkeiten, denn Vektor x (In der Skizze mit v benannt) abzubilden.

Links

$\begin{eqnarray*} x \to ||x||e_1 \end{eqnarray*}$

rechts

$\begin{eqnarray*} x \to -||x||e_1 \end{eqnarray*}$

26.07.2007, 02:24

tigerbine

Auf diesen Beitrag antworten »

7. Givens-Rotationen
Nun soll eine volle QR-Zerlegung bestimmt werden. Wieder hat A den Rang n. Es ist dann

$\begin{eqnarray*} A = QR \Leftrightarrow Q^{-1}A = R \Leftrightarrow Q^{T}A = R \end{eqnarray*}$

Ähnlich dem Prinzip beim Gaußalgorithmus ist es nun das Ziel die Matrix $\begin{eqnarray*} Q^T \end{eqnarray*}$ zu konstruieren, so dass $\begin{eqnarray*} Q^TA \end{eqnarray*}$ eine obere Dreiecksmatrix ist. Man schreibt:

$\begin{eqnarray*} Q^T = G_{ml} \cdot .... \cdot G_{21},\quad l:=\min\{m-1,n\} \end{eqnarray*}$

Dabei annulliert $\begin{eqnarray*} G_ij \end{eqnarray*}$ den Eintrag an der Stelle (i,j).

26.07.2007, 02:24

tigerbine

Auf diesen Beitrag antworten »

7a. Eigenschaften
Givens-Rotation (Definition)

$\begin{eqnarray*} G_{ij}:=\begin{pmatrix}I_{j-1} & & & \\ &c&0&\hdots &0&s \\ &0&1&&0&0\\& \vdots & \vdots & \ddots &\vdots &\vdots \\&0&0&&1&0\\ &-s&0&\hdots &0&c \\ &&&&&&I_{n-i} \end{pmatrix} \end{eqnarray*}$

Konstruktion von c und s

$\begin{eqnarray*} p:=\pm \sqrt{a_{jj}^2+a_{ij}^2} \end{eqnarray*}$

$\begin{eqnarray*} c:=\frac{a_{jj}}{p},~s:=\frac{a_{ij}}{p} \end{eqnarray*}$

Eigenschaften der Matrix G

G ist offensichtlich schiefsymmetrisch
$\begin{eqnarray*} c^2+s^2 = 1 \end{eqnarray*}$ , denn $\begin{eqnarray*} c^2+s^2 = \frac{a_{jj}^2 + a_{ij}^2}{p^2} = 1 \end{eqnarray*}$
G ist orthogonal, denn
$\begin{eqnarray*} (G_{ij}^TG_{ij})_{ij} &&=(0,...,0,c,0,...,-s,0,...,0)(0,...,0,s,0,...,c,0,...,0)^T \\&&=0 \end{eqnarray*}$

$\begin{eqnarray*} (G_{ij}^TG_{ij})_{ji} &&=(0,...,0,s,0,...,c,0,...,0)(0,...,0,c,0,...,-s,0,...,0)^T \\&&=0 \end{eqnarray*}$

$\begin{eqnarray*} ((G_{ij}^TG_{ij})_{ii} &&=(0,...,0,c,0,...,-s,0,...,0)(0,...,0,c,0,...,-s,0,...,0)^T \\&&=1 \end{eqnarray*}$

$\begin{eqnarray*} ((G_{ij}^TG_{ij})_{jj} &&=(0,...,0,s,0,...,c,0,...,0)(0,...,0,s,0,...,c,0,...,0)^T \\&&=1 \end{eqnarray*}$

Somit ist offensichtlich $\begin{eqnarray*} G_{ij}^TG_{ij} = I \end{eqnarray*}$
G ist eine Drehung, denn mit dem Entwicklungssatz von Laplace und der Determinantenregel für Blockdiagonalmatrizen folgt (mit einer kleinen Fallunterscheidung)
$\begin{eqnarray*} \det(G_{ij}) = 1 \cdot (c^2+s^2) \cdot 1 = 1 \end{eqnarray*}$
G besitzt den Eigenwert 1 und der zugehörige Eigenraum hat die Dimension (n-2). Eigenvektoren sind die (n-2) entsprechenden Standardeinheitsvektoren.
Die restlichen 2 Eigenwerte sind die (meist) komplexen Nullstellen des Polynoms $\begin{eqnarray*} \chi(\lambda) = \lambda^2-2c\cdot \lambda + 1 \end{eqnarray*}$
$\begin{eqnarray*} G_{ij} \end{eqnarray*}$ annulliert $\begin{eqnarray*} a_{ij} \end{eqnarray*}$ , denn (j < i)
$\begin{eqnarray*} (G_{ij}A)_{ij} &&= (0,...,0,-s,0,...,0,c,0,...,0)(a_{1j},...,a_{jj},...a_{ij},...,a_{nj})^T \\&& = -s\cdot a_{jj} + c \cdot a_{ij} = 0 \end{eqnarray*}$

26.07.2007, 02:25

tigerbine

Auf diesen Beitrag antworten »

7b. Konstruktion von Q
Da die Menge der Drehungen eine Gruppe, die spezielle Orthogonale Gruppe bildet, ist auch Q eine Drehung. Eine Givens-Rotation wird Q i.A. jedoch nicht mehr sein.

26.07.2007, 02:25

tigerbine

Auf diesen Beitrag antworten »

7c. Implementierung

26.07.2007, 02:25

tigerbine

Auf diesen Beitrag antworten »

7d. Geometrische Deutung
Die Wahl der Variablen c und s erklärt sich mit der folgenden Skizze. Der Vekotr v wird in der durch die Einheitsvektoren $\begin{eqnarray*} e_j,e_i \end{eqnarray*}$ aufgespannten Ebene um dem Winkel $\begin{eqnarray*} \varphi \end{eqnarray*}$ gedreht. Dabei gilt:

$\begin{eqnarray*} c=\cos(\varphi)~s=\sin(\varphi) \end{eqnarray*}$

26.07.2007, 02:26

tigerbine

Auf diesen Beitrag antworten »

8. Gauß'sches Ausgleichsverfahren (Methode der kleinsten Quadrate)
Hier möchte ich auf den Wikipedia Artikel verweisen:

http://de.wikipedia.org/wiki/Methode_der_kleinsten_Quadrate

Ähnlich wie im Workshop "Polynominterpolation" ist hier en Datengitter gegeben. Z.B. aus einer Messreihe.

jedoch sucht man nun nicht eine interpolierende Funktion (je mehr Knoten desto höher der Grad der Funktion), sondern eine Modellfunktion $\begin{eqnarray*} y_m \end{eqnarray*}$ best. Gestalt soll so durch die Datenwolke gelegt werden, dass die Quadratsumme der senkrechten Abweichungen:

$\begin{eqnarray*} e_k:=y_m(t_k)-y_k~\forall k=1,...n \end{eqnarray*}$

minimiert wird.

Dabei wirs $\begin{eqnarray*} y_m \end{eqnarray*}$ je nach zugrunde liegendem Problem unterschiedlich bestimmt. Allgemein formuliert sich das Minimierungsproblem wie folgt:

$\begin{eqnarray*} \min_{x}~\sum_{i=1}^n(y_m(x)-y)^2 \end{eqnarray*}$

Dies ist äquivalent zu Minimierung der euklidischen Norm des Differenzenvektors:

$\begin{eqnarray*} \min_{x}||y_m(x)-y||_2 \end{eqnarray*}$

26.07.2007, 02:26

tigerbine

Auf diesen Beitrag antworten »

8a. Lineare Modellfunktion
Machen wir uns einmal die erwähnte Äquivlanz bewußt. Gesucht ist eine lineare Modellfunktion $\begin{eqnarray*} y_m(t)=x_0 + x_1 \cdot t \end{eqnarray*}$ . Das Optimierungsproblem lautet nun:

$\begin{eqnarray*} &&\min_{x} \sum_{i=1}^n~\Bigl( y_m(t_i) - y_i\Bigr)^2 \\ \Leftrightarrow && \min_{x} \sum_{i=1}^n \Bigl( x_1 \cdot 1 + x_1 \cdot t_i - y_i)^2 \Bigr) \\\Leftrightarrow && \min _{x}\sum_{i=1}^n~(z_i)^2 ,\qquad z_i:=x_0 + x_1\cdot t_i -y_i \\\Leftrightarrow && \min_x (z^Tz),~\qquad z:=(z_1,...,z_n)^T\\ \Leftrightarrow && \min_{x}~||z||_2^2 \\ \Leftrightarrow && \min_{x}~||z||_2 \\ \Leftrightarrow && \min_{x}~\begin{Vmatrix} \begin{pmatrix}x_0\cdot 1 + x_1 \cdot t_1 -y_1\\ \vdots \\x_0\cdot 1 + x_1 \cdot t_n -y_n \end{pmatrix} \end{Vmatrix}_2 \\ \Leftrightarrow && \min_{x}~\begin{Vmatrix} \begin{pmatrix}1 + t_1 \\ \vdots \\ 1 + t_n \end{pmatrix} \begin{pmatrix}x_0 \\x_1 \end{pmatrix} - \begin{pmatrix}y_1 \\ \vdots \\ y_n \end{pmatrix} \end{Vmatrix}_2 \end{eqnarray*}$

Vergleiche den Wikipedia Link für die in diesem Fall direkte Berechnung des Lösungsvekotrs x.

26.07.2007, 02:28

tigerbine

Auf diesen Beitrag antworten »

8b. Allgemeiner Linearer Fall
Hängt die Modellfunktion von mehreren Variablen ab, so ergibt sich :

$\begin{eqnarray*} y_m(t) = x_0 + x_1\cdot t_1 + ...+x_q\cdot t_q \end{eqnarray*}$

Das hieraus resultierende LGS führt zu der Darstellung:

$\begin{eqnarray*} \min_x~\begin{Vmatrix}\begin{pmatrix}1&t_{11} & ...&t_{1q} \\\vdots &&&\vdots \\ 1 &t_{n1} & ... &t_{nq} \end{pmatrix}\begin{pmatrix} x_0 \\ \vdots \\x_q \end{pmatrix}-\begin{pmatrix} y_1 \\ \vdots \\ y_n\end{pmatrix} \end{Vmatrix}_2 \end{eqnarray*}$

Dies ist nun ein Problem, dessen Lösbarkeit wir bereits hier am Anfang besprochen haben. Lösungsverfahren sind das Normalengleichungsverfahren (Beachte! die Kondition der Matrizen) oder die QR-Zerlegung.

16.08.2007, 14:07

tigerbine

Auf diesen Beitrag antworten »

Ausblick
[Workshop - Lineare Ausgleichsprobleme - Beispiele]

Neue Frage »

Antworten »

[WS] Lineare Ausgleichprobleme

Verwandte Themen