% ┌                                                ┐
% │     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.
