3 Matrix-Zerlegungen

Die Lösung \(w\) des Problems der linearen Ausgleichsrechnung war entweder als Lösung eines Optimierungsproblems \[\begin{equation*} \min_{w} \| Aw - y \|^2 \end{equation*}\] oder als Lösung des linearen Gleichungssystems \[\begin{equation*} A^TAw=y \end{equation*}\] gegeben. Hierbei steht nun \(A\in \mathbb R^{N\times n}\) für die Matrix \(\Phi(\mathbf x)\) der Daten und Basisfunktionen. Wir hatten uns überlegt, dass in den meisten Fällen

  • die Matrix mehr Zeilen als Spalten hat (\(N>n\)) und
  • die Spalten linear unabhängig sind.

3.1 QR Zerlegung

Wir betrachten nochmal das Optimierungsproblem \(\min_{w} \| Aw - y \|^2\). Gäbe es eine Lösung des Systems \(Aw=y\), wäre das sofort eine Lösung des Optimierungsproblems. Da \(A\) aber mehr Zeilen als Spalten hat, ist das System \(Aw=y\) überbestimmt und eine Lösung in der Regel nicht gegeben.

Die Überlegung ist nun, die Gleichung \(Aw=y\) so gut wie möglich zu erfüllen, indem wir die relevanten Gleichungen indentifizieren und wenigstens diese lösen. Ein systematischer (und wie wir später sehen werden auch zum Optimierungsproblem passender) Zugang bietet die QR Zerlegung.

Theorem 3.1 (QR Zerlegung) Sei \(A\in \mathbb R^{m\times n}\), \(m>n\). Dann existiert eine orthonormale Matrix \(Q\in \mathbb R^{m\times m}\) und eine obere Dreiecksmatrix \(\hat R\in \mathbb R^{n\times n}\) derart dass \[\begin{equation*} A = QR =: Q \begin{bmatrix} \hat R \\ 0 \end{bmatrix}. \end{equation*}\] Hat \(A\) vollen (Spalten)Rang, dann ist \(\hat R\) invertierbar.

Hier heißt orthonormale Matrix \(Q\), dass die Spalten von \(Q\) paarweise orthogonal sind. Insbesondere gilt \[\begin{equation*} Q^T Q = I. \end{equation*}\]

Für unser zu lösendes Problem ergibt sich dadurch die Umformung \[\begin{equation*} Aw = y \quad \Leftrightarrow \quad QRw=y \quad \Leftrightarrow \quad Q^TQRw=Q^T y \quad \Leftrightarrow \quad Rw = Q^Ty \end{equation*}\] oder auch \[\begin{equation*} \begin{bmatrix} \hat R \\ 0 \end{bmatrix}w = Q^Ty = \begin{bmatrix} Q_1^T \\ Q_2^T \end{bmatrix} y \end{equation*}\] (wobei wir die \(Q_1\in \mathbb R^{m\times n}\) die Matrix der ersten \(n\) Spalten von \(Q\) ist) und als Kompromiss der Vorschlag, das Teilsystem \[\begin{equation*} \hat R w = Q_1^Ty \end{equation*}\] nach \(w\) zu lösen und in Kauf zu nehmen, dass der Rest, nämlich das \(Q_2^Ty\), nicht notwendigerweise gleich null ist.

Wir halten zunächst mal fest, dass

  • Obwohl \(Q\) eine reguläre Matrix ist, bedarf der Übergang von \(Aw=y\) zu \(Q^TAw=Q^Ty\) einer genaueren Analyse.
  • Wir bemerken, dass für eine hypothetische komplette Lösung \(Aw-y\), diese Transformation keine Rolle spielt.
  • Für die Kompromisslösung jedoch schon, weil beispielsweise verschiedene Konstruktionen eines invertierbaren Teils, verschiedene Residuen bedeuten und somit Optimalität im Sinne von \(\min_w \|Aw-y\|^2\) nicht garantiert ist.

Allerdings, wie Sie als Übungsaufgabe nachweisen werden, löst dieser Ansatz tatsächlich das Optimierungsproblem.

3.2 Singulärwertzerlegung

Eine weitere Matrix Zerlegung, die eng mit der Lösung von Optimierungsproblemen oder überbestimmten Gleichungssystemen zusammenhängt ist die Singulärwertzerlegung (SVD – von english: Singular Value Decomposition).

Theorem 3.2 (Singulärwertzerlegung) Sei \(A\in \mathbb C^{m\times n}\), \(m\geq n\). Dann existieren orthogonale Matrizen \(U \in \mathbb C^{m\times m}\) und \(V\in \mathbb C^{n\times n}\) und eine Matrix \(\Sigma \in \mathbb R^{m\times n}\) der Form \[\begin{equation*} \Sigma = \begin{bmatrix} \sigma_1 & 0 & \dots & 0\\ 0 & \sigma_2 &\ddots & \vdots\\ 0 & \ddots & \ddots &0\\ 0 & \dots&0 & \sigma_n \\ 0 & 0 & \dots & 0 \\ \vdots & \ddots & & \vdots\\ 0 & 0 & \dots & 0 \end{bmatrix} \end{equation*}\] mit reellen sogenannten Singulärwerten \[\begin{equation*} \sigma_1 \geq \sigma_2 \geq \dots \geq \sigma_n \geq 0 \end{equation*}\] sodass gilt \[\begin{equation*} A = U \Sigma V^* \end{equation*}\] wobei gilt \(V^* = \overline{V^T}\) (transponiert und komplex konjugiert).

Ein paar Bemerkungen.

  • Ist \(A\) reell, können auch \(U\) und \(V\) reell gewählt werden.
  • Die Annahme \(m \geq n\) war nur nötig um für die Matrix \(\Sigma\) keine Fallunterscheidung zu machen. (Für \(m\leq n\) “steht der Nullblock rechts von den Singulärwerten”). Insbesondere gilt \(A^* = V\Sigma U^*\) ist eine SVD von \(A^*\).
  • Eine Illustration der Zerlegung ist hier zu sehen.

Wir machen einige Überlegungen im Hinblick auf große Matrizen. Sei dazu \(m>n\), \(A\in \mathbb C^{m\times n}\) und \(A=U\Sigma V^*\) eine SVD wie in Theorem 3.2. Sei nun \[\begin{equation*} U = \begin{bmatrix} U_1 & U_2 \end{bmatrix} % = \begin{bmatrix} V_1^* & V_2^* \end{equation*}\] partitioniert sodass \(U_1\) die ersten \(n\) Spalten von \(U\) enthält.

Dann gilt (nach der Matrix-Multiplikations Regel Zeile mal Spalte die Teile \(U_2\) und \(V_2\) immer mit dem Nullblock in \(\Sigma\) multipliziert werden) dass \[\begin{equation*} A = U\Sigma V = \begin{bmatrix} U_1 & U_2 \end{bmatrix} \begin{bmatrix} \hat \Sigma \\ 0 \end{bmatrix} V^* % \begin{bmatrix} V_1^* \\ V_2^* \end{bmatrix} = U_1 \hat \Sigma V^* % \begin{bmatrix} V_1^* \\ V_2^* \end{bmatrix} \end{equation*}\] Es genügt also nur die ersten \(m\) Spalten von \(U\) zu berechnen. Das ist die sogenannte slim SVD.

Hat, darüberhinaus, die Matrix \(A\) keinen vollen Rang, also \(\operatorname{Rg}(A) = r < n\), dann:

  • ist \(\sigma_i=0\), für alle \(i=r+1, \dotsc, n\), (wir erinnern uns, dass die Singulärwerte nach Größe sortiert sind)
  • die Matrix \(\hat \Sigma\) hat \(n-r\) Nullzeilen
  • für die Zerlegung sind nur die ersten \(r\) Spalten von \(U\) und \(V\) relevant – die sogenannte Kompakte SVD.

In der Datenapproximation ist außerdem die truncated SVD von Interesse. Dazu sei \(\hat r<r\) ein beliebig gewählter Index. Dann werden alle Singulärwerte, \(\sigma_i=0\), für alle \(i=\hat r+1, \dotsc, n\), abgeschnitten – das heißt null gesetzt und die entsprechende kompakte SVD berechnet.

Hier gilt nun nicht mehr die Gleichheit in der Zerlegung. Vielmehr gilt \[\begin{equation*} A \approx A_{\hat r} \end{equation*}\] wobei \(A_{\hat r}\) aus der truncated SVD von \(A\) erzeugt wurde. Allerdings ist diese Approximation von \(A\) durch optimal in dem Sinne, dass es keine Matrix vom Rang \(\hat r \geq r=\operatorname{Rg}(A)\) gibt, die \(A\) (in der induzierten euklidischen Norm1) besser approximiert. Es gilt \[\begin{equation*} \min_{B\in \mathbb C^{m\times n}, \operatorname{Rg}(B)=\hat r} \|A-B\|_2 = \|A-A_{\hat r}\|_2 = \sigma_{\hat r + 1}; \end{equation*}\] (vgl. Satz 14.15, Bollhöfer and Mehrmann 2004).

Zum Abschluss noch der Zusammenhang zum Optimierungsproblem. Ist \(A=U\Sigma V^*\) “SV-zerlegt”, dann gilt \[\begin{equation*} A^*Aw = V\Sigma^*U^*U\Sigma V^*w = V\hat \Sigma^2 V^* \end{equation*}\] und damit \[\begin{equation*} A^*Aw = A^*y \quad \Leftrightarrow \quad V\hat \Sigma^2 V^* = V\Sigma^*U^*y \quad \Leftrightarrow \quad w = V(\Sigma^+)^*U^*y \end{equation*}\] wobei \[\begin{equation*} \Sigma^+ = \begin{bmatrix} \hat \Sigma^{-1} \\ 0_{m-n \times n} \end{bmatrix} \end{equation*}\] aus \(\Sigma \hat \Sigma^{-1}\hat \Sigma^{-1}\) herrührt.

Bemerkung: \(\Sigma^+\) kann auch definiert werden, wenn \(\hat \Sigma\) nicht invertierbar ist (weil manche Diagonaleinträge null sind). Dann wird \(\hat \Sigma^+\) betrachtet, bei welcher nur die \(\sigma_i>0\) invertiert werden und die anderen \(\sigma_i=0\) belassen werden. Das definiert eine sogenannte verallgemeinerte Inverse und löst auch das Optimierungsproblem falls \(A\) keinen vollen Rang hat.

Illustration der SVD. Bitte beachten der \(*\) bedeutet hier transponiert und komplex konjugiert. By Cmglee - Own work, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=67853297

3.3 Aufgaben

Erklärung: (T) heißt theoretische Aufgabe, (P) heißt programmieren.

3.3.1 Norm und Orthogonale Transformation (T)

Sei \(Q\in \mathbb R^{n\times n}\) eine orthogonale Matrix und sei \(y\in \mathbb R^{n}\). Zeigen Sie, dass \[\begin{equation*} \|y\|^2 = \|Qy \|^2 \end{equation*}\] gilt. (2 Punkte)

3.3.2 Kleinste Quadrate und Mittelwert

Zeigen sie, dass der kleinste Quadrate Ansatz zur Approximation einer Datenwolke \[\begin{equation*} (x_i, y_i), \quad i=1,2,\dotsc,N, \end{equation*}\] mittels einer konstanten Funktion \(f(x)=w_1\) auf \(w_1\) auf den Mittelwert der \(y_i\) führt. (6 Punkte)

3.3.3 QR Zerlegung und Kleinstes Quadrate Problem (T)

Sei \(A\in \mathbb R^{m,n}\), \(m>n\), \(A\) hat vollen Rank und sei \[\begin{equation*} \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} \hat R \\ 0 \end{bmatrix} = A \end{equation*}\] eine QR-Zerlegung von \(A\). Zeigen sie, dass die Lösung von \[\begin{equation*} \hat R w = Q_1^T y \end{equation*}\] ein kritischer Punkt (d.h. der Gradient \(\nabla_w\) verschwindet) von \[\begin{equation*} w \mapsto \frac 12 \| Aw - y \|^2 \end{equation*}\] ist. Hinweis: Die Formel für den Gradienten wurde in der Vorlesung 02 hergeleitet. (6 Punkte)

3.3.4 Eigenwerte Symmetrischer Matrizen (T)

Zeigen Sie, dass Eigenwerte symmetrischer reeller Matrizen \(A\in \mathbb R^{n\times n}\) immer reell sind. (3 Punkte)

3.3.5 Singulärwertzerlegung und Eigenwerte (T)

Zeigen Sie, dass die quadrierten Singulärwerte einer Matrix \(A\in \mathbb R^{m\times n}\), \(m>n\), genau die Eigenwerte der Matrix \(A^TA\) sind und in welcher Beziehung sie mit den Eigenwerten von \(AA^T\) stehen. Hinweis: hier ist “\(m>n\)” wichtig. (6 Punkte)

3.3.6 Python – Laden und Speichern von Arrays

Speichern sie die Covid-Daten aus obiger Tabelle zur späteren Verwendung als ein numpy.array. Beispielsweise so:

import numpy as np
import matplotlib.pyplot as plt

days = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
        12, 13, 14, 15,  16, 17, 18, 19,  20,   21,
        22,  23,  24,   25,    26, 27,  28,   29,    30,     31]
case = [352, 347, 308, 151, 360, 498, 664, 686, 740, 418, 320,
        681, 691, 1154, 1284,  127, 984, 573, 1078,  1462,   2239,
        2236,  2119,  1663,   1413,    2283,
        2717,  3113,   2972,    3136,    2615]

data = np.vstack([days, case])
print('Shape of data: ', data.shape)

datafilestr = 'coviddata.npy'
np.save(datafilestr, data)

lddata = np.load(datafilestr)

plt.figure(1)
plt.plot(lddata[0, :], lddata[1, :], 's', label='Cases/Day')
plt.title('Covid Faelle in Bayern im Oktober 2020')
plt.legend()
plt.show()

3.3.7 Lineare Regression für Covid Daten (P)

Führen sie auf den Covid-daten eine lineare Regression zum Fitten

  • einer konstanten Funktion \(f(x)=c\)
  • einer linearen Funktion \(f(x) = ax+b\)
  • einer quadratischen Funktion \(f(x) = w_1 + w_2 x + w_3x^2\)

durch. Berechnen Sie die mittlere quadratische Abweichung \(\frac{1}{N}\sum_{i=1}^N\|y_i - f(x_i)\|^2\) für alle drei Approximationen und plotten Sie die Covid-Daten zusammen mit der aus der Regression erhaltenen Funktion. Beispielsweise so:

import numpy as np
import matplotlib.pyplot as plt

datafilestr = 'coviddata.npy'
cvdata = np.load(datafilestr)

# ## Definieren der Basis Funktion(en)


def b_zero(x):
    ''' eine konstante Funktion '''
    return 1.


ybf = cvdata[1, :]  # die y-werte
xbf = cvdata[0, :]  # die x-werte

N = xbf.size   # Anzahl Datenpunkte
ybf = ybf.reshape((N, 1))  # reshape als Spaltenvektor

bzx = [b_zero(x) for x in xbf]  # eine Liste mit Funktionswerten
bzx = np.array(bzx).reshape((N, 1))  # ein Spalten Vektor

Phix = bzx  # hier nur eine Spalte

# Das LGS AtA w = At y
rhs = Phix.T @ ybf
AtA = Phix.T @ Phix

# ACHTUNG: das hier geht nur weil AtA keine Matrix ist in diesem Fall
w = 1./AtA * rhs
# ACHTUNG: das hier ging nur weil AtA keine Matrix ist in diesem Fall


def get_regfunc(weights, basfunlist=[]):
    ''' Eine Funktion, die eine Funktion erzeugt

    Eingang: die Gewichte, eine Liste von Basisfunktionen
    Ausgang: die entsprechende Funktion `f = w_1b_1 + ... `
    '''
    def regfunc(x):
        fval = 0
        for kkk, basfun in enumerate(basfunlist):
            fval = fval + weights[kkk]*basfun(x)
        return fval

    return regfunc


const_regfunc = get_regfunc([w], [b_zero])

const_approx = [const_regfunc(x) for x in xbf]
const_approx = np.array(const_approx).reshape((N, 1))

plt.figure(1)
plt.plot(cvdata[0, :], cvdata[1, :], 's', label='Cases/Day')
plt.plot(cvdata[0, :], const_approx, '-', label='constant fit')
plt.legend()
plt.show()

Hinweise:

  • liste = [f(x) for x in iterable] ist sehr bequem um Vektoren von Funktionswerten zu erzeugen aber nicht sehr pythonesque (und auch im Zweifel nicht effizient). Besser ist es Funktionen zu schreiben, die vektorisiert sind. Zum Beispiel können die meisten built-in Funktionen wie np.exp ein array als Eingang direkt in ein array der Funktionswerte umsetzen.
  • Eine Funktion, die eine Funktion erzeugt finde ich sehr hilfreich für viele Anwendungen (ist aber manchmal nicht so gut nachvollziehbar).

3.3.8 Truncated SVD (P+T)

  1. (P) Berechnen und plotten sie die Singulärwerte einer \(4000\times 1000\) Matrix mit zufälligen Einträgen und die einer Matrix mit “echten” Daten (hier Simulationsdaten einer Stroemungssimulation)2. Berechnen sie den Fehler der truncated SVD \(\|A-A_{\hat r}\|\) für \(\hat r = 10, 20, 40\) für beide Matrizen.
  2. (T) Was lässt sich bezüglich einer Kompression der Daten mittels SVD für die beiden Matrizen sagen. (Vergleichen sie die plots der Singulärwerte und beziehen sie sich auf die gegebene Formel für die Differenz).
  3. (P+T) Für die “echten” Daten: Speichern sie die Faktoren der bei \(\hat r=40\) abgeschnittenen SVD und vergleichen Sie den Speicherbedarf der Faktoren und der eigentlichen Matrix.

Beispielcode:

import numpy as np
import scipy.linalg as spla
import matplotlib.pyplot as plt

randmat = np.random.randn(4000, 1000)

rndU, rndS, rndV = spla.svd(randmat)

print('U-dims: ', rndU.shape)
print('V-dims: ', rndV.shape)
print('S-dims: ', rndS.shape)

plt.figure(1)
plt.semilogy(rndS, '.', label='Singulaerwerte (random Matrix)')

realdatamat = np.load('velfielddata.npy')

# # Das hier ist eine aufwaendige Operation
rlU, rlS, rlV = spla.svd(realdatamat, full_matrices=False)
# # auf keinen Fall `full_matrices=False` vergessen

print('U-dims: ', rlU.shape)
print('V-dims: ', rlV.shape)
print('S-dims: ', rlS.shape)

plt.figure(1)
plt.semilogy(rlS, '.', label='Singulaerwerte (Daten Matrix)')

plt.legend()
plt.show()

Hinweise:

  • Es gibt viele verschiedene Normen für Vektoren und Matrizen. Sie dürfen einfach mit np.linalg.norm arbeiten. Gerne aber mal in die Dokumentation schauen welche Norm berechnet wird.
  • Die (T) Abschnitte hier bitte mit den anderen (T) Aufgaben oder als Bildschirmausgabe im Programm.

Referenzen

Bollhöfer, M., Mehrmann, V.: Numerische Mathematik. Eine projektorientierte Einführung für Ingenieure, Mathematiker und Naturwissenschaftler. Vieweg (2004)

  1. Auf Matrixnormen kommen wir noch in der Vorlesung zu sprechen.↩︎

  2. Download bitte hier – Achtung das sind 370MB↩︎