~2025-43445-68: Möglicher Tippfehler im Beispiel. Sollte "a" sein, da dort gerade "v" ausgerechnet wird

2025-12-27T15:04:10Z

Möglicher Tippfehler im Beispiel. Sollte "a" sein, da dort gerade "v" ausgerechnet wird

Neue Seite

In der [[Mathematik]] beschreibt die '''Householdertransformation''' die [[Spiegelung (Geometrie)|Spiegelung]] eines [[Vektor]]s an einer [[Hyperebene]] durch Null im [[Euklidischer Raum|euklidischen Raum]]. Im dreidimensionalen Raum ist sie somit eine Spiegelung an einer [[Ebene (Mathematik)|Ebene]] (durch den Ursprung). Die Darstellung dieser [[lineare Abbildung|linearen Abbildung]] durch eine Matrix wird als '''Householder-Matrix''' bezeichnet. Verwendung findet sie vor allem in der [[Numerische Mathematik|numerischen Mathematik]], wenn mittels [[Orthogonale Matrix|orthogonaler Transformationen]] Matrizen so gezielt umgeformt werden, dass bestimmte Spaltenvektoren auf das Vielfache des ersten [[Einheitsvektor]]s 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 ==

[[Datei:Householdertransformation.svg|miniatur|rechts|Illustration der Householder-Transformation in 2D. Der Schnittpunkt der Spiegelebene mit <math>\vec{v}</math> ist gegeben durch <math>\vec{x}-\langle \vec{v},\, \vec{x} \rangle\vec{v}</math>.]]
Die Spiegel-Hyperebene kann durch einen [[Normalenvektor]] <math>v</math>, also einen Vektor, der [[orthogonal]] zur Hyperebene ist, definiert werden. Ist <math>v</math> als [[Spaltenvektor]] gegeben und <math>I</math> die [[Einheitsmatrix]], dann wird die oben beschriebene lineare Abbildung durch die folgende Matrix dargestellt:

:<math>H = I - \frac {2} {v^T v} v v^T.</math>

Dabei bezeichnet <math>v^T</math> die [[Transponierte]] des Spaltenvektors <math>v</math>, also einen Zeilenvektor.
Der Nenner <math>v^T v</math> ist das [[Skalarprodukt]] von <math>v</math> mit sich selbst, <math>v v^T</math> das [[dyadisches Produkt|dyadische Produkt]]. Die Matrix <math>\tfrac {1} {v^T v} v v^T</math> beschreibt die [[Orthogonalprojektion]] auf die durch <math>v</math> gegebene Richtung.
Ist <math>v</math> auf die Länge eins normiert, also <math>v^T v=1</math>, so vereinfacht sich die Formel zu

:<math>H = I - 2 v v^T.</math>

Die Spiegelungseigenschaft ersieht man daraus, dass
:<math>
Hx
= x - 2 v v^Tx
= x - 2 \langle v,\, x \rangle\, v
= (x-\langle v,\, x \rangle\, v)-\langle v,\, x \rangle\, v
</math>,
wobei <math>\langle \cdot, \cdot \rangle</math> das [[Standardskalarprodukt]] bezeichnet. Der Term <math>\langle v,\, x \rangle</math> entspricht dabei dem Abstand des Punktes <math>x</math> zur Hyperebene <math>v^\perp</math>, die orthogonal zu <math>v</math> ist. Der Vektor <math>x</math> wird also zuerst um <math>\langle v,\, x \rangle\, v</math> verschoben, d. h. es wird auf die Hyperebene <math>v^\perp</math> projiziert und anschließend erneut um <math>\langle v,\, x \rangle\, v</math> auf den gespiegelten Wert verschoben. Unter der Spiegelung wird der Anteil in der Ebene invariant gelassen, der Anteil in Richtung <math>v</math>, also senkrecht zur Ebene, wird „umgeklappt“, also nun abgezogen statt addiert.

Die Householder-Matrix hat folgende Eigenschaften:
* Sie ist [[Symmetrische Matrix|symmetrisch]]: <math>H = H^T</math>
* Sie ist [[Orthogonale Matrix|orthogonal]]: <math>H^T H = I</math>
* Sie ist involutorisch: <math>H^2 = I</math> (Dies folgt aus der Symmetrie und der Orthogonalität.)
* Sie hat den einfachen [[Eigenwert]] −1 zum Eigenvektor <math>v</math> und den <math>n-1</math>-fachen Eigenwert 1. Der Eigenraum zum Eigenwert 1 ist die Spiegelebene, also das [[Orthogonales Komplement|orthogonale Komplement]] des von <math>v</math> erzeugten eindimensionalen Unterraums.
* [[Matrix-Vektor-Multiplikation]]en mit <math>H</math> sind schnell berechenbar.

== Konstruktion einer spezifischen Spiegelung ==

Es sei ein Vektor <math>a</math> gegeben, der auf ein Vielfaches des Vektors <math>e</math> gespiegelt werden soll, das heißt, gesucht ist ein Einheitsvektor <math>v</math>, so dass mit der zugehörigen Householder-Matrix <math>H_v = I - 2 vv^T</math> gilt <math>H_va=\lambda e</math>. Geometrisch ist der Vektor <math>v</math> die Richtung einer der zwei Winkelhalbierenden der Geraden in Richtung <math>a</math> und in Richtung <math>e</math>. 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 <math>v</math>, der Vektor <math>v</math> selbst ergibt sich durch Normieren dieser Richtung. Die zweite Winkelhalbierende ergibt sich, indem die Konstruktion ausgehend von <math>a</math> und <math>-e</math> durchgeführt wird.

Der Einfachheit halber sei <math>e</math> normiert, <math>\|e\|=1</math>. Dann muss, wegen der Orthogonalität der Spiegelung, <math>\lambda=\pm\|a\|</math> gelten. Der gesuchte Spiegelungsvektor <math>v</math> ergibt sich nun durch Normieren des Differenzvektors <math>\tfrac12(a-\lambda\,e)</math>, also
:<math>v=\frac{a-\lambda\, e}{\|a-\lambda\, e\|}</math>.
Beide Vorzeichenvarianten führen zum gewünschten Ergebnis (sofern der Nenner von Null verschieden ist). Aus Gründen numerischer Stabilität wird das [[Vorzeichen (Zahl)|Vorzeichen]] von <math>\lambda</math> so gewählt, dass der Nenner am größten ist, also <math>\lambda\cdot\langle a,\,e\rangle\leq 0</math> gilt.

In der Probe ergibt sich
:<math>\begin{align}
H_va
& = a-2v\langle v,a\rangle
\\[0.5em]
& = a-2\,(a-\lambda e)\,\frac{
\|a\|^2-\lambda\langle e,a\rangle
}{
\|a\|^2-2\lambda\langle a,e\rangle+\lambda^2\|e\|^2}
\\[0.8em]
& = a-2\,(a-\lambda e)\,
\frac{\|a\|^2-\lambda\langle e,a\rangle}{\|a\|^2-2\lambda\langle a,e\rangle+\|a\|^2}
\\[0.5em]
& = \lambda e\\[0.5em]
\end{align}</math>

=== Beispiel ===

Am häufigsten wird der Fall betrachtet, in dem <math> e=e_1</math> der erste kanonische Basisvektor ist. Sei <math>a=(a_1,\bar a)</math> in erste Komponente und Restvektor zerlegt. Dann gilt für die Norm <math>\|a\|=\sqrt{|a_1|^2+\|\bar a\|^2}</math>. Als Vorzeichen von <math>\lambda</math> ist das Vorzeichen von <math>-a_1</math> zu wählen, die Richtung der Spiegelung ist dann
:<math>a-\lambda e_1= a + \operatorname{sign}(a_1)\,\|a\| e_1 = \Bigl(\operatorname{sign}(a_1)\,(|a_1|+\|a\|),\;\bar a\Bigr)</math>.
Dabei ist <math>\operatorname{sign}(a_1) := \begin{cases}1 &\quad\text{für}\quad a_1 \geq 0, \\ -1 &\quad\text{für}\quad a_1 < 0. \end{cases}</math>

Der Vektor <math>v</math> entsteht durch Normierung dieser Richtung. Nach Umformen stellt sich die Norm der Richtung als
:<math>\|a-\lambda e_1\|=\sqrt{2\,\|a\|\,(\|a\|+|a_1|)}</math>
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 ==

Householder-Spiegelungen können zur [[Stabilität (Numerik)|stabilen]] Berechnung von [[QR-Zerlegung]]en einer Matrix <math>A=QR\in\R^{m\times n}</math> verwendet werden, indem zunächst die erste Spalte der Matrix mit einer Spiegelung <math>H_1</math> 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 <math>H_1A</math> mit einer Spiegelung <math>H_2</math> 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 <math>i</math>-ten Schritt wird also die Untermatrix unter der Position <math>(i,i)</math> des [[Matrizenprodukt|Produkts]] <math>A_i=H_{i-1}\cdots H_1A</math> auf die gleiche Art reduziert, bis die Restmatrix <math>H_n\cdots H_2 H_1A = R</math> Dreiecksgestalt besitzt.

Mit <math>Q^T = H_n\cdots H_2 H_1</math> gilt <math>Q^T A = R</math>, also ergibt sich die QR-Zerlegung
:<math> A=QR</math>
mit
:<math>Q= (H_n\cdots H_2 H_1)^T = H_1 H_2\cdots H_{n} \in\R^{m\times m},\quad\quad R=\begin{pmatrix}R_1\\0\end{pmatrix}\in\R^{m\times n}.</math>

Man beachte, dass <math>Q</math> hier eine quadratische Matrix ist. Meist werden die Matrizen <math>Q</math> bzw. <math>Q^T</math> nicht explizit berechnet, sondern man nutzt direkt die Produktform. Dazu werden die Spiegelvektoren <math>v_i</math> von <math>H_i=H_{v_i}</math> im frei gewordenen Platz der Matrix <math>A</math> gespeichert.

Die Zahl der Operationen für die QR-Zerlegung einer Matrix <math>A=QR \in \R^{m\times n},{m \ge n}</math> mit dem Householder-Verfahren beträgt:
:<math>\begin{matrix} m \gg n &\qquad& 2n^2 \cdot m \\ m \approx n &\qquad& \frac{4}{3} n^3 \end{matrix}
</math>
Im Fall <math>m=n</math> genügen <math>n-1</math> Schritte, da die letzte Spalte nicht mehr transformiert werden muss.

=== Pseudocode ===
Da für die meisten Berechnungen das explizite Ausrechnen von <math>Q</math> nicht nötig ist, reicht es, nur die Matrix <math>R</math> zu berechnen. <math>z</math> ist die linke Spalte der jeweiligen [[Untermatrix]]. Bei der unten angegebenen Funktion wird das Ergebnis direkt in <math>A</math> geschrieben, so dass nach Abarbeitung des Algorithmus das <math>R</math> in <math>A</math> steht. Die Zeile <math>R=A</math> 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 <math>Q</math> 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 ==
Das folgende Beispiel in der [[Programmiersprache]] [[C++]] zeigt die Implementierung der Householdertransformation. Bei der Ausführung des Programms wird die [[Funktion (Mathematik)|Funktion]] ''main'' verwendet, die das Ergebnis auf der Konsole ausgibt.<ref>Rosetta Code: [https://rosettacode.org/wiki/QR_decomposition QR decomposition]</ref><syntaxhighlight lang="c++">
#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
}
}
</syntaxhighlight>

== Siehe auch ==
* [[Spiegelungsmatrix]]

== Literatur ==
* Alston Householder: ''Unitary triangularization of a nonsymmetric matrix''. In: ''Journal of the ACM''. Band 5, Nr. 4, 1958, S. 339–342.
* [[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 (Mathematiker)|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 ==
<references />

[[Kategorie:Numerische lineare Algebra]]
[[Kategorie:Transformation]]

Householdertransformation - Versionsgeschichte

~2025-43445-68: Möglicher Tippfehler im Beispiel. Sollte "a" sein, da dort gerade "v" ausgerechnet wird