mirror of
https://github.com/janishutz/eth-summaries.git
synced 2026-01-11 19:48:27 +00:00
[NumCS] Finish householder and givens
This commit is contained in:
Binary file not shown.
@@ -33,7 +33,6 @@
|
||||
\vspace{4cm}
|
||||
\begin{center}
|
||||
\begin{Large}
|
||||
\quote{Wer in Python Type annotation benötigt, der soll kein Python verwenden} (2025-10-09T10:43) % FIXME: Marked for removal
|
||||
\quote{Wenn ich keine Lust habe, das zu berechnen, dann wende ich einfach Gewalt an}
|
||||
\end{Large}
|
||||
|
||||
@@ -77,17 +76,6 @@ We have also taken some explanations and code examples from the slides of our TA
|
||||
\color{MidnightBlue}\fbox{\href{https://n.ethz.ch/~muellerni/courses/numcs25.php}{here}}\color{black}
|
||||
% TODO: Update this when n.ethz is taken offline completely
|
||||
|
||||
|
||||
To add to the one quote regarding Python and type annotation: This is objectively wrong and a really hot take.
|
||||
Yes, this applies for small projects, but libraries \textit{DO} need type annotation, as you can't possibly read the entire library's code to use it.
|
||||
The reason this quote was even included here is that his coding style is really awful (yes, there were semicolons in his Python-code sometimes)
|
||||
and he rambled about bad coding style for about 10 minutes in this lecture.
|
||||
|
||||
Meanwhile his code has variable names that neither future him, nor anybody else can make much sense of intuitively.
|
||||
You can get away without type annotation in Python, even in larger projects, but only if you give variables proper names!
|
||||
|
||||
Moral of the story: Use descriptive variable names and do NOT use $t$, $tt$, $ttt$, \dots
|
||||
|
||||
% ────────────────────────────────────────────────────────────────────
|
||||
|
||||
% ╭────────────────────────────────────────────────╮
|
||||
|
||||
@@ -1,8 +1,6 @@
|
||||
% ┌ ┐
|
||||
% │ Author: Robin Bacher │
|
||||
% └ ┘
|
||||
% TODO: If you want your email to be in there, note it down here.
|
||||
% I also did not touch the unedited files to avoid conflicts
|
||||
\subsection{Interpolation und Polynome}
|
||||
Bei der Interpolation versuchen wir eine Funktion $\tilde{f}$ durch eine Menge an Datenpunkten einer Funktion $f$ zu finden.\\
|
||||
Die $x_i$ heissen Stützstellen/Knoten, für welche $\tilde{f}(x_i) = y_i$ gelten soll. (Interpolationsbedingung)
|
||||
@@ -30,9 +28,6 @@ Dann lässt sich der Bezug zwischen $f$ und $\tilde{f} = f_n(x)$ so ausdrücken:
|
||||
\inlineremark Unterräume $\mathcal{V}_n$ existieren nicht nur für Polynome, wir beschränken uns aber auf $b_j(x) = x^{i-1}$.
|
||||
Andere Möglichkeiten: $b_j = \cos((j-1)\cos^-1(x))$ \textit{(Chebyshev)} oder $b_j = e^{i2\pi j x}$ \textit{(Trigonometrisch)}
|
||||
|
||||
% FIXME: This could go into a special "maths theory" section -> GOOD
|
||||
|
||||
\setLabelNumber{all}{4}
|
||||
\fancytheorem{Peano} $f$ stetig $\implies \exists p(x)$ welches $f$ in $||\cdot||_\infty$ beliebig gut approximiert.
|
||||
|
||||
\setLabelNumber{all}{5}
|
||||
|
||||
@@ -29,9 +29,7 @@ Das Inverse davon nimmt eine Funktion der Frequenz und transformiert diese in ei
|
||||
%
|
||||
\inlineremark $p_m : \R \rightarrow \C$ ist periodisch mit Periode $1$.
|
||||
Falls $\gamma_{-j} = \overline{\gamma_j}$ für alle $j$, dann ist $p_m$ reellwertig und
|
||||
% NOTE: Uhh... do we want to use the fancy symbols for real and imaginary part or just use $\text{Re}$?
|
||||
% RESPONE: whatever he uses in the script, preferably \text{Re}() etc.
|
||||
$p_m$ kann folgendermassen dargestellt werden ($a_0 = 2\gamma_0, a_j = 2\Re(\gamma_j)$ und $b_j = -2\Im(\gamma_j)$):
|
||||
$p_m$ kann folgendermassen dargestellt werden ($a_0 = 2\gamma_0, a_j = 2\text{Re}(\gamma_j)$ und $b_j = -2\text{Im}(\gamma_j)$):
|
||||
\rmvspace
|
||||
\begin{align*}
|
||||
p_m(t) = \frac{a_0}{2} + \sum_{j = 1}^{m} (a_j \cos(2\pi jt) + b_j \sin(2\pi jt))
|
||||
@@ -123,8 +121,6 @@ Die dargestellte Funktion ist die Fourier-Reihe der charakteristischen Funktion
|
||||
\end{align*}
|
||||
Mit $c = \pi(a + b)$ und $d = \pi(b - a)$
|
||||
|
||||
% TODO: Replace with rendered image from matplotlib (will be higher quality than screenshot from script and can tweak it to our liking)
|
||||
% we will have it anyway after solving the exercises, so might as well
|
||||
\begin{figure}[h!]
|
||||
\begin{center}
|
||||
\includegraphics[width=0.95\textwidth]{assets/01_interpolation/01_trigonometric/overarcing.png}
|
||||
@@ -140,27 +136,3 @@ wobei $l = 0, 1 \ldots, N - 1$ und $N$ die Anzahl der Intervalle ist:
|
||||
\begin{align*}
|
||||
\hat{f}_N(k) := \frac{1}{N} \sum_{l = 0}^{N - 1} f(t_l) e^{-2\pi ikt_l} \approx \hat{f}(k)
|
||||
\end{align*}
|
||||
|
||||
% TODO: Consider if we should use the below
|
||||
|
||||
% \begin{tikzpicture}
|
||||
% \begin{axis}[
|
||||
% legend pos=outer north east,
|
||||
% title=Function plot of $f(x)$ (parts coloured),
|
||||
% axis lines = box,
|
||||
% xlabel = $x$,
|
||||
% ylabel = $y$,
|
||||
% variable = t,
|
||||
% trig format plots = rad,
|
||||
% ]
|
||||
% \addplot [
|
||||
% domain=1:4,
|
||||
% samples=70,
|
||||
% color=blue,
|
||||
% ]
|
||||
% {log2(x)};
|
||||
% \addlegendentry{$ y=x^2 - x - 0.5$}
|
||||
% \end{axis}
|
||||
% \node (0) at (0, 0) {};
|
||||
% \end{tikzpicture}
|
||||
|
||||
|
||||
@@ -18,13 +18,11 @@
|
||||
\numberingOff
|
||||
\inlineex Zur Fehlerbetrachtung verwenden wir drei Funktionen $f : [0, 1] \rightarrow \R$, welche wir mit trigonometrischer Interpolation an den Punkten $\frac{k}{N}$ approximieren:
|
||||
\begin{enumerate}[label=(\Roman*)]
|
||||
% FIXME: Possibly wrong function definition in script
|
||||
\item Stufenfunktion (periodische Fortsetzung von $f$) $f : [0, 1] \rightarrow \R$ mit $f(t) = \begin{cases}
|
||||
0 & \text{für } \left| t - \frac{1}{2} \right| > \frac{1}{4} \\
|
||||
1 & \text{für } \left| t - \frac{1}{2} \right| \leq \frac{1}{4}
|
||||
\end{cases}$
|
||||
\item Periodische, glatte Funktion $h: \R \rightarrow \R$ mit $h(t) = \displaystyle \frac{1}{\sqrt{1 + \frac{1}{2} \sin(2\pi t)}}$
|
||||
% TODO: Is it $h$ or $g$ here? $g$ makes little sense imho
|
||||
\item Hutfunktion (periodische Fortsetzung von $h$) $g : [0, 1] \rightarrow \R$ mit $g(t) = \left| t - \frac{1}{2} \right|$
|
||||
\end{enumerate}
|
||||
Die untenstehende Abbildung \ref{fig:interpolation-error-examples} beinhaltet einen Plot, auf dem die Konvergenz in Abhängigkeit des Grades des Interpolationspolynoms aufgetragen ist.
|
||||
|
||||
@@ -97,7 +97,6 @@ Die zugehörige Quadraturformel ist exakt für Polynome $2s - 1$-ten Grades. Die
|
||||
% Yeah, there is *no way* we would have figured that out... and he even provides a proof for it... yikes...
|
||||
\inlinetheorem Die Gewichte der Gauss-Legendre-Quadraturformel sind positiv.
|
||||
|
||||
% TODO: Consider adding code
|
||||
\begin{tables}{cccc}{Algorithmus & Laufzeit & Genauigkeit Knoten & Genauigkeit Gewichte}
|
||||
GW (1969) & $\tco{s^3}$ / $\tco{s^2}$ & $\tco{1}$ & $\tco{s^2}$ \\
|
||||
Bogaert-Townsend & $\tco{s}$ & $\tco{1}$ & $\tco{1}$ \\
|
||||
|
||||
@@ -17,5 +17,3 @@ Auf Seiten 150 - 151 im Skript findet sich Code, um eine adaptive Quadratur durc
|
||||
Mit \texttt{scipy.integrate.quadrature} können wir die Gauss-Quadratur verwenden.
|
||||
|
||||
Für $x \in \R^d$, also eine mehrdimensionale Funktion der Dimension $d$ können wir \texttt{scipy.integrate.nquad} verwenden. Mehr dazu im nächsten Kapitel
|
||||
|
||||
% TODO: Possibly explain the graphs and / or add code for computation using scipy
|
||||
|
||||
@@ -6,7 +6,6 @@
|
||||
\subsection{Fixpunktiteration}
|
||||
Ein $1$-Punkt-Verfahren benötigt nur den vorigen Wert: $x^{(k + 1)} = \phi(x^{(k)})$
|
||||
|
||||
% FIXME: Below konsistent is probably wrong, but is what is in the script
|
||||
\inlinedef Eine Fixpunktiteration heisst konsistent mit $F(x) = 0$ falls $F(x) = 0 \Leftrightarrow \phi(x) = x$
|
||||
|
||||
\inlineex Für $F(x) = xe^x - 1$ mit $x \in [0, 1]$ liefert $\phi_1(x) = e^{-x}$ lineare Konvergenz,
|
||||
@@ -28,8 +27,6 @@ $\phi_2(x) = \frac{1 + x}{1 + e^x}$ quadratische Konvergenz und $\phi_3(x) = x +
|
||||
Dieser ist der Grenzwert der Folge $x^{(k + 1)} = \phi(x^{(k)})$.
|
||||
\end{theorem}
|
||||
|
||||
% NOTE: If need be, we can switch to theorem here, or I can add a new environment for "support theorem" or the like,
|
||||
% I however feel like a Lemma suits the idea of "Hilfstheorem" quite well
|
||||
\inlinelemma Für $U \subseteq \R^n$ konvex und $\phi : U \rightarrow \R^n$ stetig differenzierbar mit $L := \sup_{x \in U} ||D_\phi(x)|| < 1$
|
||||
($D_\phi(x)$ ist die Jacobi-Matrix von $\phi(x)$).
|
||||
Wenn $\phi(x^*) = x^*$ für $x^* \in U$, dann konvergiert $x^{(k + 1)} = \phi(x^{(k)})$ gegen $x^*$ lokal mindestens linear.
|
||||
|
||||
@@ -43,6 +43,7 @@ Wenn sie zudem normiert sind (also $||q_i||_2 = 1 \smallhspace \forall i \leq n$
|
||||
|
||||
\drmvspace
|
||||
\shade{gray}{Perturbierte LGS}
|
||||
\setLabelNumber{all}{18}
|
||||
|
||||
Statt $Ax = b$ ist das LGS ungenau gegeben: $(A + \Delta A)(\tilde{x} - x) = \Delta b - \Delta Ax$.
|
||||
|
||||
@@ -79,7 +80,6 @@ Passen oft nicht (direkt) in den Speicher: effizientere Speicherung nötig, mög
|
||||
% 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$
|
||||
|
||||
@@ -80,6 +80,24 @@ Wir haben jedoch das Problem, dass die Berechnung von $r$ überlaufen kann. Dies
|
||||
\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
|
||||
|
||||
Man kann nun mit der Givens-Rotation die QR-Zerlegung durchführen:
|
||||
\begin{algorithm}
|
||||
\begin{spacing}{1.2}
|
||||
\caption{\textsc{GivensQRDecomposition}(A)}
|
||||
\begin{algorithmic}[1]
|
||||
\State $m \gets \texttt{A.shape[0]}$
|
||||
\State $n \gets \texttt{A.shape[1]}$
|
||||
\State $q \gets$ Initialisiere ein $n \times n$ array
|
||||
\For{$j = 1, 2, \ldots, n$}
|
||||
\For{$i = m, \ldots, 2$}
|
||||
\State Nullsetze $a_{m, j}, a_{m - 1, j}, \ldots, a_{2, j}$ durch Givens-Rotationen $G_{m, j}, G_{m - 1, j}, \ldots, G_{2, j}$
|
||||
\EndFor
|
||||
\EndFor
|
||||
\State \Return $l$
|
||||
\end{algorithmic}
|
||||
\end{spacing}
|
||||
\end{algorithm}
|
||||
|
||||
|
||||
\newpage
|
||||
\bg{purple}{Gram-Schmidt}
|
||||
@@ -131,4 +149,47 @@ Es wurden zwei Algorithmen behandelt, beide unten in Pseudocode:
|
||||
Falls $R$ nicht benötigt wird, kann viel Speicher gespart werden, indem man das $r_{ik}$ als eine scoped variable verwendet.
|
||||
|
||||
|
||||
\newpage
|
||||
\bg{purple}{Householder-Reflektor}
|
||||
|
||||
Wir konstruieren eine Matrix $H = I - 2 \displaystyle\frac{vv^\top}{v^\top v} = I - 2uu^\top \text{ mit } u = \frac{v}{||v||}$.
|
||||
|
||||
Dabei ist die Matrix $H$ orthogonal ($H^\top H = I$) und symmetrisch ($H^\top = H$).
|
||||
|
||||
Um nun die QR-Zerlegung durchzuführen mit der Householder-Reflektion fehlen uns die Householder-Reflektoren.
|
||||
Um diese zu erstellen wollen wir das $v$ so wählen, dass $Hx = ||x|| e_1$ gilt. So werden also $m - 1$ Elemente auf einmal auf Null gesetzt.
|
||||
|
||||
Der Ansatz dazu ist entsprechend $v = x - \alpha e_1$ mit $\alpha = -\text{sign}(x_1)||x||$ (minus, um numerische Stabilität zu erhalten) und wir haben dann:
|
||||
\drmvspace
|
||||
\begin{align*}
|
||||
Hx = \alpha e_1 \Longleftrightarrow Hx = x - 2 \frac{v^\top x}{v^\top v} v
|
||||
\end{align*}
|
||||
|
||||
\drmvspace
|
||||
Dann müssen wir nur noch $v^\top x$ und $v^\top v$ berechnen und auflösen.
|
||||
Der vollständige QR-Algorithmus lautet:
|
||||
\begin{algorithm}
|
||||
\begin{spacing}{1.2}
|
||||
\caption{\textsc{HouseholderQR}(A)}
|
||||
\begin{algorithmic}[1]
|
||||
\State $n \gets \texttt{A.shape[0]}$
|
||||
\State $H \gets$ Initialisiere ein $n \times n$ array
|
||||
\For{$k = 1, 2, \ldots, n$}
|
||||
\State $x \gets A(k : m, k$ \Comment{Wähle subvektor der $k$-ten Spalte}
|
||||
\State $H_k \gets$ Konstruiere Householder-Reflektor für $x$
|
||||
\State $A(k : m, k : n) \gets H_k A(k : m, k : n)$ \Comment{Update}
|
||||
\EndFor
|
||||
\State $Q = H_1 H_2 \cdots H_n$
|
||||
\State $R = H_n H_{n - 1} \cdots H_1 A$
|
||||
\State \Return $Q$, $R$
|
||||
\end{algorithmic}
|
||||
\end{spacing}
|
||||
\end{algorithm}
|
||||
|
||||
Die Laufzeiten der verschiedenen Methoden im Vergleich:
|
||||
\begin{itemize}
|
||||
\item Householder-QR: $\approx 2mn^2$ Flops
|
||||
\item Gram-Schmidt: $\approx 2mn^2$ Flops
|
||||
\item Givens: $\approx 3mn^2$ Flops
|
||||
\end{itemize}
|
||||
Jedoch ist die Householder-Methode bedeutend stabiler als die anderen beiden.
|
||||
|
||||
@@ -37,6 +37,4 @@ $B_a =
|
||||
A^H & 0
|
||||
\end{bmatrix}
|
||||
$
|
||||
% TODO: What the f does the kappa mean???
|
||||
ersetzen, wobei wir $a$ so wählen, dass $\kappa(B_a)$ minimal wird.
|
||||
% TODO: Consider adding the code and the example, but likely skip
|
||||
ersetzen, wobei wir $a$ so wählen, dass $\kappa(B_a)$ minimal wird (Zur Erinnerung, $\kappa$ ist die Konditionszahl der Matrix).
|
||||
|
||||
Reference in New Issue
Block a user