diff --git a/latex b/latex index 565b600..4be2f5e 160000 --- a/latex +++ b/latex @@ -1 +1 @@ -Subproject commit 565b600ef0a1c71ceb05e925078f038b55d401e4 +Subproject commit 4be2f5ed0d92c58b0ca00fa475f0e11cd9569246 diff --git a/semester3/numcs/numcs-summary.pdf b/semester3/numcs/numcs-summary.pdf index 3b556e6..f2d71a3 100644 Binary files a/semester3/numcs/numcs-summary.pdf and b/semester3/numcs/numcs-summary.pdf differ diff --git a/semester3/numcs/numcs-summary.tex b/semester3/numcs/numcs-summary.tex index 09a9ef2..9be121d 100644 --- a/semester3/numcs/numcs-summary.tex +++ b/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 │ diff --git a/semester3/numcs/parts/introduction/matrix-multiplication.tex b/semester3/numcs/parts/introduction/matrix-multiplication.tex index 705f412..6f53e59 100644 --- a/semester3/numcs/parts/introduction/matrix-multiplication.tex +++ b/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}$ \\ + 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. diff --git a/semester3/numcs/parts/introduction/rounding-errors.tex b/semester3/numcs/parts/introduction/rounding-errors.tex index 5fe92e4..6a0c869 100644 --- a/semester3/numcs/parts/introduction/rounding-errors.tex +++ b/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}