mirror of
https://github.com/janishutz/eth-summaries.git
synced 2025-11-25 18:44:24 +00:00
107 lines
6.8 KiB
TeX
107 lines
6.8 KiB
TeX
% ┌ ┐
|
|
% │ AUTHOR: Janis Hutz<info@janishutz.com> │
|
|
% └ ┘
|
|
|
|
\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{code}{python}
|
|
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{code}
|
|
|
|
|
|
\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{code}{python}
|
|
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{code}
|
|
|
|
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.
|