Gram-Schmidt-Orthogonalisierung

Wähle, wie nah die Spalten beieinanderliegen. Lies dann die beiden Raster. Sie zeigen dir den Winkel zwischen jedem Vektorpaar, das die beiden Versionen zurückgeben.

Rechenbeispiel — Gram-Schmidt-Orthogonalisierung

Konditionszahl von A	1,7321 × 10⁸
Rest von Spalte 3	1,4142 × 10⁻⁸
Was q₂ von Spalte 3 abzieht	0
Größtes Skalarprodukt	0,5000
Winkel zwischen q₂ und q₃	60,00°
Orthogonalitätsverlust klass.	0,7071
Orthogonalitätsverlust modif.	1,1547 × 10⁻⁸
Orthogonalitätsverlust 2x	2,6158 × 10⁻²⁴

Das sind die Werte, die der Rechner beim Öffnen der Seite anzeigt. Ändere oben eine Eingabe, um sie neu zu berechnen.

Interaktive Simulation wird geladen...

Klassisches Gram-Schmidt liefert zwei Vektoren im Winkel von 60° 🖖

Drücke auf ε = 10⁻⁸, klassisch. Im klassischen Raster zeigen zwei Zellen 60,00°. Das ist das Paar q₂ und q₃ in beiden Ecken. Jeder andere Winkel in beiden Rastern liegt bei 90,00°. Das Skalarprodukt dahinter beträgt 0,5000. Ein Halb, bis auf eine Einheit an der letzten Stelle bei doppelter Genauigkeit.

Die Matrix besteht aus drei Spalten. Jede davon ist der Vektor (1, 0, 0, 0), bei dem ein einzelnes 10⁻⁸ in einem anderen Eintrag platziert wurde. Peter Läuchli nutzte sie 1961, und seitdem trägt sie seinen Namen. Beinahe identische Spalten sind der unbequeme Fall für jede Fehlerquadrat-Ausgleichsrechnung. Sie tauchen immer dann auf, wenn zwei deiner Messungen fast dieselbe Messung sind.

Das modifizierte Gram-Schmidt-Verfahren erhält dieselben drei Spalten. Jeder zurückgegebene Winkel beträgt 90,00°, bei exakt gleicher Anzahl an Multiplikationen.

Schiebe den Abstand eine Stufe zurück auf 10⁻⁷. Diese Zelle zeigt nun 89,24°. Ein einziger Klick auf den Schieberegler macht den Unterschied, ob du um ein Dreiviertelgrad oder um dreißig Grad danebenliegst.

Die beiden Versionen berechnen eine einzige Zahl unterschiedlich 🖖

Stell den Durchlauf auf Spalte 3 und klick auf Klassisch, dann auf Modifiziert. Die Zeile Was q₂ aus Spalte 3 abzieht springt von 0 auf 7,0711 × 10⁻⁹. Bei drei Spalten ist das der einzige Koeffizient, den beide Algorithmen unterschiedlich berechnen. Der Diagonaleintrag darunter und jeder Winkel im Raster ergeben sich unmittelbar daraus.

Das klassische Verfahren berechnet diesen Koeffizienten als q₂ᵀa₃, anhand der Spalte, die es erhalten hat. Das modifizierte Verfahren berechnet q₂ᵀv, anhand des Rests der Spalte, nachdem der Anteil von q₁ abgezogen wurde. Auf dem Papier sind beide Zahlen gleich. v unterscheidet sich von a₃ nur um ein Vielfaches von q₁, und q₂ steht senkrecht auf q₁. In der Rechenpraxis ist q₂ jedoch nicht exakt orthogonal zu q₁, und a₃ hat immer noch einen Anteil von 7,0711 × 10⁻⁹, der in Richtung q₂ zeigt. Das klassische Verfahren schaut dort gar nicht erst hin. Also behält q₃ diesen Anteil, und genau das sind die 60°.

Zieh den Schieberegler von 10⁻¹ auf 10⁻⁸ und beobachte, wie sich die beiden Kurven im unteren Feld voneinander trennen. Legt man eine Ausgleichskurve durch die 113 Punkte in diesem Bereich, wächst der Orthogonalitätsverlust auf diesen drei Spalten beim klassischen Verfahren mit κ^1,98 und beim modifizierten mit κ^0,98. Das Feld berechnet die Kurvenanpassung für jede gewählte Spaltenzahl neu. Diese beiden Exponenten verschieben sich also, wenn du die Matrix vergrößerst. Eine einzige umgestellte Subtraktion ist eine ganze Potenz der Konditionszahl wert.

Die Basis ist ruiniert, doch die Faktorisierung ist exakt 🖖

Die letzte Zeile der Ergebnistabelle lautet ‖A − QR‖ / ‖A‖. Bei ε = 10⁻⁸ steht dort unter Klassisch eine 0. Multipliziere das ruinierte Q mit dem R, das daneben herauskam. Die ursprüngliche Matrix taucht wieder auf, Eintrag für Eintrag, ohne jeden Rest.

Geht man alle Zustände durch, die mit den Reglern erreichbar sind – vier Spaltenzahlen, zehn Abstände und beide Algorithmen –, steigt dieses relative Residuum nie über 5,5511 × 10⁻¹⁷. Der Abstand zwischen 1 und der nächsten Double-Zahl beträgt 2,2204 × 10⁻¹⁶. Es bleibt also bei einem Viertel einer Einheit in der letzten Stelle.

Gram-Schmidt hat A tatsächlich als Produkt der zurückgegebenen Vektoren und der notierten Koeffizienten aufgeschrieben. Was es nicht mehr tat, war, Vektoren im rechten Winkel zueinander zurückzugeben. Welche der beiden Eigenschaften du brauchst, entscheidet dein eigener Code: Die Projektion eines Vektors mit Qᵀb stützt sich auf jene Eigenschaft, die versagt hat. Die Rekonstruktion von A aus Q und R stützt sich auf jene, die intakt blieb.

Aufgabe vollständig gelöst

Wo das klassische Gram-Schmidt-Verfahren seinen rechten Winkel verliert 8 Schritte
Das Werkzeug zeigt drei Spalten im Abstand von 10⁻⁸ und die Lösung des klassischen Gram-Schmidt-Verfahrens dafür. Rechne den Algorithmus von Hand durch. Finde jene einzelne Rundung, die das Verfahren dreißig Grad kostet. Erkläre dann, wie weit jede der beiden Versionen getrieben werden kann, bevor ein Einsatz nicht mehr lohnt.

Öffne diese Aufgabe im Tool
1. Jede Spalte ist der Vektor (1, 0, 0, 0), bei dem ein ε in einen anderen Eintrag gesetzt wurde. AᵀA ist die 3 × 3-Einsmatrix plus ε²I. Ihre Eigenwerte sind einmal 3 + ε² und zweimal ε². Die Singulärwerte sind folglich √(3 + ε²) und ε. Die Konditionszahl steht damit exakt fest und ist nicht nur geschätzt. Die oberste Zeile der Tabelle zeigt 1,7321 × 10⁸.
  
  $A = 1 ε 00 10 ε 0 100 ε, ε = 1 0^{- 8}, κ (A) = \frac{3 + ε ^{2}}{ε} = 1.7321 \times 1 0^{8}$
2. ‖a₁‖² ist 1 + ε², und ε² ist 10⁻¹⁶. Ein Double-Wert fasst 1, und die nächste darstellbare Zahl liegt 2,2204 × 10⁻¹⁶ weiter hinten. Alles unterhalb der Hälfte dieser Lücke, also 1,1102 × 10⁻¹⁶, verschwindet bei der Addition zu 1. Zehn hoch minus sechzehn liegt darunter. Die Quadratwurzel ist 1. Ein Teilen dadurch ändert nichts. Somit kommt q₁ als die Spalte selbst heraus, bei der das ε² einfach aus den Akten gestrichen wurde. Jede andere Rundung in diesem Problem kostet eine Einheit an der letzten Stelle. Diese hier kostet dreißig Grad.
  
  $∥ a_{1} ∥^{2} = 1 + ε^{2} = 1 + 1 0^{- 16}, 1 0^{- 16} < 1.1102 \times 1 0^{- 16} ⟹ fl (1 + ε^{2}) = 1, q_{1} = a_{1}$
3. Als Nächstes Spalte 2. Hier sind sich beide Versionen einig, da nur ein einziger Vektor davor liegt. Der Projektionskoeffizient ist q₁ᵀa₂ = 1. Zieht man ihn ab, bleibt (0, −ε, ε, 0). Das hat die Länge ε√2 = 1,4142 × 10⁻⁸. Spalte 2 kam mit der Länge 1 an, soweit ein Double das beurteilen kann. Bis auf 1,4142 × 10⁻⁸ war sie bereits vollständig durch q₁ abgedeckt.
  
  $r_{12} = q_{1}^{T} a_{2} = 1, v_{2} = a_{2} - q_{1} = (0, - ε, ε, 0)^{T}, ∥ v_{2} ∥ = ε 2 = 1.4142 \times 1 0^{- 8}$
4. Spalte 3, auf klassische Art. Der Koeffizient zu q₁ ist wieder 1. Der Koeffizient zu q₂ ist q₂ᵀa₃. Weil q₂ komplett in den Einträgen 2 und 3 lebt, während a₃ sein ε im Eintrag 4 behält, beträgt dieses Skalarprodukt exakt 0. Genau das besagt die Zeile Was q₂ von Spalte 3 abzieht. Zieht man nur q₁ ab, bleibt (0, −ε, 0, ε), abermals mit der Länge ε√2. Daher zeigt Rest von Spalte 3 ebenfalls 1,4142 × 10⁻⁸.
  
  $r_{13} = q_{1}^{T} a_{3} = 1, r_{23} = q_{2}^{T} a_{3} = 0, v_{3} = a_{3} - q_{1} = (0, - ε, 0, ε)^{T}, q_{3} = \frac{1}{2} (0, - 1, 0, 1)^{T}$
5. Sowohl q₂ als auch q₃ wurden gebildet, indem man q₁ von einer Spalte mit der gleichen führenden 1 abzog. Beide tragen also ein −1/√2 in Eintrag 2 und haben sonst nichts gemeinsam. Ihr Skalarprodukt ist ein Halb. Das Werkzeug druckt 0,5000 und 60,00°. Der intern gespeicherte Wert ist 0,4999999999999999. Er weicht um eine Einheit an der letzten Stelle von dem Halb ab, das die Algebra vorgibt, weil ±1/√2 kein Binärbruch ist. Das ist die Größenordnung des Fehlers, zu dem diese Arithmetik fähig ist. Die dreißig Grad sind es nicht. Die klassische Verlustzeile zeigt 0,7071, was √(2 × 0,5²) entspricht. Das ist dasselbe Halb, das in jeder Ecke von QᵀQ einmal gezählt wird.
  
  $q_{2}^{T} q_{3} = \frac{1}{2} [(- 1) (- 1) + (1) (0)] = \frac{1}{2}, arccos \frac{1}{2} = 6 0^{\circ}, 2 \cdot 0 . 5^{2} = 0.7071$
6. Das modifizierte Verfahren an derselben Spalte. Nachdem der Anteil von q₁ abgezogen wurde, lautet der Rest (0, −ε, 0, ε). Diese Version misst nun hiergegen anstatt gegen a₃. Bildet man das Skalarprodukt mit q₂, ergibt das ε/√2 = 7,0711 × 10⁻⁹, also nicht null. Zieht man das ab, verbleibt (0, −ε/2, −ε/2, ε) mit der Länge ε√(3/2) = 1,2247 × 10⁻⁸. Damit kommt q₂ᵀq₃ auf 1,1102 × 10⁻¹⁶ statt eines Halben. Die modifizierte Verlustzeile zeigt 1,1547 × 10⁻⁸. Darin steckt nur noch das, was q₁ nicht entfernen konnte. Die Subtraktion lässt −ε² in Richtung q₁ liegen. Teilt man dies durch einen Rest, der nur ε√2 lang ist, wird daraus wieder ε/√2.
  
  $r_{23} = q_{2}^{T} v_{3} = \frac{ε}{2} = 7.0711 \times 1 0^{- 9}, v_{3} \mapsto (0, - \frac{ε}{2}, - \frac{ε}{2}, ε)^{T}, ∥ v_{3} ∥ = ε \frac{3}{2} = 1.2247 \times 1 0^{- 8}, q_{2}^{T} q_{3} = 0$
7. Dasselbe Argument greift bei jeder Spaltenanzahl, die das Werkzeug bietet. Das klassische Verfahren bildet jedes spätere q, indem es q₁ von einer Spalte mit der gleichen führenden 1 abzieht. Somit landet jedes Paar davon bei 60°. Das sind exakt (n − 1)(n − 2)/2 Paare. Der Orthogonalitätsverlust entspricht der Quadratwurzel aus der Hälfte dieser Anzahl. Bei zwei Spalten existieren solche Paare nicht. Übrig bleibt lediglich die ε-große Neigung, die beide Versionen teilen.
  
  $n 2345 κ (A) 1.4142 \times 1 0^{8} 1.7321 \times 1 0^{8} 2.0000 \times 1 0^{8} 2.2361 \times 1 0^{8} # {6 0^{\circ}} 0136 ∥ Q_{C}^{T} Q_{C} - I ∥_{F} 1.0000 \times 1 0^{- 8} 0.7071 1.2247 1.7321 ∥ Q_{M}^{T} Q_{M} - I ∥_{F} 1.0000 \times 1 0^{- 8} 1.1547 \times 1 0^{- 8} 1.2247 \times 1 0^{- 8} 1.2649 \times 1 0^{- 8}$
8. Die klassische Ausgleichsgerade erreicht einen Verlust von 1 bei κ = 2,35 × 10⁸. Genau deshalb verläuft die klassische Kurve am rechten Rand des Feldes bereits flach. Die modifizierte Version erhält eine gemessene Schranke anstelle einer extrapolierten. Ihre eigene Kurve knickt hinter 10⁻⁸ nämlich um. Eine so weit hinausgezogene Ausgleichsgerade würde durch das Feld darunter sofort widerlegt. Über 724 Punkte des Durchlaufs hinweg überstieg ihr Verlust nie das 1,25-Fache der Maschinengenauigkeit mal κ. Diese Schranke erreicht die 1 erst bei κ = 3,60 × 10¹⁵. Eins geteilt durch die Maschinengenauigkeit ergibt 4,50 × 10¹⁵.
  
  $lo g_{10} κ = \frac{16 . 5755}{1 . 9801} 1.25 \cdot (2.2204 \times 1 0^{- 16}) κ (2.2204 \times 1 0^{- 16})^{- 1} = = = 8.3710 ⟹ κ = 2.35 \times 1 0^{8} 1 ⟹ κ = 3.60 \times 1 0^{15} 4.50 \times 1 0^{15}$
Antwort

Die dreißig Grad stammen aus einer einzigen Zeile. Klassisch berechnet q₂ᵀa₃ = 0, während modifiziert q₂ᵀv₃ = 7,0711 × 10⁻⁹ liest. Auf dem Papier sind das dieselben Zahlen.

Der Abstand von 10⁻⁸ ist die erste Stufe, bei der ε² aus 1 + ε² verschwindet. Schließlich ist die Wurzel aus der halben Maschinengenauigkeit exakt 1,0537 × 10⁻⁸. Von dort an abwärts zeigt die Zelle 60,00°, egal wie viel weiter du den Regler ziehst. Eine Stufe zurück, bei 10⁻⁷, misst sie 89,24° und das klassische Verfahren hat bereits 1,88 × 10⁻² an Orthogonalität eingebüßt. Der Kollaps beginnt also lange bevor der Schieberegler die Stufe erreicht, die für glatte Werte sorgt.

Åke Björck veröffentlichte die Analyse der modifizierten Version 1967. Seitdem ist die geänderte Reihenfolge der Standard. Beide führen dieselbe Anzahl an Multiplikationen aus. Bei diesem Abstand weicht der Orthogonalitätsverlust zwischen ihnen um den Faktor 6,1 × 10⁷ ab. Die modifizierte Variante behält ihre rechten Winkel bis zu dem Punkt bei, an dem ein Double die Spalten überhaupt nicht mehr voneinander unterscheiden kann. Das klassische Verfahren gibt sie sieben Größenordnungen früher auf.

$q_{2}^{T} a_{3} = 0 \neq = q_{2}^{T} v_{3} = \frac{ε}{2} = 7.0711 \times 1 0^{- 9}$

Quellen (5)

The matrix the tool is built on, and the name it carries: P. Läuchli, "Jordan-Elimination und Ausgleichung nach kleinsten Quadraten", Numerische Mathematik 3(1), 226–240, 1961. doi:10.1007/BF01386022
The analysis of the modified variant that made it the default: Åke Björck, "Solving linear least squares problems by Gram-Schmidt orthogonalization", BIT 7(1), 1–21, 1967. doi:10.1007/BF01934122
The corresponding rounding-error analysis of the classical variant: Luc Giraud, Julien Langou, Miroslav Rozložník and Jasper van den Eshof, "Rounding error analysis of the classical Gram-Schmidt orthogonalization process", Numerische Mathematik 101(1), 87–100, 2005. doi:10.1007/s00211-005-0615-4
Gram's paper of 1883, where the process is set out for least squares: J. P. Gram, "Ueber die Entwickelung reeller Functionen in Reihen mittelst der Methode der kleinsten Quadrate", Journal für die reine und angewandte Mathematik 1883(94), 41–73. doi:10.1515/crll.1883.94.41
And Schmidt's of 1907, where the form used today appears: Erhard Schmidt, "Zur Theorie der linearen und nichtlinearen Integralgleichungen", Mathematische Annalen 63, 433–476, 1907. doi:10.1007/BF01449770

Vom 42Math-Team · Anhand veröffentlichter Quellen geprüft. Korrekturen willkommen.