[NumCS] Finish up matrix section

semester3/numcs/numcs-summary.tex

@@ -2,7 +2,7 @@
 
 \newcommand{\dir}{../../latex}
 \input{\dir/include.tex}
-\load{recommended}
+\load{full}
 
 \renewcommand{\authorTitle}{Robin Bacher, Janis Hutz\\\url{https://github.com/janishutz/eth-summaries}}
 \renewcommand{\authorHeaders}{Robin Bacher, Janis Hutz}
@@ -14,6 +14,7 @@
 \startDocument
 \usetcolorboxes
 \setcounter{numberingConfig}{3}
+\setcounter{numberSubsections}{1}
 
 %          ╭────────────────────────────────────────────────╮
 %          │                   Title page                   │

semester3/numcs/parts/introduction/matrix-multiplication.tex

@@ -1,9 +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}
-Das Matrixprodukt kann mit dem Strassen Algorithmus mithilfe der Block-Partitionierung in $\tco{n^{\log_2(7)}}$
+\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.

semester3/numcs/parts/introduction/rounding-errors.tex

@@ -44,6 +44,7 @@ Falls jedoch hier die Auswertung von $\text{Im}(f(x_0 + ih))$ nicht exakt ist, s
             & = -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}