[NumCS] Restructure

semester3/numcs/parts/00_introduction/00_rounding-errors.tex

@@ -0,0 +1,64 @@
+\subsection{Rundungsfehler}
+
+\begin{definition}[]{Absoluter \& Relativer Fehler}
+    \begin{multicols}{2}
+        \begin{itemize}
+            \item \bi{Absoluter Fehler}: $||\tilde{x} - x||$
+            \item \bi{Relativer Fehler}: $\displaystyle \frac{||\tilde{x} - x||}{||x||}$ für $||x|| \neq 0$
+        \end{itemize}
+    \end{multicols}
+    wobei $\tilde{x}$ eine Approximation an $x \in \R$ ist
+\end{definition}
+
+Rundungsfehler entstehen durch die (verhältnismässig) geringe Präzision die man mit der Darstellung von Zahlen auf Computern erreichen kann.
+Zusätzlich kommt hinzu, dass durch Unterläufe (in diesem Kurs ist dies eine Zahl die zwischen $0$ und der kleinsten darstellbaren, positiven Zahl liegt) Präzision verloren gehen kann.
+
+Überläufe hingegen sind konventionell definiert, also eine Zahl, die zu gross ist und nicht mehr dargestellt werden kann.
+
+
+\setcounter{all}{9}
+\begin{remark}[]{Auslöschung}
+    Bei der Subtraktion von zwei ähnlich grossen Zahlen kann es zu einer Addition der Fehler der beiden Zahlen kommen, was dann den relativen Fehler um einen sehr grossen Faktor vergrössert.
+    Die Subtraktion selbst hat einen vernachlässigbaren Fehler
+\end{remark}
+
+\setcounter{all}{18}
+\fancyex{Ableitung mit imaginärem Schritt} Als Referenz in Graphen wird hier oftmals die Implementation des Differenzialquotienten verwendet.
+
+Der Trick hier ist, dass wir mit Komplexen Zahlen in der Taylor-Approximation einer glatten Funktion in $x_0$ einen rein imaginären Schritt durchführen können:
+\begin{align*}
+    f(x_0 + ih) = f(x_0) + f'(x_0)ih - \frac{1}{2} f''(x_0)h^2 - iC \cdot h^3 \text{ für } h \in \R \text{ und } h \rightarrow 0
+\end{align*}
+Da $f(x_0)$ und $f''(x_0)h^2$ reell sind, verschwinden die Terme, wenn wir nur den Imaginärteil des Ausdruckes weiterverwenden. Nach weiteren Vereinfachungen und Umwandlungen erhalten wir
+\begin{align*}
+    f'(x_0) \approx \frac{\text{Im}(f(x_0 + ih))}{h}
+\end{align*}
+
+Falls jedoch hier die Auswertung von $\text{Im}(f(x_0 + ih))$ nicht exakt ist, so kann der Fehler beträchtlich sein.
+
+
+\setcounter{all}{20}
+\fancyex{Konvergenzbeschleunigung nach Richardson}
+\begin{align*}
+    y f'(x) & = y d\left(\frac{h}{2}\right) + \frac{1}{6} f'''(x) h^2 + \frac{1}{480}f^{(s)} h^4 + \ldots - f'(x) \\
+            & = -d(h) - \frac{1}{6}f'''(x) h^2 + \frac{1}{120} f^{(s)}(x) h^n \Leftrightarrow 3 f'(x)             \\
+            & = 4 d\left(\frac{h}{2}\right)  d(h) + \tco{h^4} \Leftrightarrow
+\end{align*}
+% TODO: Need to finish with notes from the exercise sessions
+
+
+\fhlc{Cyan}{Schema}
+
+\begin{align*}
+    d(h) = \frac{f(x + h) - f(x - h)}{2h}
+\end{align*}
+
+wobei im Schema dann
+\begin{align*}
+    R_{l, 0} = d\left( \frac{h}{2^l} \right)
+\end{align*}
+und
+\begin{align*}
+    R_{l, k} = \frac{4^k \cdot R_{l, k - 1} - R_{l - 1, k - 1}}{4^k - 1}
+\end{align*}
+und $f'(x) = R_{l, k} + C \cdot \left( \frac{h}{2^l} \right)^{2k + 2}$

semester3/numcs/parts/00_introduction/01_time-complexity.tex

@@ -0,0 +1,10 @@
+\newsection
+\subsection{Rechenaufwand}
+In NumCS wird die Anzahl elementarer Operationen wie Addition, Multiplikation, etc benutzt, um den Rechenaufwand zu beschreiben. 
+Wie in Algorithmen und * ist auch hier wieder $\tco{\ldots}$ der Worst Case.
+Teilweise werden auch andere Funktionen wie $\sin, \cos, \sqrt{\ldots}, \ldots$ dazu gezählt.
+
+Die Basic Linear Algebra Subprograms (= BLAS), also grundlegende Operationen der Linearen Algebra, wurden bereits stark optimiert und sollten wann immer möglich verwendet werden und man sollte auf keinen Fall diese selbst implementieren.
+
+Dieser Kurs verwendet \texttt{numpy}, \texttt{scipy}, \texttt{sympy} (collection of implementations for symbolic computations) und \texttt{matplotlib}.
+Dieses Ecosystem ist eine der Stärken von Python und ist interessanterweise zu einem Grossteil nicht in Python geschrieben, da dies sehr langsam wäre.

semester3/numcs/parts/00_introduction/02_matrix-multiplication.tex

@@ -0,0 +1,106 @@
+\newsection
+\subsection{Rechnen mit Matrizen}
+Wie in Lineare Algebra besprochen, ist das Resultat der Multiplikation einer Matrix $A \in \C^{m \times n}$ und einer Matrix $B \in \C^{n \times p}$ ist eine Matrix $AB = \in \C^{m \times p}$
+
+\fhlc{Cyan}{In NumPy} haben wir folgende Funktionen:
+\begin{itemize}
+    \item \verb|b @ a| (oder \verb|np.dot(b, a)| oder \verb|np.einsum('i,i', b, a)| für das Skalarprodukt
+    \item \verb|A @ B| (oder \verb|np.einsum('ik,kj->ij', )|) für das Matrixprodukt
+    \item \verb|A @ x| (oder \verb|np.einsum('ij,j->i', A, x)|) für Matrix $\times$ Vektor
+    \item \verb|A.T| für die Transponierung
+    \item \verb|A.conj()| für die komplexe Konjugation (kombiniert mit \verb|.T| = Hermitian Transpose)
+    \item \verb|np.kron(A, B)| für das Kroneker Produkt
+    \item \verb|b = np.array([4.j, 5.j])| um einen Array mit komplexen Zahlen zu erstellen (\verb|j| ist die imaginäre Einheit, aber es muss eine Zahl direkt daran geschrieben werden)
+\end{itemize}
+
+
+\setcounter{all}{4}
+\fancyremark{Rang der Matrixmultiplikation} $\text{Rang}(AX) = \min(\text{Rang}(A), \text{Rang}(X))$
+
+\setcounter{all}{7}
+\fancyremark{Multiplikation mit Diagonalmatrix $D$} $D \times A$ skaliert die Zeilen von $A$ während $A \times D$ die Spalten skaliert
+
+\stepcounter{all}
+\inlineex \verb|D @ A| braucht $\tco{n^3}$ Operationen, wenn wir jedoch \verb|D.diagonal()[:, np.newaxis] * A| verwenden, so haben wir nur noch $\tco{n^2}$ Operationen, da wir die vorige Bemerkung Nutzen und also nur noch eine Skalierung vornehmen.
+So können wir also eine ganze Menge an Speicherzugriffen sparen, was das Ganze bedeutend effizienter macht
+
+\setcounter{all}{14}
+\inlineremark Wir können bestimmte Zeilen oder Spalten einer Matrix skalieren, in dem wir einer Identitätsmatrix im unteren Dreieck ein Element hinzufügen.
+Wenn wir nun diese Matrix $E$ (wie die in der $LU$-Zerlegung) linksseitig mit der Matrix $A$ multiplizieren (bspw. $E^{(2, 1)}A$), dann wird die zugehörige Zeile skaliert.
+Falls wir aber $AE^{(2, 1)}$ berechnen, so skalieren wir die Spalte
+
+\fancyremark{Blockweise Berechnung} Man kann das Matrixprodukt auch Blockweise berechnen.
+Dazu benutzen wir eine Matrix, deren Elemente andere Matrizen sind, um grössere Matrizen zu generieren.
+Die Matrixmultiplikation funktioniert dann genau gleich, nur dass wir für die Elemente Matrizen und nicht Skalare haben.
+
+% ────────────────────────────────────────────────────────────────────
+\hspace{1mm}
+\hrule
+\hspace{1mm}
+Untenstehend eine Tabelle zum Vergleich der Operationen auf Matrizen
+
+\begin{tables}{lcccc}{Name           & Operation & Mult  & Add         & Komplexität}
+              Skalarprodukt          & $x^H y$   & $n$   & $n - 1$     & $\tco{n}$    \\
+              Tensorprodukt          & $x y^H$   & $nm$  & $0$         & $\tco{mn}$   \\
+              Matrix $\times$ Vektor & $Ax$      & $mn$  & $(n - 1)m$  & $\tco{mn}$   \\
+              Matrixprodukt          & $AB$      & $mnp$ & $(n - 1)mp$ & $\tco{mnp}$  \\
+\end{tables}
+\inlineremark Das Matrixprodukt kann mit Strassen's Algorithmus mithilfe der Block-Partitionierung in $\tco{n^{\log_2(7)}} \approx \tco{n^{2.81}}$ berechnet werden.
+
+\fancyremark{Rang 1 Matrizen} Können als Tensorprodukt von zwei Vektoren geschrieben werden.
+Dies ist beispielsweise hierzu nützlich:
+
+Sei $A = ab^\top$. Dann gilt $y = Ax \Leftrightarrow y = a(b^\top x)$, was dasselbe Resultat ergibt, aber nur $\tco{m + n}$ Operationen und nicht $\tco{mn}$ benötigt wie links.
+
+\inlineex Für zwei Matrizen $A, B \in \R^{n \times p}$ mit geringem Rang $p \ll n$, dann kann mithilfe eines Tricks die Rechenzeit von \verb|np.triu(A @ B.T) @ x| von $\tco{pn^2}$ auf $\tco{pn}$ reduziert werden.
+Die hier beschriebene Operation berechnet $\text{Upper}(AB^\top) x$ wobei $\text{Upper}(X)$ das obere Dreieck der Matrix $X$ zurück gibt.
+Wir nennen diese Matrix hier $R$.
+Wir können in NumPy den folgenden Ansatz verwenden, um die Laufzeit zu verringern:
+Da die Matrix $R$ eine obere Dreiecksmatrix ist, ist das Ergebnis die Teilsummen von unserem Umgekehrten Vektor $x$, also können wir mit \verb|np.cumsum(x[::-1], axis=0)[::-1]| die Kummulative Summe berechnen.
+Das \verb|[::-1]| dient hier lediglich dazu, den Vektor $x$ umzudrehen, sodass das richtige Resultat entsteht.
+Die vollständige Implementation sieht so aus:
+\begin{minted}{python}
+import numpy as np
+
+def low_rank_matrix_vector_product(A: np.ndarray, B: np.ndarray, x: np.ndarray):
+    n, _ = A.shape
+    y = np.zeros(n)
+
+    # Compute B * x with broadcasting (x needs to be reshaped to 2D)
+    v = B * x[:, None]
+
+    # s is defined as the reverse cummulative sum of our vector
+    # (and we need it reversed again for the final calculation to be correct)
+    s = np.cumsum(v[::-1], axis=0)[::-1]
+
+    y = np.sum(A * s)
+\end{minted}
+
+
+\setcounter{all}{21}
+\fancydef{Kronecker-Produkt} Das Kronecker-Produkt ist eine $(ml) \times (nk)$-Matrix, für $A \in \R^{m \times n}$ und $B \in \R^{l \times k}$, die wir in NumPy einfach mit \verb|np.kron(A, B)| berechnen können (ist jedoch nicht immer ideal):
+\begin{align*}
+    A \otimes B :=
+    \begin{bmatrix}
+        (A)_{1, 1} B & (A)_{1, 2}B & \ldots & \ldots & (A)_{1, n} B \\
+        (A)_{2, 1} B & (A)_{2, 2}B & \ldots & \ldots & (A)_{2, n} B \\
+        \vdots       & \vdots      & \ddots & \ddots & \vdots       \\
+        (A)_{m, 1} B & (A)_{m, 2}B & \ldots & \ldots & (A)_{m, n} B \\
+    \end{bmatrix}
+\end{align*}
+
+\fancyex{Multiplikation des Kronecker-Produkts mit Vektor} Wenn man $A \otimes B \cdot x$ berechnet, so ist die Laufzeit $\tco{m \times n \times l \times k}$, aber wenn wir den Vektor $x$ in $n$ gleich grosse Blöcke aufteilen (was man je nach gewünschter nachfolgender Operation in NumPy in $\tco{1}$ machen kann mit \verb|x.reshape(n, x.shape[0] / n)|), dann ist es möglich das Ganze in $\tco{m \cdot l \cdot k}$ zu berechnen. 
+
+Die vollständige Implementation ist auch hier nicht schwer und sieht folgendermassen aus:
+\begin{minted}{python}
+import numpy as np
+
+def fast_kron_vector_product(A: np.ndarray, B: np.ndarray, x: np.ndarray):
+    # First multiply Bx_i, (and define x_i as a reshaped numpy array to save cost (as that will create a valid array))
+    # This will actually crash if x.shape[0] is not divisible by A.shape[0]
+    bx = B * x.reshape(A.shape[0], round(x.shape[0] / A.shape[0]))
+    # Then multiply a with the resulting vector
+    y = A * bx
+\end{minted}
+
+Um die oben erwähnte Laufzeit zu erreichen muss erst ein neuer Vektor berechnet werden, oben im Code \verb|bx| genannt, der eine Multiplikation von \verb|Bx_i| als Einträge hat.