Householder-Transformation schrittweise

import numpy as np
from numpy.linalg import norm

Die allgemeine Schreibweise der Householder-Transformation für einen beliebigen Vektor \(w\) ist gegeben durch \begin{equation}\label{eq:householdertransformation}H(w) = \text{id} - 2\,\frac{w\cdot w^T}{\langle w, w\rangle} \end{equation} wobei \(w\cdot w^T\) das Kroneckerprodukt

\[w\cdot w^T = (w_i\,w_j)_{i,j=1\ldots m}\in\mathbb{R}^{m\times m}\]

sei.

# Selber implementieren
def HouseholderTransformation(w):
    return <<snipp>>

Gesucht ist der geeignete Normalenvektor so, dass der gespiegelte Spaltenvektor auf die \(e_1\) Achse zu liegen kommt. Sei mit \(y\) der Spaltenvektor bezeichnet, so kann man zeigen (siehe Skript), dass der Vektor \begin{equation} w = y \pm \|y\|_2 e_1 \end{equation} die gewünschte Eigenschaft hat. Um Auslöschung in der Berechnung von \(w\) zu vermeiden, wählt man \begin{equation} w = y + \text{sign}(y_1) \|y\|_2 e_1 \end{equation} mit \begin{equation} \text{sign}(s) = \begin{cases} 1 & \quad \text{für} s \ge 0\\ -1 & \quad \text{sonst}.\end{cases} \end{equation}

def mysign(x): # numpy sign liefert 0 für 0
    if x >= 0:
        return 1
    else:
        return -1

Funktion für den n-dimensionalen Einheitsvektor

def e(n):
    return np.array([1]+[0 for k in range(n-1)])

Mit Hilfe der Householder-Transformation soll nun die Matrix \(A\) in eine orthogonale Matrix \(Q\) und reguläre obere Dreiecksmatrix \(R\) zerlegt werden. Im Beispiel wählen wir eine zufällige Matrix \(A \in \mathbb{R}^{10\times5}\).

A = np.array([[-1,  7, -8, -9,  6],
       [-6, -8,  0,  3,  8],
       [-4, -2,  8,  0, -2],
       [-1, -9,  4, -8,  2],
       [-3, -5, -5,  7, -4],
       [-7, -4,  7, -1,  5],
       [-9, -7,  6, -5, -8],
       [-4, -3, -5,  3, -6],
       [ 5,  7,  5, -4, -5],
       [ 4, -6, -8, -2, -5]],dtype=float)
m,n = A.shape

1. Spalte

k = 0

Die Hyperebene ist definiert durch

w = <<snipp>>

Für die Householder-Transformationsmatrix angewand auf \(A\) erhalten wir

Q1 = HouseholderTransformation(w)
A1 = Q1@A

In der ersten Spalte der Zwischenmatrix \(A_1\) stehen nun abgesehen vom ersten Eintrag Nullen:

print(np.round(A1,4))

2. Spalte

k = 1

Die Hyperebene ist definiert durch

w = <<snipp>>

wobei nun das letzte Resultat \(A_1\) benutzt wird. Die Householder-Transformationsmatrix wird nun nur auf die Submatrix von \(A_1\) angewand und in der Submatrix von \(A_1\) wiederum gespeichert:

Q2 = HouseholderTransformation(w)
A1[k:,k:] = Q2@A1[k:,k:]

Die Dimension der zweiten Householder-Transformationsmatrix \(Q_2\) ist

Q2.shape

In dem ersten beiden Spalte der Zwischenmatrix \(A_1\) steht:

print(np.round(A1,4))

3. - 5. Spalte

Wir automatisieren nun den Prozess und überschreiben die Submatrizen der Matrix \(A_1\) sukzessive:

for k in range(2,n):
    print('Spalte '+str(k+1))
    w = <<snipp>>
    Qk = HouseholderTransformation(w)
    A1[k:,k:] = <<snipp>>
    print(np.round(A1,4))

Q, R berechnen

• Berechnen sie abschliessend \(Q,R\) so, dass \(Q\cdot R = A\) gilt.

• Vergleichen Sie Ihr Resultat mit der Funktion von NumPy.

Q,R = np.linalg.qr(A)

Q.shape

R.shape

Frage: Warum reicht diese reduzierte \(Q\) und \(R\) Matrix?