diff --git a/semester3/numcs/format.py b/semester3/numcs/format.py index e69de29..6512e81 100644 --- a/semester3/numcs/format.py +++ b/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 diff --git a/semester3/numcs/numcs-summary.pdf b/semester3/numcs/numcs-summary.pdf index 642b060..f6f8e20 100644 Binary files a/semester3/numcs/numcs-summary.pdf and b/semester3/numcs/numcs-summary.pdf differ diff --git a/semester3/numcs/numcs-summary.tex b/semester3/numcs/numcs-summary.tex index d422199..2b8ff24 100644 --- a/semester3/numcs/numcs-summary.tex +++ b/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} diff --git a/semester3/numcs/parts/03_zeros/08_quasi-newton.tex b/semester3/numcs/parts/03_zeros/08_quasi-newton.tex index 855743d..dba29ff 100644 --- a/semester3/numcs/parts/03_zeros/08_quasi-newton.tex +++ b/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} diff --git a/semester3/numcs/parts/04_linalg/00_intro.tex b/semester3/numcs/parts/04_linalg/00_intro.tex index 9347fcf..affe42c 100644 --- a/semester3/numcs/parts/04_linalg/00_intro.tex +++ b/semester3/numcs/parts/04_linalg/00_intro.tex @@ -1,5 +1,45 @@ -% FIXME: Not true anymore! -Das Skript-Kapitel hierzu existiert nicht. Die Vorlesung orientiert sich an dem äquivalenten Teil des NumCSE Dokuments. +\subsection{Grundlagen} +\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) diff --git a/semester3/numcs/parts/04_linalg/01_LU.tex b/semester3/numcs/parts/04_linalg/01_LU.tex index e69de29..94e84fa 100644 --- a/semester3/numcs/parts/04_linalg/01_LU.tex +++ b/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