diff --git a/semester3/numcs/numcs-summary.pdf b/semester3/numcs/numcs-summary.pdf index 638a3a1..72d9636 100644 Binary files a/semester3/numcs/numcs-summary.pdf and b/semester3/numcs/numcs-summary.pdf differ diff --git a/semester3/numcs/parts/04_linalg/00_intro.tex b/semester3/numcs/parts/04_linalg/00_intro.tex index affe42c..a450613 100644 --- a/semester3/numcs/parts/04_linalg/00_intro.tex +++ b/semester3/numcs/parts/04_linalg/00_intro.tex @@ -41,27 +41,52 @@ Wenn sie zudem normiert sind (also $||q_i||_2 = 1 \smallhspace \forall i \leq n$ \end{align*} -\textbf{Perturbierte LGS} +\drmvspace +\shade{gray}{Perturbierte LGS} Statt $Ax = b$ ist das LGS ungenau gegeben: $(A + \Delta A)(\tilde{x} - x) = \Delta b - \Delta Ax$. -$\text{cond}(A) := \left\lvert\left\lvert A^{-1} \right\rvert\right\rvert \cdot \lvert\lvert A \rvert\rvert \in \mathbb{R}$ (Konditionszahl von $A$) +\fancydef{Konditionszahl} $\text{cond}(A) := ||A^{-1}|| \cdot ||A||$. Manchmal auch mit $\kappa(A)$ notiert + +Auch hier gibt es sie wieder für verschiedene Normen: +\begin{itemize}[noitemsep] + \item $\kappa_2(A) = \frac{\sigma_\text{max} (A)}{\sigma_\text{min}(A)}$ (Spektralnorm mit Singulärwerten) + \item $\kappa_\infty(A) = ||A||_\infty \cdot ||A^{-1}||_\infty$ + \item $\kappa_1(A) = ||A||_1 \cdot ||A^{-1}||_1$ +\end{itemize} $\text{cond}(A) \gg 1$ bedeutet intuitiv: kleine Änderung der Daten $\mapsto$ grosse Änderung in der Lösung -\textbf{Grosse Matrizen} +Zudem haben wir folgende Eigenschaften: +\drmvspace +\begin{multicols}{2} + \begin{itemize}[noitemsep] + \item $\kappa(A) \geq 1$ + \item $\kappa(cA) = \kappa(A) \smallhspace \forall c \neq 0$ + \item $\kappa(A) = \kappa(A^{-1})$ + \item Für orthogonale und unitäre Matrizen $Q$: $\kappa_2(Q) = 1$ + \end{itemize} +\end{multicols} + + + +\drmvspace +\shade{gray}{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. + +% TODO: Update with notes from TA and script \textbf{Dünnbesetzte Matrizen} + $\text{nnz}(A) := |\{ (i,j) \ |\ a_{ij} \in A, a_{ij} \neq 0 \}| \ll m\cdot n$ $\limit{l}{\infty} \frac{\text{nnz}(A^{(l)})}{n_l m_l} = 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) +Es gibt viele Formate, je nach Anwendung sind gewisse sinnvoller als andere. (Siehe Tabelle, NumCSE) % TODO: Insert here \verb|scipy.sparse.csr_matrix(A)| $\mapsto$ Dramatische Speichereinsparung.\\ Deprecated: \verb|bsr_array| und \verb|coo_array| verwenden, kompatibel mit \verb|numpy| arrays. diff --git a/semester3/numcs/parts/04_linalg/01_LU.tex b/semester3/numcs/parts/04_linalg/01_LU.tex index 94e84fa..73d4c8a 100644 --- a/semester3/numcs/parts/04_linalg/01_LU.tex +++ b/semester3/numcs/parts/04_linalg/01_LU.tex @@ -18,6 +18,25 @@ Jedoch ist dies nicht sinnvoll, wenn wir die Dreiecksmatrizen gar nicht benötig $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. +Diese Zerlegung kann $Ax = b$ potenziell schneller lösen als LU und wir verwenden nur halb so viel Speicher. +Zudem ist keine Pivotierung nötig, also ist das Verfahren für symmetrisch positiv definite Matrizen numerisch stabil. + +Im Folgenden ist der Cholesky algorithmus in Pseudocode beschrieben: +\begin{algorithm} + \begin{spacing}{1.2} + \caption{\textsc{cholesky}(A)} + \begin{algorithmic}[1] + \State $n \gets \texttt{A.shape[0]}$ + \State $l \gets$ Initialisiere ein $n \times n$ array + \For{$j = 1, 2, \ldots, n$} + \State $l_{jj} \gets \sqrt{A_{jj} - \sum_{k = 1}^{j - 1} l_{jk}^2}$ + \For{$i = j + 1, \ldots, n$} + \State $l_{ij} \gets \displaystyle\frac{1}{l_{jj}} \left( A_{ij} - \sum_{k = 1}^{j - 1} l_{ik} l_{jk} \right)$ + \EndFor + \EndFor + \State \Return $l$ + \end{algorithmic} + \end{spacing} +\end{algorithm} \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 diff --git a/semester3/numcs/parts/04_linalg/02_qr.tex b/semester3/numcs/parts/04_linalg/02_qr.tex index e69de29..148ad83 100644 --- a/semester3/numcs/parts/04_linalg/02_qr.tex +++ b/semester3/numcs/parts/04_linalg/02_qr.tex @@ -0,0 +1,134 @@ +\subsubsection{QR-Zerlegung} +Wir können eine Matrix $A \in \R^{m \times n}$ mit $m \geq n$ als $A = QR$ zerlegen, wobei $Q \in \R^{m \times n}$ +orthonormale Spalten besitzt und $R \in \R^{n \times n}$ ist eine obere Dreiecksmatrix. + +Die QR-Zerlegung ist numerisch stabiler bei schlecht konditionierten Problemen, ist die Basis vieler Eigenwertverfahren und ist ideal für Least Squares. + +Wir können mit der QR-Zerlegung auch LGS $Ax = b$ lösen: +\rmvspace +\begin{align*} + Ax = b \Longleftrightarrow QRx = b \Longleftrightarrow Rx = Q^\top b +\end{align*} + +Da $Q$ orthogonal ist, haben wir $Q^{-1} = Q^\top$. +Um das Ganze einfacher zu machen, lösen wir das System $Rx = y$, wobei $y := Q^\top b$. +In numpy können wir direkt mit \texttt{np.linalg.solve()} dies lösen (nutzt automatisch Rückwärtssubstitution) + + +\newpage +\bg{purple}{Givens-Rotations} + +Bei der Givens-Rotation generiert man eine Rotationsmatrix, die die $(i, j)$-Ebene um einen Winkel $\theta$ rotiert. +Die dazu konstruierte Matrix hat dabei die folgende Form (rechts eine Beschreibung des Eintrags $(k, l)$): +\rmvspace +\begin{align*} + G(i, j, \theta) = + \begin{bmatrix} + 1 & 0 & & \cdots & & & 0 \\ + 0 & \ddots & & & & & \\ + & & c & \cdots & s & & \\ + \vdots & & \vdots & & \vdots & & \vdots \\ + & & -s & \cdots & c & & \\ + 0 & & & & & \ddots & \\ + 0 & & & \cdots & & & 1 \\ + \end{bmatrix} + \text{ oder } + G(i, j, \theta)_{k, l} = + \begin{cases} + k = l \land k \neq i, j & 1 \\ + k = l \land (k = i \lor k = j) & c \\ + k = i \land l = j & -s \\ + k = j \land l = i & s \\ + \text{sonst} & 0 + \end{cases} +\end{align*} + +\drmvspace +Dabei ist $c = \cos(\theta)$ und $s = \sin(\theta)$. +Diese Matrix hat einige nützliche Eigenschaften: $G^\top G = I$ (also ist $G$ orthogonal), also gilt auch $G^{-1} = G^\top$ und $G$ modifiziert nur Zeilen $i$ und $j$ + +Im Zweidimensionalen Raum können wir die Werte für $c$ und $s$ so bestimmen: +\rmvspace +\begin{align*} + \begin{bmatrix} + c & s \\ + -s & c + \end{bmatrix} + \begin{bmatrix} + a \\ b + \end{bmatrix} + = + \begin{bmatrix} + r \\ 0 + \end{bmatrix} \\ + r = \sqrt{a^2 + b^2}, \smallhspace c = \frac{a}{r}, \smallhspace s = \frac{b}{r} +\end{align*} +Wir haben jedoch das Problem, dass die Berechnung von $r$ überlaufen kann. Dies lösen wir, indem wir skalieren: +\drmvspace +\begin{multicols}{2} + Falls $|b| > |a|$: + \rmvspace + \begin{align*} + t = \frac{a}{b}, \smallhspace s = \frac{1}{\sqrt{1 + t^2}}, \smallhspace c = s \cdot t + \end{align*} + + Falls $|a| \geq |b|$: + \drmvspace + \begin{align*} + t = \frac{b}{a}, \smallhspace s = c \cdot t, \smallhspace c = \frac{1}{\sqrt{1 + t^2}} + \end{align*} +\end{multicols} +Es ist wichtig, dass wir das $r = \text{sign}(a) \sqrt{a^2 + b^2}$ mit Vorzeichen berechnen, um Auslöschung zu vermeiden + + +\newpage +\bg{purple}{Gram-Schmidt} + +Die Idee des Gram-Schmidt-Algorithmus ist es, orthonormale Vektoren zu konstruieren und diese dann zur Matrix $Q$ zusammenzubasteln. + +Es wurden zwei Algorithmen behandelt, beide unten in Pseudocode: +\begin{algorithm} + \begin{spacing}{1.2} + \caption{\textsc{classicalGramSchmidt}(A)} + \begin{algorithmic}[1] + \State $n \gets \texttt{A.shape[0]}$ + \State $q \gets$ Initialisiere ein $n \times n$ array + \State $r \gets$ Initialisiere ein $n \times n$ array + \For{$k = 1, 2, \ldots, n$} + \State $v_k \gets a_k$ \Comment{Der $k$-te Spaltenvektor} + \For{$i = 1, \ldots, k - 1$} + \State $r_{ik} \gets q_i^\top a_k$ + \State $v_k \gets a_k - \sum_{i = 1}^{k - 1} r_{ik} q_i$ + \EndFor + \State $r_{kk} = ||v_k||$ + \State $q_k = v_k / r_{kk}$ \Comment{Vektor normieren} + \EndFor + \State \Return $q$, $r$ + \end{algorithmic} + \end{spacing} +\end{algorithm} + +\drmvspace +\begin{algorithm} + \begin{spacing}{1.2} + \caption{\textsc{modifiedGramSchmidt}(A)} + \begin{algorithmic}[1] + \State $n \gets \texttt{A.shape[0]}$ + \State $q \gets$ Initialisiere ein $n \times n$ array + \State $r \gets$ Initialisiere ein $n \times n$ array + \For{$k = 1, 2, \ldots, n$} + \State $v_k \gets a_k$ \Comment{Der $k$-te Spaltenvektor} + \For{$i = 1, \ldots, k - 1$} + \State $r_{ik} \gets q_i^\top v_k$ + \State $v_k \gets v_k - r_{ik} q_i$ + \EndFor + \State $q_k = v_k / r_{kk}$ \Comment{Vektor normieren} + \EndFor + \State \Return $q$, $r$ + \end{algorithmic} + \end{spacing} +\end{algorithm} +Falls $R$ nicht benötigt wird, kann viel Speicher gespart werden, indem man das $r_{ik}$ als eine scoped variable verwendet. + + +\bg{purple}{Householder-Reflektor} diff --git a/semester3/numcs/parts/04_linalg/03_svd.tex b/semester3/numcs/parts/04_linalg/03_svd.tex index e69de29..d1254b9 100644 --- a/semester3/numcs/parts/04_linalg/03_svd.tex +++ b/semester3/numcs/parts/04_linalg/03_svd.tex @@ -0,0 +1 @@ +\subsubsection{Singulärwertzerlegung}