Householdertransformation

(Weitergeleitet von Householder-Spiegelung)

In der Mathematik beschreibt die Householdertransformation die Spiegelung eines Vektors an einer Hyperebene durch Null im euklidischen Raum. Im dreidimensionalen Raum ist sie somit eine Spiegelung an einer Ebene (durch den Ursprung). Die Darstellung dieser linearen Abbildung durch eine Matrix wird als Householder-Matrix bezeichnet. Verwendung findet sie vor allem in der numerischen Mathematik, wenn mittels orthogonaler Transformationen Matrizen so gezielt umgeformt werden, dass bestimmte Spaltenvektoren auf das Vielfache des ersten Einheitsvektors abgebildet werden, insbesondere beim QR-Verfahren und der QR-Zerlegung.

Die Householdertransformation wurde 1958 durch den amerikanischen Mathematiker Alston Scott Householder eingeführt.

Definition und Eigenschaften

Bearbeiten
 
Illustration der Householder-Transformation in 2D

Die Spiegel-Hyperebene kann durch einen Normalenvektor  , also einen Vektor, der orthogonal zur Hyperebene ist, definiert werden. Ist   als Spaltenvektor gegeben und   die Einheitsmatrix, dann wird die oben beschriebene lineare Abbildung durch die folgende Matrix dargestellt:

 

Dabei bezeichnet   die Transponierte des Spaltenvektors  , also einen Zeilenvektor. Der Nenner   ist das Skalarprodukt von   mit sich selbst,   das dyadische Produkt. Die Matrix   beschreibt die Orthogonalprojektion auf die durch   gegebene Richtung. Ist   auf die Länge eins normiert, also  , so vereinfacht sich die Formel zu

 

Die Spiegelungseigenschaft ersieht man daraus, dass

 ,

wobei   das Standardskalarprodukt bezeichnet. Der Term   entspricht dabei dem Abstand des Punktes   zur Hyperebene  . Der Vektor   wird also in zwei zueinander orthogonale Anteile zerlegt, wobei der erste Anteil in der Hyperebene liegt und der zweite ein Vielfaches des Vektors   ist. Unter der Spiegelung wird der Anteil in der Ebene invariant gelassen, der Anteil in Richtung  , also senkrecht zur Ebene, wird „umgeklappt“, also nun abgezogen statt addiert.

Die Householder-Matrix hat folgende Eigenschaften:

  • Sie ist symmetrisch:  
  • Sie ist orthogonal:  
  • Sie ist involutorisch:   (Dies folgt aus der Symmetrie und der Orthogonalität.)
  • Sie hat den einfachen Eigenwert −1 zum Eigenvektor   und den  -fachen Eigenwert 1. Der Eigenraum zum Eigenwert 1 ist die Spiegelebene, also das orthogonale Komplement des von   erzeugten eindimensionalen Unterraums.
  • Matrix-Vektor-Multiplikationen mit   sind schnell berechenbar.

Konstruktion einer spezifischen Spiegelung

Bearbeiten

Es sei ein Vektor   gegeben, der auf ein Vielfaches des Vektors   gespiegelt werden soll, das heißt, gesucht ist ein Einheitsvektor  , so dass mit der zugehörigen Householder-Matrix   gilt  . Geometrisch ist der Vektor   die Richtung einer der zwei Winkelhalbierenden der Geraden in Richtung   und in Richtung  . Die Winkelhalbierende ergibt sich, indem man auf beiden Geraden Punkte mit demselben Abstand zum Nullpunkt wählt und auf der Verbindungsstrecke dieser zwei Punkte den Mittelpunkt konstruiert. Die Gerade durch Nullpunkt und Mittelpunkt hat dann die gesuchte Richtung  , der Vektor   selbst ergibt sich durch Normieren dieser Richtung. Die zweite Winkelhalbierende ergibt sich, indem die Konstruktion ausgehend von   und   durchgeführt wird.

Der Einfachheit halber sei   normiert,  . Dann muss, wegen der Orthogonalität der Spiegelung,   gelten. Der gesuchte Spiegelungsvektor   ergibt sich nun durch Normieren des Differenzvektors  , also

 .

Beide Vorzeichenvarianten führen zum gewünschten Ergebnis (sofern der Nenner von Null verschieden ist). Aus Gründen numerischer Stabilität wird das Vorzeichen von   so gewählt, dass der Nenner am größten ist, also   gilt.

In der Probe ergibt sich

 

Beispiel

Bearbeiten

Am häufigsten wird der Fall betrachtet, in dem   der erste kanonische Basisvektor ist. Sei   in erste Komponente und Restvektor zerlegt. Dann gilt für die Norm  . Als Vorzeichen von   ist das Vorzeichen von   zu wählen, die Richtung der Spiegelung ist dann

 .

Dabei ist  

Der Vektor   entsteht durch Normierung dieser Richtung. Nach Umformen stellt sich die Norm der Richtung als

 

dar, wobei in dieser Form nur bereits berechnete Zwischenergebnisse benutzt werden. In der unnormierten Variante der Spiegelung ergeben sich weitere Einsparungen an Rechenschritten.

Anwendung: QR-Zerlegung

Bearbeiten

Householder-Spiegelungen können zur stabilen Berechnung von QR-Zerlegungen einer Matrix   verwendet werden, indem zunächst die erste Spalte der Matrix mit einer Spiegelung   auf das Vielfache des ersten Einheitsvektors gespiegelt wird, wie im letzten Abschnitt erläutert (jetzt bezeichnet der Index aber die Nummer der Spiegelung).

Danach behandelt man   mit einer Spiegelung   analog, wobei die Spiegelung so konstruiert wird, dass erste Zeile und Spalte von der Transformation unberührt bleiben. Dies wird erreicht, indem die erste Komponente des Spiegelungsvektors zu Null gesetzt wird. Zur Bestimmung des dritten Schrittes geht analog nur die Hauptuntermatrix unter dem dritten Diagonalelement ein, der Spiegelungsvektor ist Null in den ersten zwei Komponenten etc. Im  -ten Schritt wird also die Untermatrix unter der Position   des Produkts   auf die gleiche Art reduziert, bis die Restmatrix   Dreiecksgestalt besitzt.

Mit   gilt  , also ergibt sich die QR-Zerlegung

 

mit

 

Man beachte, dass   hier eine quadratische Matrix ist. Meist werden die Matrizen   bzw.   nicht explizit berechnet, sondern man nutzt direkt die Produktform. Dazu werden die Spiegelvektoren   von   im frei gewordenen Platz der Matrix   gespeichert.

Die Zahl der Operationen für die QR-Zerlegung einer Matrix   mit dem Householder-Verfahren beträgt:

 

Im Fall   genügen   Schritte, da die letzte Spalte nicht mehr transformiert werden muss.

Pseudocode

Bearbeiten

Da für die meisten Berechnungen das explizite Ausrechnen von   nicht nötig ist, reicht es, nur die Matrix   zu berechnen.   ist die linke Spalte der jeweiligen Untermatrix. Bei der unten angegebenen Funktion wird das Ergebnis direkt in   geschrieben, so dass nach Abarbeitung des Algorithmus das   in   steht. Die Zeile   könnte also auch weggelassen werden.

  function GetR(A)
     for k=1…n
         z=A(k…m,k)
         uk=z
         uk(1)+=sign(z(1))*norm(z)
         uk=uk/norm(uk)
         vk=zeros(m)
         vk(k…m)= uk
         A=A-(2*vk)*(vk'*A)
     R=A
     return R

Sollte   dennoch benötigt werden, lässt sich das obere Beispiel einfach erweitern:

  function GetR(A)
    Q=eye(m) 
    for k=1…n
         z=A(k…m,k)
         uk=z
         uk(1)+=sign(z(1))*norm(z)
         uk=uk/norm(uk)
         vk=zeros(m)
         vk(k…m)= uk
         A=A-(2*vk)*(vk'*A)
         Q=Q-Q*vk*(2*vk')
     R=A
     return R

Programmierung

Bearbeiten

Das folgende Beispiel in der Programmiersprache C++ zeigt die Implementierung der Householdertransformation. Bei der Ausführung des Programms wird die Funktion main verwendet, die das Ergebnis auf der Konsole ausgibt.[1]

#include <iostream>
#include <vector>
using namespace std;

// Diese Funktion berechnet c = a + b * s
void vmadd(const Vector& a, const Vector& b, double s, Vector& c)
{
    if (c.size != a.size || c.size != b.size)
    {
        return;
    }
    for (int i = 0; i < c.size; i++)
    {
        c(i) = a(i) + s * b(i);
    }
}

// Diese Funktion berechnet matrix = I - 2 * v * v^T
void computeHouseholderFactor(Matrix& matrix, const Vector& v)
{
    int n = v.size;
    matrix.allocate(n, n);
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            matrix(i, j) = -2 * v(i) * v(j);
        }
    }
    for (int i = 0; i < n; i++)
    {
        matrix(i, i) += 1;
    }
}

void Matrix::extractColumn(Vector& vector, int c)
{
    if (m != vector.size)
    {
        cerr << "[Matrix::extract_column]: Matrix and Vector sizes don't match\n";
        return;
    }
    for (int i = 0; i < m; i++)
    {
        vector(i) = (*this)(i, c);
    }
}

// Diese Funktion gibt die Matrix auf der Konsole aus
void matrixShow(const Matrix& matrix, const string& str = "")
{
    cout << str << "\n";
    for (int i = 0; i < matrix.m; i++)
    {
        for (int j = 0; j < matrix.n; j++)
        {
            printf(" %8.3f", matrix(i, j));
        }
        printf("\n");
    }
    printf("\n");
}

// Diese Funktion berechnet die L2-norm ||A - B||^2
double matrixCompare(const Matrix& A, const Matrix& B)
{
    if (A.m != B.m or A.n != B.n)
    {
        return numeric_limits<double>::max();
    }
    double norm = 0;
    for (int i = 0; i < A.m; i++)
    {
        for (int j = 0; j < A.n; j++)
        {
            norm += (A(i, j) - B(i, j)) * (A(i, j) - B(i, j));
        }
    }
    norm /= A.m * A.n;
    return norm;
}

// Diese Funktion bestimmt die QR-Zerlegung der gegebenen Matrix
void householder(Matrix& matrix, Matrix& R, Matrix& Q)
{
    int m = matrix.m;
    int n = matrix.n;
    vector<Matrix> qArray(m);
    Matrix z(matrix);
    Matrix z1;
    for (int k = 0; k < n && k < m - 1; k++)
    {
        Vector e(m), x(m);
        double a;
        z1.computeMinor(z, k);
        z1.extractColumn(x, k);
        a = x.calculateNorm();
        if (matrix(k, k) > 0)
        {
            a = -a;
        }
        for (int i = 0; i < e.size; i++)
        {
            e(i) = i == k;
        }
        vmadd(x, e, a, e);
        e.rescaleUnit();
        computeHouseholderFactor(qArray[k], e);
        z.multiply(qArray[k], z1);
    }
    Q = qArray[0];
    for (int i = 1; i < n && i < m - 1; i++)
    {
        z1.multiply(qArray[i], Q);
        Q = z1;
    }
    R.multiply(Q, matrix);
    Q.transpose();
}

// Hauptfunktion die das Programm ausführt
int main()
{
    double in[][3] = {
      { 12, -51,   4},
      {  6, 167, -68},
      { -4,  24, -41},
      { -1,   1,   0},
      {  2,   0,   3},
    };
    Matrix A(in);
    Matrix Q, R;
    matrixShow(A, "A"); // Aufruf der Funktion, gibt die Matrix A auf der Konsole aus
    // Berechnet die QR-Zerlegung
    householder(A, R, Q);
    matrixShow(Q, "Q"); // Aufruf der Funktion, gibt die Matrix Q auf der Konsole aus
    matrixShow(R, "R"); // Aufruf der Funktion, gibt die Matrix R auf der Konsole aus
    Matrix ACheck;
    ACheck.multiply(Q, R); // Berechnet das Matrixprodukt Q * R
    // Berechnet die L2-norm ||A - ACheck||^2
    double l2norm = matrixCompare(A, ACheck); // Aufruf der Funktion, vergleicht das Matrixprodukt Q * R mit der ursprünglichen Matrix A
    if (l2norm < 1e-12)
    {
        cout << "Die QR-Zerlegung ist korrekt." << endl; // Ausgabe auf der Konsole
    }
    else
    {
        cout << "Die QR-Zerlegung ist nicht korrekt." << endl; // Ausgabe auf der Konsole
    }
}

Siehe auch

Bearbeiten

Literatur

Bearbeiten
  • Gene H. Golub, Charles F. van Loan: Matrix Computations. 2nd Edition. The Johns Hopkins University Press, 1989.
  • Gerhard Opfer: Numerische Mathematik für Anfänger. Eine Einführung für Mathematiker, Ingenieure und Informatiker, 5. Aufl., Vieweg + Teubner, Wiesbaden 2008, ISBN 978-3-8348-0413-6
  • Martin Hermann: Numerische Mathematik, Band 1: Algebraische Probleme. 4., überarbeitete und erweiterte Auflage, Walter de Gruyter Verlag, Berlin und Boston 2020, ISBN 978-3-11-065665-7.

Einzelnachweise

Bearbeiten
  1. Rosetta Code: QR decomposition