Mandalex

Vorsicht, der Inhalt dieser Seite ist mathematisch nicht 100%-ig korrekt. Besonders die Variablen sind nur gerade so gewählt, dass man die Herleitungen versteht.

Die Methode der kleinsten Quadrate (auf Englisch "linear least squares") ist dazu da, um mittels Messwerten die Parameter einer Formel so gut wie nur möglich anzunähern. Mathematisch ausgedrückt ist dies eine orthogonale Projektion eines Kolonnenvektors auf einen Vektorraum. Dies bedeutet, dass der Abstand dieser Projektion (die ermittelten Parameter) und des durch den Kolonnenvektor gegebenen Punktes (die Gesamtheit der Messwerte) minimal unter allen möglichen Lösungen ist.

Herleitung:

Man stelle sich dazu folgendes vor: Gegeben sei ein Vektorraum A (Beispielsweise ein Blatt Papier: 2 Dimensionen). Nun erhält man den Auftrag, zu einen beliebigen Punkt b irgendwo im gesamten Universum (das heisst: definiert durch beliebig viele Dimensionen, in diesem Beispiel die 3 Raum-Dimensionen) einen entsprechenden Punkt x im Vektorraum A (also auf dem Papier) zu finden, sodass der Abstand (allgemein ausgedrückt der Kern der Abbildung) des Punktes x und des Punktes p orthogonal (senkrecht) auf dem Vektorraum A steht.

Einfach ausgedrückt: Wenn man das Papier A flach auf einen Tisch legt und von dem Punkt b einen Filzstift auf dieses Papier fallen lässt, so ist der dadurch entstandene Punkt der gesuchte Punkt x.

Mit diesem Beispiel geht es nun weiter: Das Blatt Papier wird definiert durch eine Matrix (in diesem Fall eine 3x2-Matrix), wodurch genau festgelegt wird, wie sie im dreidimensionalen Raum liegt. Normalerweise sind die einzelnen Vektoren linear unabhängig, weslhalb der andere Fall nun nicht weiter betrachtet wird. Der Punkt b wird definiert durch einen 3-dimensionalen Vektor, welcher angibt, an welchem absoluten Punkt des gesamten Universums sich dieser Punkt befindet. Der Punkt x ist wiederum ein Vektor, welcher jedoch nur 2 Dimensionen hat und relativ zum Ursprung des Papiers genau angibt, wo auf diesem Papier sich dieser Punkt x befindet. Der absolute Punkt, welcher durch x beschrieben wird, errechnet sich durch die einfache Formel A·x.

Der einfache Fall

Wie man nun zur korrekten Lösung kommt: Man stelle sich vor, der Punkt x befindet sich bereits auf dem Papier A (Wie im folgenden Beispiel). Dann lautet die korrekte Formel wie bereits angetönt:

A·x = b

Dies bedeutet, dass b die Formel A·x löst, und somit bedeutet dies auch, dass der Abstand der beiden Punkte x und b gleich 0 ist. Die korrekte Formel für den Abstand (oder Fehler) E lautet:

E = ‖A·x − b‖

wobei A·x − b der Fehlervektor genannt wird, welcher senkrecht zum Papier steht. Zurück zum Beispiel. Wir geben konkrete Werte an und erhalten folgende Matrizen:

⎛ 1 4 ⎞ ⎛ 7 ⎞ ⎛ 7 ⎞ A = ⎜ 2 5 ⎟ b = ⎜ 8 ⎟ x = ⎛−1 ⎞ A·x = ⎜ 8 ⎟ ⎝ 3 6 ⎠ ⎝ 9 ⎠ ⎝ 2 ⎠ ⎝ 9 ⎠

Der weniger einfache Fall

Wenn sich nun der Punkt b nicht zwingend auf dem Papier befindet, so wird E grössergleich 0 sein, da ja die Projektion A·x und der ursprüngliche Punkt b nicht gleichauf liegen. Dieser Fehler E muss minimiert werden. Um diese Minimierung zu erreichen, erinnert man sich daran, dass das innere Produkt zweier Vektoren, die senkrecht aufeinander stehen, 0 sein muss. Wir nehmen also einen beliebigen Vektor A·y, und bilden mit dem Fehlervektor A·x − b das innere Produkt, welches naturgemäss 0 ergibt:

(A·y)ʹ · (A·x − b) = 0

Daraus ergibt sich in der Folge:

yʹ·Aʹ · (A·x − b) = 0 yʹ · (Aʹ·A·x − Aʹ·b) = 0

Da y beliebig, muss gelten:

(Aʹ·A)·x − (Aʹ·b) = 0 (Aʹ·A)·x = Aʹ·b

und somit:

x = (Aʹ·A)⁻¹ · Aʹ·b

Diese letzte Formel nennt man die Methode der kleinsten Quadrate und bringt in jedem Fall (solange A linear unabhängig ist) das Ergebnis mit der minimalsten Abweichung zustande.

Projektionsmatrix und endgültige Lösung:

Die Lösung x ist nun erst der Punkt relativ zum Papier, welcher dem ursprünglichen Punkt b am nächsten liegt. A·x ist somit der absolute (angenäherte) Punkt, der gesucht wird. Er wird auch mit p bezeichnet. Die vollständige Formel für p = A·x lautet:

p = A·x = A · (Aʹ·A)⁻¹ · Aʹ·b

wobei P = A·(Aʹ·A)⁻¹·Aʹ die Projektionsmatrix genannt wird. Durch die Projektionsmatrix kann die absolute Lösung einfach beschrieben werden:

p = P · b

Wir kommen zurück zum Beispiel und lassen die Matrix A so, wie sie ist. Der Punkt b wird nun jedoch so gewählt, dass er nicht mehr auf dem Blatt Papier zu liegen kommt. Wir suchen somit die absoluten Koordinaten desjenigen Punktes p auf dem Papier, der diesem neuen Punkt b am nächsten kommt.

⎛ 1 4 ⎞ ⎛ 4 ⎞ A = ⎜ 2 5 ⎟ b = ⎜ 5 ⎟ ⎝ 3 6 ⎠ ⎝ 4 ⎠

Man kann den Punkt p direkt berechnen, indem man zuerst x berechnet:

Aʹ= ⎛ 1 2 3 ⎞ Aʹ·A = ⎛ 14 32 ⎞ (Aʹ·A)⁻¹ = ⎛ 1.4259 −0.5926 ⎞ ⎝ 4 5 6 ⎠ ⎝ 32 77 ⎠ ⎝−0.5926 0.2593 ⎠ Aʹ · b = ⎛ 26 ⎞ ⎝ 65 ⎠ x = (Aʹ·A)⁻¹ ·Aʹ·b = ⎛−1.4444 ⎞ ⎝ 1.4444 ⎠ ⎛ 4.3333 ⎞ p = A·x = ⎜ 4.3333 ⎟ ⎝ 4.3333 ⎠

Oder aber man berechnet direkt die Projektionsmatrix P und kommt darauf ebenfalls auf das korrekte Resultat:

⎛ 0.8333 0.3333 −0.1667 ⎞ P = A·(Aʹ·A)⁻¹ ·Aʹ = ⎜ 0.3333 0.3333 0.3333 ⎟ ⎝−0.1667 0.3333 0.8333 ⎠ ⎛ 4.3333 ⎞ p = P·b = ⎜ 4.3333 ⎟ ⎝ 4.3333 ⎠

Zum Schluss noch die Fehlerberechnung (wie weit ist der Punkt p vom ursprünglichen Punkt b entfernt?):

⎛ 0.3333 ⎞ ‖A·x − b‖ = 0.8165 A·x − b = ⎜−0.6667 ⎟ ⎝ 0.3333 ⎠

QR-Zerlegung:

Die Methode der kleinsten Quadrate kann auch mittels einer QR-Zerlegung der Matrix A erreicht werden (A = Q·R). Wenn diese Zerlegung bereits gemacht wurde, so kann die Lösung einfach errechnet werden durch die Formel:

x = R⁻¹ ·Qʹ·b

Pseudoinverse:

Grundsätzlich löst die Methode der kleinsten Quadrate nichts anderes als die Gleichung A·x = b. Da es sich bei Anwendung der kleinsten Quadrate jedoch oft um überbestimmte Gleichungssysteme handelt, kann man nicht einfach schreiben x = A⁻¹·b, denn A ist in einem überbestimmten Gleichungssystem nicht invertierbar. Wie jedoch oben hergeleitet wurde, ist x = (Aʹ·A)⁻¹·Aʹ·b, woraus man nun folgen kann, dass (Aʹ·A)⁻¹·Aʹ eine Art Inverse von A ist. Man nennt dies deshalb auch die Pseudoinverse:

A⁺ = (Aʹ·A)⁻¹ ·Aʹ

und somit: x = A⁺·b

Dies kann in Matlab so nicht geschrieben werden. Die Pseudoinverse kann mit pinv(A) errechnet werden (oder manuell). Dazu gerade noch die Anmerkung, dass Matlab die Pseudoinverse nicht mittels der oben aufgeführten Formel berechnet. Aufgrund stark singulärer Matrizen kann es vorkommen, dass bereits die Berechnung von (Aʹ·A)⁻¹ fehlschlägt. Matlab umgeht dies, indem es von vorne herein die Singulärwertzerlegung anwendet.

Singulärwert-Zerlegung:

Eine Matrix A kann stets folgendermassen unterteilt werden:

A = Q₁·Σ·Q₂ʹ

Q₁ ist eine mxm-orthogonal-Matrix
Σ ist eine Diagonalmatrix, welche mit den Singulärwerten von A gefüllt ist
Q₂ ist eine nxn-orthogonal-Matrix

Erinnerung: orthogonale Matrizen haben folgende Eigenschaften:

Q·Qʹ = I Qʹ = Q⁻¹

somit gilt auch Q·Q⁻¹ = I

Wie genau diese Zerlegung funktioniert, sei hier übersprungen. Sie ist in Matlab durch folgende Funktion gegeben: [Q1, S, Q2]=svd(A). Im Beispiel:

Q1 = ⎛ 0.4287 0.8060 0.4082 ⎞ ⎜ 0.5663 0.1124 −0.8165 ⎟ ⎝ 0.7039 −0.5812 0.4082 ⎠ S = ⎛ 9.5080 0 ⎞ ⎜ 0 0.7729 ⎟ ⎝ 0 0 ⎠ Q2 = ⎛ 0.3863 −0.9224 ⎞ ⎝ 0.9224 0.3863 ⎠

Ausgehend von der oben hergeleiteten Formel x = A⁺·b mit A⁺ als Pseudoinverse kommt man nun auf folgende überlegung: Wenn man für die Pseudoinverse (Q₁·Σ·Q₂ʹ)⁺ einsetzt, so ergibt sich folgendes:

(Q₁·Σ·Q₂ʹ)⁺ = ( (Q₁·Σ·Q₂ʹ)ʹ · (Q₁·Σ·Q₂ʹ) )⁻¹ · (Q₁·Σ·Q₂ʹ)ʹ = ( Q₂·Σʹ·Q₁ʹ · Q₁·Σ·Q₂ʹ )⁻¹ · Q₂·Σʹ·Q₁ʹ = ( Q₂·Σʹ · Σ·Q₂ʹ )⁻¹ · Q₂·Σʹ·Q₁ʹ = ( Q₂ʹ⁻¹·Σ⁻¹ · Σʹ⁻¹·Q₂⁻¹ ) · Q₂·Σʹ·Q₁ʹ = Q₂ · Σ⁻¹ · Σʹ⁻¹ · Q₂ʹ · Q₂ · Σʹ · Q₁ʹ = Q₂ · (Σʹ · Σ)⁻¹ · Σʹ · Q₁ʹ = Q₂ · Σ⁺ · Q₁ʹ

Somit entsteht folgende Formel: (Q₁·Σ·Q₂ʹ)⁺ = (Q₂·Σ⁺·Q₁ʹ)

Nun ist es glücklichen Umständen zu verdanken, dass S eine Diagonalmatrix ist, und dass die Pseudoinversen von Diagonalmatrizen sehr einfach zu berechnen sind: Von jedem Diagonaleintrag wird einfach der Kehrwert genommen: ⁻¹. Allerdings muss dabei noch folgendes definiert werden: Irgendeine Zahl durch 0 ist gleich 0. Somit fährt das obige Beispiel folgendermassen fort:

S⁺ = ⎛ 0.1052 0 0 ⎞ ⎝ 0 1.2939 0 ⎠ A⁺ = Q₂ * S⁺ * Q₁' = ⎛ −0.9444 −0.1111 0.7222 ⎞ ⎝ 0.4444 0.1111 −0.2222 ⎠

und somit:

x = A⁺ · b = ⎛−1.4444 ⎞ ⎝ 1.4444 ⎠

Und wiederum stimmt das Resultat mit dem bisher gesehenen überein.

Für die exakte Anwendung in Matlab, folgenden Befehl eingeben: help matfun