Matrizen mit dem QR Verfahren zerlegen

Existenzsatz der Zerlegung mit QR

Die $\mathrm{QR}$ -Zerlegung faktorisiert eine Matrix in einen orthogonalen und einen oberen dreieckigen Faktor. $Orthogonale Multiplikationen erhalten Längen und Winkel$ und eignen sich deshalb für stabile numerische Rechnungen.

Sei $A\in\mathbb{R}^{m\times n}$ mit $m\geq n$ und $\operatorname{rang}(A)=n$ . Dann gibt es eine orthogonale Matrix $Q\in\mathbb{R}^{m\times m}$ mit $Q^TQ=I$ und eine obere Dreiecksmatrix $R\in\mathbb{R}^{m\times n}$ , deren erste $n$ $Diagonaleinträge$ positiv sind, sodass

A = Q \cdot R

Diese Faktorisierung heißt $\mathrm{QR}$ -Zerlegung von $A$ .

Beweis mit Orthogonalisierung nach Gram und Schmidt

Sei $a_j$ die $j$ -te Spalte von $A$ . Das Gram-Schmidt-Verfahren $entfernt die Anteile$ in bereits bestimmten Richtungen und normiert anschließend den Rest:

\tilde{q}_j := a_j - \sum_{k=1}^{j-1} \langle a_j, q_k \rangle \cdot q_k

Der volle Spaltenrang garantiert $\widetilde q_j\neq0$ in jedem Schritt. Die normierten Vektoren bilden die Spalten von $Q_1\in\mathbb{R}^{m\times n}$ und erfüllen $Q_1^TQ_1=I_n$ . Jede ursprüngliche Spalte lässt sich darstellen als

a_j = \sum_{k=1}^{j-1} \langle a_j, q_k \rangle \cdot q_k + \tilde{q}_j

Damit gilt $A=Q_1R_1$ , wobei $R_1\in\mathbb{R}^{n\times n}$ eine obere Dreiecksmatrix mit $r_{ii}=\|\widetilde q_i\|_2>0$ ist.

Ergänzt man $Q_1$ um $m-n$ orthonormale Spalten, entsteht $Q=(Q_1\quad Q_2)\in\mathbb{R}^{m\times m}$ .

Durch Anhängen von $m-n$ Nullzeilen an $R_1$ entsteht:

R=\begin{pmatrix}R_1\\0\end{pmatrix}\in\mathbb{R}^{m\times n}

Die ergänzten Matrizen erfüllen weiterhin $A=Q_1R_1=QR$ .

Vollständige und reduzierte Zerlegung mit QR

Für $A\in\mathbb{R}^{m\times n}$ mit $m\geq n$ beschreiben die vollständige und die reduzierte Form dieselbe Faktorisierung mit $unterschiedlichen Matrixgrößen$ . Beide ergeben beim Zurückmultiplizieren wieder $A$ , sodass die Wahl davon abhängt, ob du die vollständige Orthonormalbasis brauchst oder nur die Spalten, die die Daten tragen.

Vollständige Zerlegung mit QR

Die vollständige Form $behält alle$ $m$ Spalten von $Q$ :

A = Q \cdot R

Dabei ist $Q\in\mathbb{R}^{m\times m}$ orthogonal. Die Matrix $R\in\mathbb{R}^{m\times n}$ enthält oben einen dreieckigen $n\times n$ -Block und darunter $m-n$ Nullzeilen.

Reduzierte Zerlegung mit QR

Die reduzierte Form $entfernt diese Nullzeilen$ und die dazugehörigen ungenutzten Spalten von $Q$ :

\begin{aligned} A &= Q \cdot R \\ &= (Q_1 \quad Q_2) \cdot \begin{pmatrix} R_1 \\ 0 \end{pmatrix} \\ &= Q_1 \cdot R_1 \end{aligned}

Hier enthält $Q_1\in\mathbb{R}^{m\times n}$ die ersten $n$ Spalten von $Q$ , und $R_1\in\mathbb{R}^{n\times n}$ ist eine obere Dreiecksmatrix. Das Speichern von $Q_1$ spart Speicher und vermeidet Rechnungen mit Blöcken, die bereits als null bekannt sind.

Die ungenutzten Spalten $Q_2\in\mathbb{R}^{m\times(m-n)}$ bilden eine Orthonormalbasis des Kerns von $A^T$ :

A^T Q_2 = R_1^T Q_1^T Q_2 = 0

Der Kern von $A^T$ enthält genau die Vektoren $x$ , die $A^T x = 0$ erfüllen.

Eindeutigkeit

Für eine Matrix mit vollem Spaltenrang ist die reduzierte Faktorisierung $A=Q_1R_1$ eindeutig, wenn alle $Diagonaleinträge$ von $R_1$ $positiv$ sind. Diese Vorzeichenbedingung beseitigt die verbleibende Freiheit, denn multiplizierst du eine Spalte von $Q_1$ und die zugehörige Zeile von $R_1$ mit $-1$ , bleibt das Produkt unverändert.

Beispiel einer dünnen Zerlegung mit QR

Sei

A=\begin{pmatrix}1&1\\1&0\\0&1\end{pmatrix}=(a_1\ \,a_2)

$Gram-Schmidt liefert$ zunächst

q_1=\frac{a_1}{\|a_1\|_2}=\frac1{\sqrt2}\begin{pmatrix}1\\1\\0\end{pmatrix}

Damit gilt

Q_1=\begin{pmatrix}1/\sqrt2&1/\sqrt6\\1/\sqrt2&-1/\sqrt6\\0&2/\sqrt6\end{pmatrix}

Householderverfahren zur Zerlegung mit QR

Das klassische Gram-Schmidt-Verfahren lässt sich direkt herleiten, kann aber bei $nahezu linear abhängigen Spalten$ in Gleitkommaarithmetik Orthogonalität verlieren. Das modifizierte Verfahren verringert diesen Fehler. Numerische Software verwendet für dichte Matrizen meist Householder-Spiegelungen. Givens-Rotationen eignen sich, wenn nur einzelne Einträge beseitigt werden sollen.

Householdertransformation

Für einen Einheitsvektor $v\in\mathbb{R}^m$ ist die Householder-Matrix definiert durch

S = I - 2vv^T \in \mathbb{R}^{m \times m}

Das $äußere Produkt$ $vv^T$ multipliziert eine $m\times1$ -Spalte mit einer $1\times m$ -Zeile und ergibt daher eine $m\times m$ -Matrix vom Rang eins. Das umgekehrte Produkt $v^Tv$ ist der Skalar $1$ .

Eigenschaften der Householdertransformation

Für $S=I-2vv^T$ und $\|v\|_2=1$ gelten folgende Eigenschaften. Symmetrie und Orthogonalität folgen aus dem Einheitsvektor in der Definition, die Spiegelung erklärt die Wirkung auf einen Vektor, und $Konditionszahl eins$ bedeutet, dass eine exakte Spiegelung die euklidische Norm eines Störvektors erhält.

$S$ ist symmetrisch: $S^T = S$
$S$ ist eine orthogonale Matrix: $S^T S = I$
$Sx$ spiegelt $x\in\mathbb{R}^m$ an der zu $v$ senkrechten Hyperebene
$\text{cond}_2(S) = 1$

Householderverfahren

Sei $A\in\mathbb{R}^{m\times n}$ mit $m\geq n$ und $\operatorname{rang}(A)=n$ . $In jedem Householder-Schritt$ werden die Einträge unter einer Diagonalposition null, während die bereits erzeugten Nullen erhalten bleiben.

Wir beginnen mit $A_1=A$ und spiegeln die erste Spalte auf ein Vielfaches des ersten Koordinatenvektors:

\tilde{a}_1 \text{ erste Spalte von } A_1

u_1 = \tilde{a}_1 \mp \|\tilde{a}_1\|_2 \cdot e_1

Iterativer Transformationsprozess

Mit $S_1=I_m-2v_1v_1^T$ werden alle Einträge der ersten Spalte unterhalb der Diagonale null:

A_2 = S_1 \cdot A_1 = \begin{pmatrix} r_{11} & * \\ 0 & \tilde{A}_2 \end{pmatrix}

In der Matrix $A_2$ gilt $r_{11} = \pm\|\tilde{a}_1\|_2$ und $\tilde{A}_2 \in \mathbb{R}^{m-1 \times n-1}$ .

Im nächsten Schritt spiegeln wir die erste Spalte der $verbleibenden Untermatrix$ :

\tilde{a}_2 \text{ erste Spalte von } \tilde{A}_2

u_2 = \tilde{a}_2 \mp \|\tilde{a}_2\|_2 \cdot \tilde{e}_1

Der zweite Reflektor wird in den rechten unteren Block der Einheitsmatrix eingebettet:

\tilde{S}_2 = I_{m-1} - 2v_2v_2^T \in \mathbb{R}^{m-1 \times m-1}

Danach sind auch die Einträge der zweiten Spalte unterhalb der Diagonale null:

A_3 = S_2 \cdot A_2 = \begin{pmatrix} r_{11} & * & * \\ 0 & r_{22} & * \\ 0 & 0 & \tilde{A}_3 \end{pmatrix}

Dabei gilt $r_{22}=\pm\|\widetilde a_2\|_2$ und $\widetilde A_3\in\mathbb{R}^{(m-2)\times(n-2)}$ . Durch Wiederholung bleibt schließlich die letzte aktive Spalte

\tilde{A}_n = \begin{pmatrix} \tilde{a}_{nn} \\ \vdots \\ \tilde{a}_{mn} \end{pmatrix} \in \mathbb{R}^{m-(n-1) \times 1}

Nach der letzten Spiegelung ist die transformierte Matrix oberdreieckig:

\begin{aligned} A_{n+1} &= R \\ &= \begin{pmatrix} r_{11} & & * \\ & \ddots & \\ 0 & & r_{nn} \\ 0 & \cdots & 0 \\ \vdots & & \vdots \\ 0 & \cdots & 0 \end{pmatrix} \\ &= S_n \cdot \ldots \cdot S_1 \cdot A \\ &= Q^T \cdot A \end{aligned}

Da jeder Reflektor orthogonal und symmetrisch ist, liefert ihr Produkt

A = Q \cdot R

Die Matrix $Q$ ist als Produkt orthogonaler Reflektoren ebenfalls orthogonal.

Implementierungsalgorithmus

Zuerst werden die aktive Matrix und der gesammelte orthogonale Faktor initialisiert:

A_1 := A

Im Schritt $i$ wird der aktive Teil der $i$ -ten Spalte entnommen:

\begin{aligned} \tilde{a}_i &= \begin{pmatrix} a_{ii} \\ \vdots \\ a_{mi} \end{pmatrix} \\ &:= \begin{pmatrix} (A_i)_{ii} \\ \vdots \\ (A_i)_{mi} \end{pmatrix} \in \mathbb{R}^{m-i+1} \end{aligned}

Danach werden Norm, Vorzeichen und Householder-Vektor berechnet. Das Vorzeichen verhindert die $Subtraktion fast gleicher Zahlen$ :

\sigma := \|\tilde{a}_i\|_2

Der zugehörige Skalenfaktor ist

\beta := \frac{2}{u_i^Tu_i}=\frac{1}{\sigma(\sigma+|\tilde a_{ii}|)}

Der normierte Vektor und der vollständige Reflektor lauten

v_i=\frac{u_i}{\|u_i\|_2}

Die aktive Matrix und der gesammelte Faktor werden ohne ein zusätzliches dichtes Matrixprodukt aktualisiert:

\begin{aligned} A_{i+1} &:= S_i A_i \\ &= (I_m - \beta \hat{u}_i \hat{u}_i^T) A_i \\ &= A_i - \hat{u}_i (\beta \hat{u}_i^T A_i) \end{aligned}

Nach $n$ Schritten gilt

R := A_{n+1}

Numerische Eigenschaften und Aufwand

Die Umsetzung der Zerlegung in Gleitkommaarithmetik wirft zwei praktische Fragen auf, nämlich wie viele Operationen sie kostet und wie sich der $Rundungsfehler$ verhält. Beide Antworten entscheiden, ob das Verfahren zum Problem passt: Die Operationszahl bestimmt die Laufzeit bei einer großen Matrix, und das Rundungsverhalten bestimmt, wie viel vom berechneten Ergebnis du verwerten kannst.

Rechenkomplexität

Die folgende Zählung speichert $Q$ $implizit als Reflektoren$ und berechnet $R$ . Das explizite Bilden der vollständigen $m\times m$ -Matrix $Q$ benötigt zusätzliche Arbeit.

Für $m\ge n$ benötigt eine $\mathrm{QR}$ -Zerlegung mit Householder-Spiegelungen ungefähr:

2mn^2-\frac23n^3

Die ungefähre Anzahl an $Gleitkommaoperationen$ für $m=n$ beträgt:

\frac43n^3

Numerische Eigenschaften

Orthogonale Multiplikationen $erhalten die Singulärwerte$ . Daher gilt für eine Matrix mit vollem Spaltenrang $\operatorname{cond}_2(R)=\operatorname{cond}_2(A)$ .

Im Schritt $i$ entsteht der neue Diagonaleintrag $\pm\|\widetilde a_i\|_2$ . Das folgende Vorzeichen verhindert beim Berechnen des Reflektors die Subtraktion fast gleicher Zahlen:

u_i = \tilde{a}_i + s_i\|\tilde{a}_i\|_2e_1

Nach der Faktorisierung können die Vorzeichen einer Spalte von $Q$ und der zugehörigen Zeile von $R$ gemeinsam gewechselt werden, damit die Diagonale von $R$ positiv ist. Gilt $\|\widetilde a_i\|_2=0$ , fehlt der volle Spaltenrang und diese Eindeutigkeitsbedingung entfällt. Eine $\mathrm{QR}$ -Zerlegung mit Spaltenpivotisierung eignet sich besser, wenn der numerische Rang bestimmt werden soll.

Eine kompakte Implementierung skaliert den Householder-Vektor so, dass sein erster aktiver Eintrag 1 ist. Die übrigen Einträge werden unterhalb der Diagonalen gespeichert, die Reflektorstärke in einem separaten Skalar. Diagonale und oberes Dreieck enthalten weiterhin $R$ . Die gespeicherten Reflektoren lassen sich ohne explizites $Q$ anwenden.

Die reduzierte Faktorisierung behält die ersten $n$ Spalten von $Q$ und die ersten $n$ Zeilen von $R$ .