eth-summaries/semester3/numcs/parts/04_linalg/02_qr.tex

\subsubsection{QR-Zerlegung}
Wir können eine Matrix $A \in \R^{m \times n}$ mit $m \geq n$ als $A = QR$ zerlegen, wobei $Q \in \R^{m \times n}$
orthonormale Spalten besitzt und $R \in \R^{n \times n}$ ist eine obere Dreiecksmatrix.

Die QR-Zerlegung ist numerisch stabiler bei schlecht konditionierten Problemen, ist die Basis vieler Eigenwertverfahren und ist ideal für Least Squares.

Wir können mit der QR-Zerlegung auch LGS $Ax = b$ lösen:
\rmvspace
\begin{align*}
    Ax = b \Longleftrightarrow QRx = b \Longleftrightarrow Rx = Q^\top b
\end{align*}

Da $Q$ orthogonal ist, haben wir $Q^{-1} = Q^\top$.
Um das Ganze einfacher zu machen, lösen wir das System $Rx = y$, wobei $y := Q^\top b$.
In numpy können wir direkt mit \texttt{np.linalg.solve()} dies lösen (nutzt automatisch Rückwärtssubstitution)


\newpage
\bg{purple}{Givens-Rotations}

Bei der Givens-Rotation generiert man eine Rotationsmatrix, die die $(i, j)$-Ebene um einen Winkel $\theta$ rotiert.
Die dazu konstruierte Matrix hat dabei die folgende Form (rechts eine Beschreibung des Eintrags $(k, l)$):
\rmvspace
\begin{align*}
    G(i, j, \theta) =
    \begin{bmatrix}
        1      & 0      &        & \cdots &        &        & 0      \\
        0      & \ddots &        &        &        &        &        \\
               &        & c      & \cdots & s      &        &        \\
        \vdots &        & \vdots &        & \vdots &        & \vdots \\
               &        & -s     & \cdots & c      &        &        \\
        0      &        &        &        &        & \ddots &        \\
        0      &        &        & \cdots &        &        & 1      \\
    \end{bmatrix}
    \text{ oder }
    G(i, j, \theta)_{k, l} =
    \begin{cases}
        k = l \land k \neq i, j        & 1  \\
        k = l \land (k = i \lor k = j) & c  \\
        k = i \land l = j              & -s \\
        k = j \land l = i              & s  \\
        \text{sonst}                   & 0
    \end{cases}
\end{align*}

\drmvspace
Dabei ist $c = \cos(\theta)$ und $s = \sin(\theta)$.
Diese Matrix hat einige nützliche Eigenschaften: $G^\top G = I$ (also ist $G$ orthogonal), also gilt auch $G^{-1} = G^\top$ und $G$ modifiziert nur Zeilen $i$ und $j$

Im Zweidimensionalen Raum können wir die Werte für $c$ und $s$ so bestimmen:
\rmvspace
\begin{align*}
    \begin{bmatrix}
        c  & s \\
        -s & c
    \end{bmatrix}
    \begin{bmatrix}
        a \\ b
    \end{bmatrix}
    =
    \begin{bmatrix}
        r \\ 0
    \end{bmatrix} \\
    r = \sqrt{a^2 + b^2}, \smallhspace c = \frac{a}{r}, \smallhspace s = \frac{b}{r}
\end{align*}
Wir haben jedoch das Problem, dass die Berechnung von $r$ überlaufen kann. Dies lösen wir, indem wir skalieren:
\drmvspace
\begin{multicols}{2}
    Falls $|b| > |a|$:
    \rmvspace
    \begin{align*}
        t = \frac{a}{b}, \smallhspace s = \frac{1}{\sqrt{1 + t^2}}, \smallhspace c = s \cdot t
    \end{align*}

    Falls $|a| \geq |b|$:
    \drmvspace
    \begin{align*}
        t = \frac{b}{a}, \smallhspace s = c \cdot t, \smallhspace c = \frac{1}{\sqrt{1 + t^2}}
    \end{align*}
\end{multicols}
Es ist wichtig, dass wir das $r = \text{sign}(a) \sqrt{a^2 + b^2}$ mit Vorzeichen berechnen, um Auslöschung zu vermeiden

Man kann nun mit der Givens-Rotation die QR-Zerlegung durchführen:
\begin{algorithm}
    \begin{spacing}{1.2}
        \caption{\textsc{GivensQRDecomposition}(A)}
        \begin{algorithmic}[1]
            \State $m \gets \texttt{A.shape[0]}$
            \State $n \gets \texttt{A.shape[1]}$
            \State $q \gets$ Initialisiere ein $n \times n$ array
            \For{$j = 1, 2, \ldots, n$}
                \For{$i = m, \ldots, 2$}
                \State Nullsetze $a_{m, j}, a_{m - 1, j}, \ldots, a_{2, j}$ durch Givens-Rotationen $G_{m, j}, G_{m - 1, j}, \ldots, G_{2, j}$
                \EndFor
            \EndFor
            \State \Return $l$
        \end{algorithmic}
    \end{spacing}
\end{algorithm}


\newpage
\bg{purple}{Gram-Schmidt}

Die Idee des Gram-Schmidt-Algorithmus ist es, orthonormale Vektoren zu konstruieren und diese dann zur Matrix $Q$ zusammenzubasteln.

Es wurden zwei Algorithmen behandelt, beide unten in Pseudocode:
\begin{algorithm}
    \begin{spacing}{1.2}
        \caption{\textsc{classicalGramSchmidt}(A)}
        \begin{algorithmic}[1]
            \State $n \gets \texttt{A.shape[0]}$
            \State $q \gets$ Initialisiere ein $n \times n$ array
            \State $r \gets$ Initialisiere ein $n \times n$ array
            \For{$k = 1, 2, \ldots, n$}
                \State $v_k \gets a_k$ \Comment{Der $k$-te Spaltenvektor}
                \For{$i = 1, \ldots, k - 1$}
                    \State $r_{ik} \gets q_i^\top a_k$
                    \State $v_k \gets a_k - \sum_{i = 1}^{k - 1} r_{ik} q_i$
                \EndFor
                \State $r_{kk} = ||v_k||$
                \State $q_k = v_k / r_{kk}$ \Comment{Vektor normieren}
            \EndFor
            \State \Return $q$, $r$
        \end{algorithmic}
    \end{spacing}
\end{algorithm}

\drmvspace
\begin{algorithm}
    \begin{spacing}{1.2}
        \caption{\textsc{modifiedGramSchmidt}(A)}
        \begin{algorithmic}[1]
            \State $n \gets \texttt{A.shape[0]}$
            \State $q \gets$ Initialisiere ein $n \times n$ array
            \State $r \gets$ Initialisiere ein $n \times n$ array
            \For{$k = 1, 2, \ldots, n$}
                \State $v_k \gets a_k$ \Comment{Der $k$-te Spaltenvektor}
                \For{$i = 1, \ldots, k - 1$}
                    \State $r_{ik} \gets q_i^\top v_k$
                    \State $v_k \gets v_k - r_{ik} q_i$
                \EndFor
                \State $q_k = v_k / r_{kk}$ \Comment{Vektor normieren}
            \EndFor
            \State \Return $q$, $r$
        \end{algorithmic}
    \end{spacing}
\end{algorithm}
Falls $R$ nicht benötigt wird, kann viel Speicher gespart werden, indem man das $r_{ik}$ als eine scoped variable verwendet.


\newpage
\bg{purple}{Householder-Reflektor}

Wir konstruieren eine Matrix $H = I - 2 \displaystyle\frac{vv^\top}{v^\top v} = I - 2uu^\top \text{ mit } u = \frac{v}{||v||}$.

Dabei ist die Matrix $H$ orthogonal ($H^\top H = I$) und symmetrisch ($H^\top = H$).

Um nun die QR-Zerlegung durchzuführen mit der Householder-Reflektion fehlen uns die Householder-Reflektoren.
Um diese zu erstellen wollen wir das $v$ so wählen, dass $Hx = ||x|| e_1$ gilt. So werden also $m - 1$ Elemente auf einmal auf Null gesetzt.

Der Ansatz dazu ist entsprechend $v = x - \alpha e_1$ mit $\alpha = -\text{sign}(x_1)||x||$ (minus, um numerische Stabilität zu erhalten) und wir haben dann:
\drmvspace
\begin{align*}
    Hx = \alpha e_1 \Longleftrightarrow Hx = x - 2 \frac{v^\top x}{v^\top v} v
\end{align*}

\drmvspace
Dann müssen wir nur noch $v^\top x$ und $v^\top v$ berechnen und auflösen.
Der vollständige QR-Algorithmus lautet:
\begin{algorithm}
    \begin{spacing}{1.2}
        \caption{\textsc{HouseholderQR}(A)}
        \begin{algorithmic}[1]
            \State $n \gets \texttt{A.shape[0]}$
            \State $H \gets$ Initialisiere ein $n \times n$ array
            \For{$k = 1, 2, \ldots, n$}
                \State $x \gets A(k : m, k$ \Comment{Wähle subvektor der $k$-ten Spalte}
                \State $H_k \gets$ Konstruiere Householder-Reflektor für $x$
                \State $A(k : m, k : n) \gets H_k A(k : m, k : n)$ \Comment{Update}
            \EndFor
            \State $Q \gets H_1 H_2 \cdots H_n$
            \State $R \gets H_n H_{n - 1} \cdots H_1 A$
            \State \Return $Q$, $R$
        \end{algorithmic}
    \end{spacing}
\end{algorithm}

\drmvspace
Die Laufzeiten der verschiedenen Methoden im Vergleich:
\begin{itemize}[noitemsep]
    \item Householder-QR: $\approx 2mn^2$ Flops
    \item Gram-Schmidt: $\approx 2mn^2$ Flops
    \item Givens: $\approx 3mn^2$ Flops
\end{itemize}
Jedoch ist die Householder-Methode bedeutend stabiler als die anderen beiden.