Add summary about $QR$-decomposition

Notes/Semester 3/HM1 - Höhere Mathematik/Zusammenfassung.md

@@ -4,7 +4,12 @@
     margin-bottom: 2px;
     color: black;
     padding: 0.5rem;
-    border: solid black 2px;
+    padding-bottom: 1px;
+    margin-bottom: 0.5rem;
+  }
+
+  .formula p:last-child, .letters p:last-child {
+    padding-bottom: 0;
   }
 
   .formula {
@@ -42,6 +47,15 @@
     - [Formelbuchstaben zu Nullstellenproblem](#formelbuchstaben-zu-nullstellenproblem)
   - [Lineare Gleichungssysteme](#lineare-gleichungssysteme)
     - [Lernziele](#lernziele)
+    - [Eigenschaften](#eigenschaften)
+    - [Dreiecks-Matrizen](#dreiecks-matrizen)
+    - [Der Gauss-Algorithmus](#der-gauss-algorithmus)
+    - [Fehlerfortpflanzung und Pivotisierung](#fehlerfortpflanzung-und-pivotisierung)
+    - [Determinanten-Bestimmung](#determinanten-bestimmung)
+    - [Die $LR$-Zerlegung](#die-lr-zerlegung)
+    - [$QR$-Zerlegung](#qr-zerlegung)
+      - [Housholder-Matrizen](#housholder-matrizen)
+      - [Vorgang](#vorgang)
   - [Formelbuchstaben](#formelbuchstaben)
   - [Glossar](#glossar)
 
@@ -54,7 +68,7 @@
 
 ### Arten von Lösungen
   - Direkte Verfahren - Exakte Lösung nach endlicher Zeit
-  - Näherungsverfahren - Approximation nach begrenzter Anzahl Rechenschritte
+  - Näherungsverfahren/Iteratives Verfahren - Approximation nach begrenzter Anzahl Rechenschritte
 
 ### Verbindung zur Informatik
   - Effiziente Berechnung numerischer Algorithmen
@@ -188,6 +202,9 @@ Es wird der korrekte Wert $x$ für eine Aufgabe gesucht.
 >     - Ein Threshold von $0$ ergibt das genaue Resultat
 
 ### Fixpunktiteration
+
+![](FixedPointIteration.png) 
+
 Ein möglicher Ansatz für ein solches Problem ist eine Fixpunktiteration.
 
 Der Vorgang für eine solche ist folgende:
@@ -305,6 +322,9 @@ Folgendermassen kann dieser aufgestellt werden:
  4. Die a-priori und die a-posteriori Abschätzung kann nun beliebig angewendet werden. Hierbei wird für $x_0$ der Wert $a$ verwendet.
 
 ### Newton-Verfahren
+
+![](NewtonMethod.png) 
+
 Das Newton-Verfahren erreicht die Konvergenz (d.h. das (approximierte) Resultat) um einiges schneller.
 
 Hierfür wird die Funktion $f$ in der Nullstellenform benötigt ($f(x) = \text{[...]} = 0$).
@@ -338,6 +358,8 @@ Das Ergebnis ist wahr, wenn mit dem gewählten $x$ eine Konvergenz erreicht werd
 
 ### Sekantenverfahren
 
+![](SecantMethod.png)
+
 <div class="formula">
 
 $$x_{n + 1} = x_n - \frac{x_n - x_{n - 1}}{f(x_n) - f(x_{n - 1})} \cdot f(x_n)$$
@@ -425,17 +447,462 @@ Vorgang:
     - [ ] Fehlerabschätzungen
     - [ ] Eigenwerte und Eigenvektoren von Matrizen berechnen
 
+<div class="letters">
+
+**Lineares Gleichungssystem:**
+
+Lineare Gleichungssysteme haben jeweils die Form $A \cdot x = b$ wobei $A$ und $b$ gegeben und $x$ gesucht ist:
+
+$$A = \left(
+  \begin{matrix}
+    a_{11} & a_{12} & \cdots & a_{1n} \\
+    a_{21} & a_{22} & \cdots & a_{2n} \\
+    \vdots & \vdots & & \vdots \\
+    a_{n1} & a_{n2} & \cdots & a_{nn}
+  \end{matrix}
+\right),
+x = \left(
+  \begin{matrix}
+    x_1 \\
+    x_2 \\
+    \vdots \\
+    x_n
+  \end{matrix}
+\right),
+b = \left(
+  \begin{matrix}
+    b_1 \\
+    b_2 \\
+    \vdots \\
+    b_n
+  \end{matrix}
+\right)$$
+
+</div>
+
+### Eigenschaften
+  - Gleich viele gesuchte Variablen $x_n$ wie Gleichungen $n$. Folglich:
+    - Die Matrix $A$ ist eine quadratische Matrix mit Dimensionen $n \times n$
+  - $A$ ist invertierbar
+  - $A$ hat eine Determinante $\det(A)$
+
+### Dreiecks-Matrizen
+
+<div class="letters">
+
+***$L$: Untere Dreiecksmatrix***
+
+Eine Matrix, die in der oberen rechten Ecke nur den Wert $0$ und auf der Diagonale nur den Wert $1$ hat. Eine Untere Dreiecksmatrix hat also folgende Form:
+
+$$L = \left(
+  \begin{matrix}
+    1 & 0 & 0 & \cdots & 0 \\
+    l_{21} & 1 & 0 & \cdots & 0 \\
+    l_{31} & l_{32} & 1 & \ddots & 0 \\
+    \vdots & \vdots & \ddots & \ddots & 0 \\
+    l_{n1} & l_{n2} & \cdots & l_{nn - 1} & 1
+  \end{matrix}
+\right)$$
+
+***$R$: Obere Dreiecksmatrix***
+
+Eine Matrix, die unten links von der Diagonale nur den Wert $0$ beinhaltet. Eine Obere Dreiecksmatrix hat dementsprechend folgende Form:
+
+$$R = \left(
+  \begin{matrix}
+    r_{11} & r_{12} & r_{13} & \cdots & r_{1n} \\
+    0 & r_{22} & r_{23} & \cdots & r_{2n} \\
+    0 & 0 & r_{33} & \cdots & r_{3n} \\
+    \vdots & \vdots & \ddots & \ddots & \vdots \\
+    0 & 0 & \cdots & 0 & r_{nn}
+  \end{matrix}
+\right)$$
+
+</div>
+
+***Code-Beispiele:***
+
+_Umwandlung in $R$-Matrix:_
+
+```py
+for i in range(n):
+    if A[i, i] == 0:
+        index = -1
+        for j in range(i + 1, n):
+            if A[j, i] > 0:
+                index = j
+        if index == -1:
+            raise Exception("Invalid Matrix")
+        else:
+            # Swap lines
+            A[[i, index]] = A[[index, i]]
+    for j in range(i + 1, n):
+        factor = A[j, i] / A[i, i]
+        A[j] = A[j] - (factor * A[i])
+```
+
+### Der Gauss-Algorithmus
+Der Gauss-Algorithmus basiert darauf, dass ein lineares Gleichungssystem leicht lösbar ist, falls $A$ eine obere Dreiecksmatrix ist. $A$ muss also hierfür die Form einer oberen Dreiecksmatrix $R$ haben.
+
+<div class="formula">
+
+***Gauss-Algorithmus:***
+
+$$x_i = \frac{b_i - \sum_{j = i + 1}^n{a_{ij} \cdot x_j}}{a_{ii}}, i = n, n - 1, \dots, 1$$
+
+</div>
+
+Um den Gauss-Algorithmus anzuwenden, muss die Matrix $A$ erst in eine $R$-Matrix umgewandelt werden. Dies funktioniert wie folgt:
+
+ 1. Mit $i$ von $1$ bis $n$
+ 2. Falls $a_{ii}$ den Wert $0$ hat:
+    1. Mit $j$ von $i + 1$ bis $n$
+    2. Prüfen, ob $a_{ji}$ einen höheren Wert als $0$ hat
+       - Falls Zeile gefunden wurde:
+         - $a_{i}$ mit $a_{j}$ tauschen
+         - $b_{i}$ mit $b_{j}$ tauschen
+       - Sonst beenden: ungültige Matrix
+ 3. Mit $j$ von $i + 1$ bis $n$
+    1. $a_k = a_k - \frac{a_{ki}}{a_{ii}} \cdot a_i$
+    2. $b_k = b_k - \frac{a_{ki}}{a_{ii}} \cdot b_i$
+
+***Code-Beispiel:***
+
+```py
+from numpy import array, zeros
+
+def gaussMethod(A, b):
+    A = array(A)
+    n = A.shape[0]
+    A = A.reshape((n, n))
+    b = array(b).reshape((n))
+    result = zeros(n)
+
+    # Convert to R-Matrix
+    for i in range(n):
+        maxIndex = i
+        for j in range(i + 1, n):
+            if A[j, i] > A[maxIndex, i]:
+                maxIndex = j
+        # Swap lines
+        A[[i, maxIndex]] = A[[maxIndex, i]]
+        b[[i, maxIndex]] = b[[maxIndex, i]]
+        for j in range(i + 1, n):
+            factor = A[j, i] / A[i, i]
+            A[j] = A[j] - (factor * A[i])
+            b[j] = b[j] - (factor * b[i])
+    # Calculate result
+    for index in range(n, 0, -1):
+        i = index - 1
+        value = b[i]
+        for j in range(i, n):
+            value = value - A[i, j] * result[j]
+        result[i] = value / A[i, i]
+    return result.reshape((n, 1))
+```
+
+### Fehlerfortpflanzung und Pivotisierung
+  - Da beim Umwandeln einer Matrix $A$ in die $R$-Form Zeilen in jedem Schritt mit dem Faktor $\lambda = \frac{a_{ji}}{a_{ii}}$ multipliziert werden, vergrössert sich der Schritt immer um $|\lambda|$
+  - $\lambda$ kann klein gehalten werden, indem Zeilen der Grösse nach sortiert werden
+  - In den Code-Beispielen ist dies bereits berücksichtigt
+
+### Determinanten-Bestimmung
+Die Determinante einer Matrix $A$ lässt sich einfach berechnen, sobald sie in die $R$-Form gebracht wurde mit folgender Formel:
+
+<div class="formula">
+
+Determinanten-Bestimmung mit Matrix $\tilde{A}$ (die Matrix $A$ in der $R$-Form):
+
+$$\det(A) =
+    (-1)^l \cdot \det(\tilde{A}) =
+    (-1)^l \cdot \prod_{i = 1}^n{\tilde{a_{ii}}}$$
+
+</div>
+
+***Code-Beispiel:***
+
+```py
+from numpy import array
+
+def det(A):
+    l = 0
+    n = A.shape[0]
+    A = A.reshape((n, n))
+
+    # Convert to R-Matrix
+    for i in range(n):
+        maxIndex = i
+        for j in range(i + 1, n):
+            if A[j, i] > A[maxIndex, i]:
+                maxIndex = j
+        # Swap lines
+        A[[i, maxIndex]] = A[[maxIndex, i]]
+        l = l + 1
+        for j in range(i + 1, n):
+            factor = A[j, i] / A[i, i]
+            A[j] = A[j] - (factor * A[i])
+
+    result = 1
+    for i in range(n):
+        result = result * A[i, i]
+    return (-1 ** l) * result
+```
+
+### Die $LR$-Zerlegung
+In der $LR$-Zerlegung wird die Matrix $A$ in die Matrizen $L$ und $R$ aufgeteilt, sodass $A = L \cdot R$ gilt.
+
+Alternative Namen dieses Vorgangs sind ***$LR$-Faktorisierung*** und $LU$-decomposition.
+
+<div class="formula">
+
+Für in $L$ und $R$ zerlegte Matrizen gilt:
+
+$$A \cdot x = b$$
+
+und
+
+$$A \cdot x = L \cdot R \cdot x = L \cdot y = b$$
+
+Aufwand: Berechnung der $LR$-Zerlegung mit Gauss-Algorithmus benötigt ca. $\frac{2}{3}n^3$ Punktoperationen.
+
+</div>
+
+Falls Zeilenvertauschungen stattfinden, entsteht bei der $LR$-Zerlegung eine zusätzliche Permutations-Matrix $P$.
+
+<div class="formula">
+
+Für $L$ und $R$ zerlegte Matrizen mit Permutation $P$ gilt:
+
+$$P \cdot A = L \cdot R$$
+
+$$L \cdot y = P \cdot b$$
+
+$$R \cdot x = y$$
+
+</div>
+
+Das Verfahren für die $LR$-Zerlegung ist identisch zu den Schritten bei der Umwandlung in eine $R$-Matrix. Jedoch wird jeweils der Wert $l_{ji}$ in der (zu Beginn) leeren Matrix $L$ mit dem im aktuellen Eliminationsschritt gesetzt. Zudem muss bei Vertauschungen die Permutations-Matrix $P$ entsprechend angepasst werden:
+
+***Code-Beispiel:***
+
+```py
+from numpy import array, identity, zeros
+
+def decomposite(A):
+    l = 0
+    n = A.shape[0]
+    R = A.reshape((n, n))
+    L = zeros((n, n))
+    P = identity((n, n))
+    
+    # Convert to LR-Matrix
+    for i in range(n):
+        maxIndex = i
+        for j in range(i + 1, n):
+            if A[j, i] > A[maxIndex, i]:
+                maxIndex = j
+        # Swap lines
+        Pn = identity((n, n))
+        A[[i, maxIndex]] = A[[maxIndex, i]]
+        Pn[[i, maxIndex]] = Pn[[maxIndex, i]]
+        P = P * Pn
+        for j in range(i + 1, n):
+            factor = R[j, i] / R[i, i]
+            L[j, i] = factor
+            R[j] = R[j] - (factor * R[i])
+
+    result = 1
+    for i in range(n):
+        result = result * R[i, i]
+    return [L, R, P]
+```
+
+Wenn die $LR$-Zerlegung, wie in diesem Code, Zeilenaustausch und das Berechnen von $P$ involviert, spricht man von einer $LR$-Zerlegung mit **Spaltenmaximum-Strategie**.
+
+***Vorgang:***
+
+ 1. Gemäss vorhergehender Beschreibung und Code-Beispiel die Matrizen $L$ und $R$ berechnen
+ 2. Mit Hilfe des Gauss-Algorithmus $L \cdot y = P \cdot b$ nach $y$ auflösen
+ 3. Mit Hilfe des Gauss-Algorithmus $R \cdot x = y$ nach $x$ auflösen
+
+### $QR$-Zerlegung
+
+  - Die Matrix $A$ wird in eine orthogonale Matrix $Q$ und eine obere Dreiecksmatrix $R$ zerlegt.
+  - Orthogonal-Matrizen beschreiben Drehungen, Spiegelungen oder Kombinationen daraus.
+  - Eine $QR$-Zerlegung erfordert ca. $\frac{5}{3}n^3$ Punktoperationen - ca. doppelt so viel wie die $LR$-Zerlegung.
+
+<div class="formula">
+
+***Orthogonal-Matrix:***
+
+Eine Matrix $Q$ ist orthogonal, wenn folgendes gilt:
+
+$$Q^T \cdot Q = I_n$$
+
+($x^T$ steht hierbei für eine **T**ransformation)
+
+</div>
+
+#### Housholder-Matrizen
+
+Im Rahmen der Berechnung der Matrizen $Q$ und $R$ werden sogenannte "Housholder-Matrizen" berechnet.
+
+<div class="formula">
+
+***Housholder-Matrizen:***
+
+Sei $u$ ein Vektor mit beliebig vielen Dimensionen, für den gilt:
+
+$$|u| = \sqrt{u_1^2 + u_2^2 + \dots + u_n^2} = 1$$
+
+Die Householder-Matrix hat folgende Eigenschaft:
+
+$$H := I_n - 2 \cdot u \cdot u^T$$
+
+Für Housholder-Matrizen gilt zudem folgendes:
+
+$$H = H^T = H^{-1}$$
+
+und
+
+$$H \cdot H = I_n$$
+
+</div>
+
+***Berechnung einer Housholder-Matrix***
+
+Beispiel der Berechnung einer Housholder-Matrix zur ersten Spalte der Matrix $A$.
+
+> Für die Berechnung wird ein Einheitsvektor $e$ benötigt, welcher genauso viele Werte hat, wie die Matrix Dimensionen. Ein Einheitsvektor hat im ersten Feld den Wert $1$ und in allen anderen Feldern der Wert $0$.
+>
+> Für eine Matrix $A$ mit der Dimension $n = 3$ lautet der Einheitsvektor $e$ also wie folgt:
+> $$e = \left(\begin{matrix}
+>   1 \\
+>   0 \\
+>   0
+> \end{matrix}\right)
+
+ 1. Vektor $v$ bestimmen
+    $$v = a_1 + sign(a_{11}) \cdot |a_1| \cdot e$$
+ 2. Vektor normieren:
+    $$u = \frac{1}{|v|} \cdot v =
+    \frac{1}{\sqrt{1^2 + 2^2 + 3^2}} \cdot
+    \left(\begin{matrix}
+      1 \\
+      2 \\
+      3
+    \end{matrix}\right)  =
+    \frac{1}{\sqrt{14}} \cdot
+    \left(\begin{matrix}
+      1 \\
+      2 \\
+      3
+    \end{matrix}\right)$$
+ 2. Die Housholder-Matrix $H = I_n - 2 \cdot u \cdot u^T$ berechnen.
+    $$H =
+    \left(\begin{matrix}
+      1 & 0 & 0 \\
+      0 & 1 & 0 \\
+      0 & 0 & 1
+    \end{matrix}\right) -
+    2 \cdot \frac{1}{\sqrt{14}} \cdot
+    \left(\begin{matrix}
+      1 \\
+      2 \\
+      3
+    \end{matrix}\right) \cdot
+    \frac{1}{\sqrt{14}} \cdot
+    \left(\begin{matrix}
+      1 & 2 & 3
+    \end{matrix}\right) \\
+    H =
+    \left(\begin{matrix}
+      1 & 0 & 0 \\
+      0 & 1 & 0 \\
+      0 & 0 & 1
+    \end{matrix}\right) -
+    2 \cdot \frac{1}{14} \cdot
+    \left(\begin{matrix}
+      1 & 2 & 3 \\
+      2 & 4 & 6 \\
+      3 & 6 & 9
+    \end{matrix}\right) =
+    -\frac{1}{7} \cdot
+    \left(\begin{matrix}
+      -6 & 2 & 3 \\
+      2 & -3 & 6 \\
+      3 & 6 & 2
+    \end{matrix}\right)$$
+
+<div class="letters">
+
+  - $H$: Housholder-Matrix
+  - $I$: Identitäts-Matrix
+  - $n$: Anzahl Dimensionen der Matrix
+
+</div>
+
+#### Vorgang
+Im Rahmen des Vorgangs entspricht $A_1$ der Matrix $A$.
+
+Die $QR$-Zerlegung kann folgendermassen durchgeführt werden:
+
+  1. $R = A$
+  2. $Q = I_n$
+  3. Für $i$ von $1$ bis $n - 1$
+     1. Gemäss vorheriger Anleitung Householder-Matrix $H_i$ für die erste Spalte von $A_i$ berechnen
+     2. Householder-Matrix um Identitäts-Matrix erweitern. Beispiel:  
+        ![](ExpandHouseholder.png)
+     3. Erweiterte Householder-Matrix als $Q_i$ speichern
+     4. $R = Q_i \cdot R$
+     5. $Q = Q \cdot Q_i^T$
+
+***Code-Beispiel:***
+
+```py
+from numpy import array, identity, sign, sqrt, square, sum, zeros
+
+def qrDecomposition(A):
+    A = array(A)
+    n = A.shape[0]
+    R = A.reshape((n, n))
+    Q = identity(n)
+
+    for i in range(n - 1):
+        I = identity(n - i)
+        Qi = identity(n)
+        e = zeros((n - i, 1))
+        e[0][0] = 1
+        a = R[i:,i:i + 1]
+        v = a + sign(a[0]) * sqrt(sum(square(a))) * e
+        u = (1 / sqrt(sum(square(v)))) * v
+        H = I - 2 * u @ u.T
+        Qi[i:,i:] = H
+        R = Qi @ R
+        Q = Q @ Qi.T
+    return [Q, R]
+```
+
 ## Formelbuchstaben
 <div class="letters">
 
   - $\alpha$: Lipschitz-Konstante (siehe Fixpunktsatz)
   - $[a,b]$: Das Untersuchungs-Interval für den Banachschen Fixpunktsatz
+  - $A$: Matrix eines linearen Gleichungssystems
+  - $\tilde{A}$: Umgewandelte Version der Matrix $A$
+  - $A^T$: Transformierte Matrix $A$
+  - $b$: Das gewünschte Resultat eines linearen Gleichungssystems
   - $B$: Basis der Maschinenzahl
   - $e$: Exponent der Maschinenzahl
+  - $H$: Housholder-Matrix (siehe $QR$-Zerlegung)
+  - $I$: Identitäts-Matrix (Matrix, überall den Wert $0$ und auf der Diagonalen den Wert $1$ hat)
   - $K$: Konditionszahl
+  - $L$: Untere Dreiecksmatrix/Normierte Matrix
   - $m$: Mantisse (Darstellbarer Bereich der Maschinenzahl)
   - $n$: Anzahl möglicher Stellen der Mantisse $m$
   - $q$: Konvergenz-Ordnung
+  - $Q$: Orthogonal-Matrix in der $QR$-Zerlegung
+  - $R$: Obere Dreiecksmatrix
   - $x$: Darzustellender Wert
   - $x_n$: Die $n$-te Approximation von $x$
   - $\tilde{x}$: Approximation/Annäherung an $x$