[NumCS] Finish zeros section, start update to linalg section

semester3/numcs/format.py

@@ -0,0 +1,33 @@
+import numpy as np
+import scipy as sp
+
+
+def fastbroyd(x0, F, J, tol=1e-12, maxit=20):
+    x = x0.copy()  # make sure we do not change the iput
+    lup = sp.linalg.lu_factor(J)  # LU decomposition of J
+    s = sp.linalg.lu_solve(lup, F(x))  # start with a Newton corection
+    sn = np.dot(s, s)  # squared norm of the correction
+    x -= s
+    f = F(x)  # start with a full Newton step
+    dx = np.zeros((maxit, len(x)))  # containers for storing corrections s and their sn:
+    dxn = np.zeros(maxit)
+    k = 0
+    dx[k] = s
+    dxn[k] = sn
+    k += 1  # the number of the Broyden iteration
+
+    # Broyden iteration
+    while sn > tol and k < maxit:
+        w = sp.linalg.lu_solve(lup, f)  # f = F (actual Broyden iteration x)
+        # Using the Sherman-Morrison-Woodbury formel
+        for r in range(1, k):
+            w += dx[r] * (np.dot(dx[r - 1], w)) / dxn[r - 1]
+            z = np.dot(s, w)
+            s = (1 + z / (sn - z)) * w
+            sn = np.dot(s, s)
+            dx[k] = s
+            dxn[k] = sn
+            x -= s
+            f = F(x)
+            k += 1  # update x and iteration number k
+    return x, k  # return the final value and the numbers of iterations needed

semester3/numcs/numcs-summary.tex

@@ -161,6 +161,9 @@ Moral of the story: Use descriptive variable names and do NOT use $t$, $tt$, $tt
 \newsection
 \section{Intermezzo: Lineare Algebra}
 \input{parts/04_linalg/00_intro.tex}
+\input{parts/04_linalg/01_LU.tex}
+\input{parts/04_linalg/02_qr.tex}
+\input{parts/04_linalg/03_svd.tex}
 
 \newsection
 \section{Ausgleichsrechnung}

semester3/numcs/parts/03_zeros/08_quasi-newton.tex

@@ -10,4 +10,43 @@ Im höherdimensionalen Raum ist dies jedoch nicht direkt möglich und wir erhalt
 
 \drmvspace
 Dabei ist $J_0$ z.B. durch $DF(x^{(0)})$ definiert.
-% Page 222
+
+\fancyremark{Broyden-Update} Das Broyden-Update ergibt bezüglich der $||\cdot||_2$-Norm die minimale korrektur der Jakobi-Matrix $J_k$ an, so dass die Sekantenbedingung erfüllt ist.
+Die Implementierung erzielt man folgendermassen mit der \bi{Sherman-Morrison-Woodbury} Update-Formel:
+\begin{align*}
+    J_{k + 1}^{-1} = \left( I - \frac{J_k^{-1} F(x^{(k + 1)}) (\Delta x^{(k)})^\top}{||\Delta x^{(k)}||_2^2 + (\Delta x^{(k)})^\top J_k^{-1} F(x^{(k + 1)})} \right) J^{-1}_k
+\end{align*}
+
+Das Broyden-Quasi-Newton-Verfahren konvergiert langsamer als das Newton-Verfahren, aber schneller als das vereinfachte Newton-Verfahren. (\texttt{sp} ist \texttt{Scipy} und \texttt{np} logischerweise \texttt{Numpy} im untenstehenden code)
+
+\begin{code}{python}
+def fastbroyd(x0, F, J, tol=1e-12, maxit=20):
+    x = x0.copy()  # make sure we do not change the iput
+    lup = sp.linalg.lu_factor(J)  # LU decomposition of J
+    s = sp.linalg.lu_solve(lup, F(x))  # start with a Newton corection
+    sn = np.dot(s, s)  # squared norm of the correction
+    x -= s
+    f = F(x)  # start with a full Newton step
+    dx = np.zeros((maxit, len(x)))  # containers for storing corrections s and their sn:
+    dxn = np.zeros(maxit)
+    k = 0
+    dx[k] = s
+    dxn[k] = sn
+    k += 1  # the number of the Broyden iteration
+
+    # Broyden iteration
+    while sn > tol and k < maxit:
+        w = sp.linalg.lu_solve(lup, f)  # f = F (actual Broyden iteration x)
+        # Using the Sherman-Morrison-Woodbury formel
+        for r in range(1, k):
+            w += dx[r] * (np.dot(dx[r - 1], w)) / dxn[r - 1]
+            z = np.dot(s, w)
+            s = (1 + z / (sn - z)) * w
+            sn = np.dot(s, s)
+            dx[k] = s
+            dxn[k] = sn
+            x -= s
+            f = F(x)
+            k += 1  # update x and iteration number k
+    return x, k  # return the final value and the numbers of iterations needed
+\end{code}

semester3/numcs/parts/04_linalg/00_intro.tex

@@ -1,5 +1,45 @@
-% FIXME: Not true anymore!
+\subsection{Grundlagen}
-Das Skript-Kapitel hierzu existiert nicht. Die Vorlesung orientiert sich an dem äquivalenten Teil des NumCSE Dokuments.
+\inlineremark Eine Tabelle mit invertierbaren und nicht invertierbaren Matrizen findet sich unten:
+\begin{tables}{ll}{Invertierbar                                             & Nicht Invertierbar}
+              $A$ ist regulär                                               & $A$ ist singulär                                                  \\
+              Spalten sind linear unabhängig                                & Spalten sind linear abhängig                                      \\
+              Zeilen sind linear unabhängig                                 & Zeilen sind linear abhängig                                       \\
+              $\det(A) \neq 0$                                              & $\det(A) = 0$                                                     \\
+              $Ax = 0$ hat eine Lösung $x = b$                              & $Ax = 0$ hat unendlich viele Lösungen                             \\
+              $Ax = b$ hat eine Lösung $x = A^{-1}b$                        & $Ax = b$ hat keine oder unendlich viele Lösungen                  \\
+              $A$ hat vollen Rang                                           & $A$ hat Rang $r < n$                                              \\
+              $A$ hat $n$ non-zero Pivots                                   & $A$ hat $r < n$ Pivots                                            \\
+              $\text{span}\{A_{:, 1}, \ldots, A_{:, n}\}$ hat Dimension $n$ & $\text{span}\{A_{:, 1}, \ldots, A_{:, n}\}$ hat Dimension $r < n$ \\
+              $\text{span}\{A_{1, :}, \ldots, A_{n, :}\}$ hat Dimension $n$ & $\text{span}\{A_{1, :}, \ldots, A_{n, :}\}$ hat Dimension $r < n$ \\
+              Alle Eigenwerte von $A$ sind nicht Null                       & $0$ ist der Eigenwert von $A$                                     \\
+              $0 \notin \sigma(A) =$ Spektrum von $A$                       & $0 \in \sigma(A)$                                                 \\
+              $A^H A$ ist symmetrisch positiv definit                       & $A^H A$ ist nur semidefinit                                       \\
+              $A$ hat $n$ (positive) Singulärwerte                          & $A$ hat $r < n$ (positive) Singulärwerte                          \\
+\end{tables}
+
+\fancydef{Orthogonale Vektoren} Vektoren $q_1, \ldots, q_n$ heissen \bi{orthogonal}, falls
+\rmvspace
+\begin{align*}
+    q_i^H \cdot q_j = 0 \smallhspace \forall i, j \leq n \text{ with } i \neq j
+\end{align*}
+
+\drmvspace
+Wenn sie zudem normiert sind (also $||q_i||_2 = 1 \smallhspace \forall i \leq n$), dann heissen sie \bi{orthonormal}
+
+\inlineremark In der vorigen Definition wird die \bi{Euklidische Norm} $||q||_2^2 = q^H \cdot q$ verwendet
+
+\setLabelNumber{all}{7}
+\fancyremark{Rotationen} Die Rotationsmatrix für eine Rotation um Winkel $\theta$ ist gegeben durch:
+\rmvspace
+\begin{align*}
+    R_\theta =
+    \begin{bmatrix}
+        \cos(\theta) & 0 & -\sin(\theta) \\
+        0            & 1 & 0             \\
+        \sin(\theta) & 0 & \cos(\theta)
+    \end{bmatrix}
+\end{align*}
+
 
 \textbf{Perturbierte LGS}
 
@@ -9,22 +49,6 @@ $\text{cond}(A) := \left\lvert\left\lvert A^{-1} \right\rvert\right\rvert \cdot
 
 $\text{cond}(A) \gg 1$ bedeutet intuitiv: kleine Änderung der Daten $\mapsto$ grosse Änderung in der Lösung
 
-\textbf{Gauss Elimination / LU Zerlegung}
-
-$A \in \mathbb{R}^{n\times m} = PLU$ wie bekannt.
-
-\textbf{Cholesky Zerlegung} ($A$ pos. def. und hermitesch)
-
-$A = LDL^\top = \underbrace{L\sqrt{D}}_{R^\top}\underbrace{\sqrt{D}L^\top}_{R} = R^\top R$
-
-Kann $Ax = b$ potenziell schneller lösen als LU.
-
-\begin{code}{python}
-    L = np.linalg.solve(A)      # A = L @ L.T
-    y = np.linalg.solve(L, b)
-    x = np.linalg.solve(L.T, y)
-\end{code}
-
 \textbf{Grosse Matrizen}
 
 Passen oft nicht (direkt) in den Speicher: effizientere Speicherung nötig, möglich für z.B. Diagonalmatrizen, Dreiecksmatrizen. Auch für Cholesky möglich.
@@ -33,9 +57,9 @@ Passen oft nicht (direkt) in den Speicher: effizientere Speicherung nötig, mög
 
 $\text{nnz}(A) := |\{ (i,j) \ |\ a_{ij} \in A, a_{ij} \neq 0 \}| \ll m\cdot n$
 
-$\underset{l \to \infty}{\lim} \frac{\text{nnz}(A^{(l)})}{n_l m_l} = 0$
+$\limit{l}{\infty} \frac{\text{nnz}(A^{(l)})}{n_l m_l} = 0$
 
-Einfacher zu speichern: \verb|val, col, row| vektoren s.d. \verb|val[k]| $ = a_{ij}$ wobei $i=$ \verb|row[k]|, $j=$ \verb|col[k]|. (nur $a_{ij} \neq 0$)
+Einfacher zu speichern: \verb|val, col, row| sind Vektoren so dass \verb|val[k]| $ = a_{ij}$, wobei $i=$ \verb|row[k]|, $j=$ \verb|col[k]|. (nur $a_{ij} \neq 0$)
 
 Viele Formate, je nach Anwendung gewisse sinnvoller als andere. (Siehe Tabelle, NumCSE)
 

semester3/numcs/parts/04_linalg/01_LU.tex

@@ -0,0 +1,23 @@
+\subsubsection{Gauss Elimination / LU Zerlegung}
+
+Das Anwenden der Gauss-Elimintation ergibt die LU-Zerlegung, gegeben durch $A \in \mathbb{R}^{n\times m} = PLU$,
+wobei $U$ eine obere Dreiecksmatrix (die resultierende Matrix der Gauss-Elimintation), $L$ eine untere Dreiecksmatrix (Matrix aller Schritte der Gauss-Elimintation)
+und $P$ eine Permutationsmatrix ist.
+
+\innumpy können wir \texttt{P, L, U = scipy.linalg.lu(A)} (\texttt{Numpy} liefert keine LU-Zerlegung). Mit \texttt{scipy.linalg.lu\_solve(P, L, U)} kann man dann das System lösen.
+Jedoch ist dies nicht sinnvoll, wenn wir die Dreiecksmatrizen gar nicht benötigen. In diesem Fall verwenden wir einfach \texttt{numpy.linalg.solve(A)}
+
+\begin{code}{python}
+    L = np.linalg.solve(A)      # A = L @ L.T
+    y = np.linalg.solve(L, b)
+    x = np.linalg.solve(L.T, y)
+\end{code}
+
+
+\shade{Cyan}{Cholesky Zerlegung} ($A$ ist positiv defefinit und hermetisch)
+
+$A = LDL^H = \underbrace{L\sqrt{D}}_{R^H}\underbrace{\sqrt{D}L^H}_{R} = R^H R$
+
+Diese Zerlegung kann $Ax = b$ potenziell schneller lösen als LU.
+
+\innumpy haben wir via \texttt{scipy.linalg} die Funktionen \texttt{cholesky}, \texttt{cho\_factor} und \texttt{cho\_solve}, wie auch bereits äquivalent für die LU-Zerlegung