eth-janishutz/eth-summaries
1
0
Fork 0
mirror of https://github.com/janishutz/eth-summaries.git synced 2026-01-11 19:48:27 +00:00
Code Issues Packages Projects Releases Wiki Activity

[NumCS] QR decomp partially done

Browse Source
This commit is contained in:
Janis Hutz
2025-12-28 13:00:55 +01:00
parent 817fba55a8
commit 259fa8459d
5 changed files with 184 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
semester3/numcs/numcs-summary.pdf
View File

Binary file not shown.

33
semester3/numcs/parts/04_linalg/00_intro.tex
View File

@@ -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.

21
semester3/numcs/parts/04_linalg/01_LU.tex
View File

@@ -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

134
semester3/numcs/parts/04_linalg/02_qr.tex
View File

@@ -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}

1
semester3/numcs/parts/04_linalg/03_svd.tex
View File

@@ -0,0 +1 @@
\subsubsection{Singulärwertzerlegung}