Methode der kleinsten Quadrate

Messpunkte und deren Abstand zu einer nach der Methode der kleinsten Quadrate bestimmten Funktion. Hier wurde eine logistische Funktion als Modellkurve gewählt.

Die Methode der kleinsten Quadrate (engl.: method of least squares) ist das mathematische Standardverfahren zur Ausgleichungsrechnung. Dabei wird zu einer Datenpunktwolke eine Kurve gesucht, die möglichst nahe an den Datenpunkten verläuft. Die Daten können physikalische Messwerte, wirtschaftliche Größen oder Ähnliches repräsentieren, während die Kurve aus einer parameterabhängigen problemangepassten Familie von Funktionen stammt. Die Methode der kleinsten Quadrate besteht dann darin, die Kurvenparameter so zu bestimmen, dass die Summe der quadratischen Abweichungen der Kurve von den beobachteten Punkten minimiert wird. Die Abweichungen werden Residuen genannt.

In der Beispielgrafik sind Datenpunkte eingetragen. In einem ersten Schritt wird eine Funktionenklasse ausgewählt, die zu dem Problem und den Daten passen sollte, hier eine logistische Funktion. Deren Parameter werden nun so bestimmt, dass die Summe der Quadrate der Abweichungen e der Beobachtungen y zu den Werten der Funktion minimiert wird. In der Grafik ist die Abweichung e an der Stelle x als senkrechter Abstand der Beobachtung y von der Kurve zu erkennen.

In der Stochastik wird die Methode der kleinsten Quadrate meistens als Schätzmethode in der Regressionsanalyse benutzt, wo sie auch als Kleinste-Quadrate-Schätzung bezeichnet wird. Angewandt als Systemidentifikation ist die Methode der kleinsten Quadrate in Verbindung mit Modellversuchen für Ingenieure ein Ausweg aus der paradoxen Situation, Modellparameter für unbekannte Gesetzmäßigkeiten zu bestimmen.

Geschichtliches

Carl Friedrich Gauß

Am Neujahrstag 1801 entdeckte der italienische Astronom Giuseppe Piazzi den Zwergplaneten Ceres. 40 Tage lang konnte er die Bahn verfolgen, dann verschwand Ceres hinter der Sonne. Im Laufe des Jahres versuchten viele Wissenschaftler erfolglos, anhand von Piazzis Beobachtungen die Bahn zu berechnen – unter der Annahme einer Kreisbahn, denn nur für solche konnten damals die Bahnelemente aus beobachteten Himmelspositionen mathematisch ermittelt werden. Der 24-jährige Gauß hingegen konnte auch elliptische Bahnen aus drei Einzelbeobachtungen berechnen. Da aber deutlich mehr Bahnpunkte vorlagen, wandte er seine Methode der kleinsten Quadrate an, um so die Genauigkeit zu erhöhen. Als Franz Xaver von Zach und Heinrich Wilhelm Olbers im Dezember 1801 den Kleinplaneten genau an dem von Gauß vorhergesagten Ort wiederfanden, war das nicht nur ein großer Erfolg für Gauß’ Methode: Piazzis Ruf, der aufgrund seiner nicht zu einer Kreisbahn passen wollenden Bahnpunkte stark gelitten hatte, war ebenfalls wieder hergestellt.^[1]

Piazzis Beobachtungen veröffentlicht in der Monatlichen Correspondenz vom September 1801

Die Grundlagen seines Verfahrens hatte Gauß schon 1795 im Alter von 18 Jahren entwickelt. Basis war eine Idee von Pierre-Simon Laplace, die Beträge von Fehlern aufzusummieren, so dass sich die Fehler zu Null addieren. Gauß nahm statt dessen die Fehlerquadrate und konnte die künstliche Zusatzanforderung an die Fehler weglassen. Unabhängig davon entwickelte der Franzose Adrien-Marie Legendre dieselbe Methode erstmalig im Jahre 1805 am Schluss eines kleinen Werkes über die Berechnung der Kometenbahnen^[2] und veröffentlichte eine zweite Abhandlung darüber im Jahr 1810. Von ihm stammt der Name Méthode des moindres carrés (Methode der kleinsten Quadrate).

1809 publizierte Gauß dann im zweiten Band seines himmelsmechanischen Werkes Theoria Motus Corporum Coelestium in sectionibus conicis solem ambientium (Theorie der Bewegung der Himmelskörper, welche in Kegelschnitten die Sonne umlaufen) sein Verfahren,^[3] inklusive der Normalgleichungen und des Gaußschen Eliminationsverfahrens^[4]. Dabei erwähnte er, dass er es schon vor Legendre entdeckt und benutzt habe, was zu einem Prioritätsstreit zwischen den beiden führte. Die Methode der kleinsten Quadrate wurde nun schnell das Standardverfahren zur Behandlung von astronomischen oder geodätischen Datensätzen.

Gauß benutzte dann das Verfahren intensiv bei seiner Vermessung des Königreichs Hannover durch Triangulation. 1821 und 1823 erschien die zweiteilige Arbeit sowie 1826 eine Ergänzung zur Theoria Combinationis observationum erroribus minimis obnoxiae (Theorie der den kleinsten Fehlern unterworfenen Kombination der Beobachtungen),^[5] in denen Gauß eine Begründung liefern konnte, wieso sein Verfahren im Vergleich zu den anderen so erfolgreich war: Die Methode der kleinsten Quadrate ist in einer breiten Hinsicht optimal, also besser als andere Methoden. Die genaue Aussage ist als der Satz von Gauß-Markow bekannt, da die Arbeit von Gauß wenig Beachtung fand und schließlich im 20. Jahrhundert von Andrei Andrejewitsch Markow wiederentdeckt und bekannt gemacht wurde. Theoria Combinationis enthält ferner wesentliche Fortschritte beim effizienten Lösen der auftretenden linearen Gleichungssysteme, wie das Gauß-Seidel-Verfahren und die LR-Zerlegung^[6].

Der französische Vermessungsoffizier André-Louis Cholesky entwickelte während des Ersten Weltkrieges die Cholesky-Zerlegung, die gegenüber den Lösungsverfahren von Gauß nochmal einen erheblichen Effizienzgewinn darstellte. In den 1960er Jahren entwickelte Gene Golub die Idee, die auftretenden linearen Gleichungssysteme mittels QR-Zerlegung zu lösen.

Das Verfahren

Voraussetzungen

Man betrachtet eine abhängige Größe y, die von einer Variablen x oder auch von mehreren Variablen beeinflusst wird. So hängt die Dehnung einer Feder nur von der aufgebrachten Kraft ab, der Gewinn eines Unternehmens jedoch von mehreren Faktoren wie Umsatz, den verschiedenen Kosten oder dem Eigenkapital. Zur Vereinfachung der Notation wird im Folgenden die Darstellung auf eine Variable x beschränkt. Der Zusammenhang zwischen y und den Variablen wird über eine Modellfunktion f, beispielsweise einer Parabel oder einer Exponentialfunktion

$y(x) = f(x;\alpha_1,\dots,\alpha_m)$,

die von x sowie von m Funktionsparametern α_j abhängt, modelliert. Diese Funktion entstammt entweder der Kenntnis des Anwenders; im Falle der Feder ist dies etwa das Hooksche Gesetz und damit eine lineare Funktion mit der Federkonstanten als einzigem Parameter. In schwierigeren Fällen, wie dem des Unternehmens kann der Wahl des Funktionstyps ein beliebig komplexer Modellierungsprozess vorausgehen, eventuell müssen auch verschiedene Modellfunktionen angesetzt werden und die Ergebnisse verglichen werden.

Die Parameter α_j dienen zur Anpassung des gewählten Funktionstyps an die beobachteten Werte y_i. Ziel ist es nun, die Parameter α_j so zu wählen, dass die Modellfunktion die Daten bestmöglich approximiert.

Um Informationen über die Parameter und damit die konkrete Art des Zusammenhangs zu erhalten, werden zu jeweils n gegebenen Werten x_i der unabhängigen Variablen x entsprechende Beobachtungswerte y_i (i = 1,...,n) erhoben.

Gauß und Legendre hatten die Idee, Verteilungsannahmen über die Messfehler dieser Beobachtungswerte zu machen. Sie sollten im Durchschnitt Null sein, eine gleichbleibende Varianz haben und von jedem anderen Messfehler stochastisch unabhängig sein. Man verlangt damit, dass in den Messfehlern keinerlei systematische Information mehr steckt, sie also rein zufällig um Null schwanken. Außerdem sollten die Messfehler normalverteilt sein, was zum einen wahrscheinlichkeitstheoretische Vorteile hat und zum anderen garantiert, dass Ausreißer in y so gut wie ausgeschlossen sind.

Um dies zu erfüllen, ist es notwendig, dass deutlich mehr Datenpunkte als Parameter vorliegen, es ist also n > m.

Minimierung der Summe der Fehlerquadrate

Das Kriterium zur Bestimmung der Approximation sollte so gewählt werden, dass große Abweichungen der Modellfunktion von den Daten stärker bestraft werden als kleine.

Dazu wird die Summe der Fehlerquadrate, die auch als Fehlerquadratsumme bezeichnet wird, als die Summe der quadrierten Differenzen zwischen den Werten der Modellkurve f(x_i) und den Daten y_i definiert. In Formelschreibweise mit $\vec{\alpha} = (\alpha_1, \alpha_2, \dots, \alpha_m) \in \mathbb{R}^m$ und $\vec{f} = (f(x_1,\vec{\alpha}), \dots, f(x_n,\vec{\alpha})) \in \mathbb{R}^n$ ergibt sich

$\sum_{i=1}^{n}(f(x_i, \vec{\alpha}) - y_i)^2 = \| \vec{f} - \vec{y} \|_2^2.$

Es sollen dann diejenigen Parameter α_j ausgewählt werden, bei denen die Summe der Fehlerquadrate minimal wird:

$\min_{\vec{\alpha}} \| \vec{f} - \vec{y} \|_2^2.$

Wie genau dieses Minimierungsproblem gelöst wird, hängt von der Art der Modellfunktion ab.

Lineare Modellfunktion

Lineare Modellfunktionen sind Linearkombinationen aus beliebigen, im Allgemeinen nicht-linearen Basisfunktionen. Für solche Modellfunktionen lässt sich das Minimierungsproblem analytisch über einen Extremwertansatz ohne iterative Annährungsschritte lösen. Zunächst werden einige einfache Spezialfälle und Beispiele gezeigt.

Spezialfall einer einfachen linearen Ausgleichsgeraden

Eine einfache Modellfunktion mit zwei linearen Parametern stellt das Polynom erster Ordnung

f(x) = α₀ + α₁x.

dar. Gesucht werden zu n gegebenen Messwerten $(x_1, y_1), \ldots, (x_n, y_n)$ die Koeffizienten α₀ und α₁ der bestangepassten Geraden. Die Abweichungen r_i zwischen der gesuchten Geraden und den jeweiligen Messwerten

$\begin{matrix} r_1 = \alpha_0 + \alpha_1 x_{1} - y_1\\ r_2 = \alpha_0 + \alpha_1 x_{2} - y_2\\ \vdots \\ r_n = \alpha_0 + \alpha_1 x_{n} - y_n\\ \end{matrix}$

nennt man Anpassungsfehler oder Residuen. Gesucht sind nun die Koeffizienten α₀ und α₁ mit der kleinsten Summe der Fehlerquadrate

$\min_{\alpha_0,\alpha_1} \sum_{i=1}^n r_i^2.$

Der große Vorteil des Ansatzes mit diesem Quadrat der Fehler wird sichtbar, wenn man diese Minimierung mathematisch durchführt: Die Summenfunktion wird als Funktion der beiden Variablen α₀ und α₁ aufgefasst (die eingehenden Messwerte sind dabei numerische Konstanten), dann die Ableitung (genauer: partielle Ableitungen) der Funktion nach diesen Variablen gebildet und von dieser Ableitung schließlich die Nullstelle gesucht. Da die Summenfunktion Quadrate der Variablen aufaddiert, wird bei der Ableitung aus diesen Quadraten einfach ein linearer Ausdruck, von dem die Nullstelle direkt und allgemein angegeben werden kann und schließlich nach den gesuchten Parametern auflösbar ist (längere Zwischenrechnung hier nicht dargestellt):

$\alpha_1 = \frac{\left(\sum\limits_{i=1}^n x_iy_i\right) - n \cdot \bar x \bar y}{\left(\sum\limits_{i=1}^n x_i^2\right) - n \cdot (\bar x)^2}$ und $\alpha_0 = \bar y - \alpha_1 \bar x$

mit $\bar x = \frac{1}{n} \sum\limits_{i=1}^n x_i$ als arithmetischem Mittel der x-Werte, $\bar y$ entsprechend. Die Lösung für α₁ kann mit Hilfe des Verschiebungssatzes auch als

$\alpha_1 = \frac{\sum\limits_{i=1}^n (x_i - \bar x)(y_i - \bar y)}{\sum\limits_{i=1}^n (x_i - \bar x)^2}$

angegeben werden.

Beispiel mit einer Ausgleichsgeraden

Streudiagramm von Längen und Breiten von 10 zufällig ausgewählten Kriegsschiffen

Folgendes Beispiel soll das Approximieren der linearen Funktion f(x) = α₀ + α₁x zeigen. Es wurden zufällig 10 Kriegsschiffe ausgewählt und bezüglich mehrerer Merkmale, darunter Länge (m) und Breite (m), analysiert. Es soll untersucht werden, ob die Breite eines Kriegsschiffs möglicherweise in einem festen Bezug zur Länge steht.

Das Streudiagramm zeigt, dass zwischen Länge und Breite eines Schiffs ein ausgeprägter linearer Zusammenhang besteht. Es wird also als Modellfunktion eine Ausgleichsgerade genommen und mit Hilfe der Methode der kleinsten Quadrate errechnet. Man erhält nun analog zum oben angegebenen Fall zunächst

$\bar x = \frac {\sum\limits_{i=1}^n x_i}{n} = \frac {1678}{10} = 167{,}8$

und entsprechend

$\bar y = 18{,}41.$

In der folgenden Tabelle sind die Daten zusammen mit den Zwischenergebnissen aufgeführt.

Nummer	Länge (m)	Breite (m)	$x_i- \overline x$	$y_i- \overline y$
i	xi	yi	$x_i^*$	$y_i^*$	$x_i^* \cdot y_i^*$	$x_i^*\cdot x_i^*$	$y_i^* \cdot y_i^*$
1	208	21,6	40,2	3,19	128,238	1616,04	10,1761
2	152	15,5	−15,8	−2,91	45,978	249,64	8,4681
3	113	10,4	−54,8	−8,01	438,948	3003,04	64,1601
4	227	31,0	59,2	12,59	745,328	3504,64	158,5081
5	137	13,0	−30,8	−5,41	166,628	948,64	29,2681
6	238	32,4	70,2	13,99	982,098	4928,04	195,7201
7	178	19,0	10,2	0,59	6,018	104,04	0,3481
8	104	10,4	−63,8	−8,01	511,038	4070,44	64,1601
9	191	19,0	23,2	0,59	13,688	538,24	0,3481
10	130	11,8	−37,8	−6,61	249,858	1428,84	43,6921
Σ	1678	184,1	0,0	0,00	3287,820	20391,60	574,8490

Damit bestimmt man α₁ als

$\alpha_1 = \frac{\sum\limits_{i=1}^n (x_i- \bar {x})(y_i - \bar y)}{\sum\limits_{i=1}^n (x_i- \bar x)^2} = \frac{3287{,}820} {20391{,}60} = 0{,}1612,$

so dass man sagen könnte, mit jedem Meter Länge wächst ein Kriegsschiff im Durchschnitt etwa 16 Zentimeter in die Breite. Das Absolutglied α₀ ergibt sich als

$\alpha_0 = \bar y - \alpha_1 \bar x = 18{,}41 - 0{,}1612 \cdot 167{,}8 = -8{,}6451.$

Die Anpassung der Punkte ist recht gut. Im Mittel beträgt die Abweichung zwischen der vorhergesagten Breite mit Hilfe des Merkmals Länge und der beobachteten Breite 2,1 m. Auch das Bestimmtheitsmaß, als normierter Koeffizient, ergibt einen Wert von ca. 92% (100% würde einer mitteleren Abweichung von 0 m entsprechen); zur Berechnung siehe das Beispiel zum Bestimmtheitsmaß.

Einfache polynomiale Ausgleichskurven

Streudiagramm: Durchschnittliches Gewicht von Männern nach Alter mit parabelförmiger Modellfunktion

Datensatz mit approximierenden Polynomen

Allgemeiner als eine lineare Ausgleichsgerade sind Ausgleichspolynome

$y(t) \approx \alpha_0 + \alpha_1 x + \alpha_2 x^2 + ... + \alpha_q x^q,$

die nun anhand eines Beispiels illustriert werden.

Als Ergebnisse der Mikrozensus-Befragung durch das statistische Bundesamt sind die durchschnittlichen Gewichte von Männern nach Altersklassen gegeben (Quelle: Statistisches Bundesamt, Wiesbaden 2009). Für die Analyse wurden die Altersklassen durch die Klassenmitten ersetzt. Es soll die Abhängigkeit der Variablen Gewicht (y) von der Variablen Alter (x) analysiert werden.

Das Streudiagramm lässt auf eine annähernd parabolische Beziehung zwischen x und y schließen, welche sich häufig gut durch ein Polynom annähern lässt. Es wird ein polynomialer Ansatz der Form

$y(x) \approx \alpha_0 + \alpha_1 x + \alpha_2 x^2 + \alpha_3 x^3 + \alpha_4 x^4$

versucht. Als Lösung ergibt sich das Polynom 4. Grades

$y(x) \approx 47{,}86 + 2{,}2 x -0{,}04809 x^2 + 0{,}0004935 x^3 -0{,}000002148 x^4$.

Die Messpunkte weichen im Mittel (Standardabweichung) 0,19 kg von der Modellfunktion ab. Reduziert man den Grad des Polynoms auf 3, erhält man die Lösung

$y(x) \approx 54{,}22 + 1{,}515 x -0{,}0226 x^2 + 0{,}0001002 x^3$

mit einer mittleren Abweichung von 0,22 kg und beim Polynomgrad 2 die Lösung

$y(x) \approx 61{,}42 + 0{,}9397 x -0{,}008881 x^2$

mit einer mittleren Abweichung von 0,42 kg. Wie zu erkennen ist, ändern sich beim Wegfallen der höheren Terme die Koeffizienten der niedrigeren Terme. Die Methode versucht, das Beste aus jeder Situation herauszuholen. Entsprechend werden die fehlenden höheren Terme mit Hilfe der niedrigeren Terme so gut wie möglich ausgeglichen, bis das mathematische Optimum erreicht ist. Mit dem Polynom zweiten Grades (Parabel) wird der Verlauf der Messpunkte noch sehr gut beschrieben (siehe Abbildung).

Spezialfall einer linearen Ausgleichsfunktion mit mehreren Variablen

Ist die Modellfunktion ein mehrdimensionales Polynom erster Ordnung, besitzt also statt nur einer Variablen x mehrere unabhängige Modellvariablen $x_1, \ldots, x_N$, erhält man eine lineare Funktion der Form

$f(x_1,\ldots,x_N;\alpha_0, \alpha_1, \dots, \alpha_N )= \alpha_0 + \alpha_1 x_1 + \cdots + \alpha_N x_N$

die auf die Residuen

$\begin{matrix} r_1 = \alpha_0 + \alpha_1 x_{1,1} + \cdots + \alpha_j x_{j,1}+ \cdots +\alpha_N x_{N,1} - y_1\\ r_2 = \alpha_0 + \alpha_1 x_{1,2} + \cdots + \alpha_j x_{j,2}+ \cdots +\alpha_N x_{N,2} - y_2\\ \vdots \\ r_i = \alpha_0 + \alpha_1 x_{1,i} + \cdots + \alpha_j x_{j,i}+ \cdots +\alpha_N x_{N,i} - y_i\\ \vdots\\ r_n = \alpha_0 + \alpha_1 x_{1,n} + \cdots + \alpha_j x_{j,n}+ \cdots +\alpha_N x_{N,n} - y_n\\ \end{matrix}$

führt und über den Minimierungsansatz

$\min_{\alpha} \sum_{i=1}^n r_i^2$

gelöst werden kann.

Der allgemeine lineare Fall

Zweidimensionale Polynomfläche zweiter Ordnung mit 3 × 3 = 9 Basisfunktionen:
f(x₁, x₂) = α₀ + α₁x₁¹ + α₂x₁² + α₃x₂¹ + α₄x₁¹x₂¹ + α₅x₁²x₂¹ + α₆x₂² + α₇x₁¹x₂² + α₈x₁²x₂²

Im Folgenden soll der allgemeine Fall von beliebigen linearen Modellfunktionen mit beliebiger Dimension gezeigt werden. Zu einer gegebenen N-dimensionalen Messwertfunktion

$y(x_1, x_2, \dots, x_N)$

mit N unabhängigen Variablen sei eine optimal angepasste lineare Modellfunktion

$f(x_1, x_2, \dots, x_N; \alpha_1, \alpha_2, \dots, \alpha_m) = \sum_{j=1}^{m}\alpha_j\varphi_j(x_1, x_2, \dots, x_N)$

gesucht, deren quadratische Abweichung dazu minimal sein soll. x_i sind dabei die Funktionskoordinaten, α_j die zu bestimmenden linear eingehenden Parameter und φ_j beliebige zur Anpassung an das Problem gewählte linear unabhängige Funktionen.

Bei n gegebenen Messpunkten $(x_{1,1},x_{2,1},\dots ,x_{N,1};y_1)$, $(x_{1,2},x_{2,2},\dots ,x_{N,2};y_2)$,… , $(x_{1,n},x_{2,n},\dots ,x_{N,n};y_n)$ erhält man die Anpassungsfehler

$\begin{matrix} r_1 &=&\alpha_0 \varphi_0(x_{1,1},\dots, x_{N,1}) & + \alpha_1 \varphi_1(x_{1,1},\dots, x_{N,1}) & + & \cdots &+ \alpha_m\varphi_m(x_{1,1},\dots, x_{N,1}) - y_1\\ r_2 &=&\alpha_0 \varphi_0(x_{1,2},\dots, x_{N,2}) & + \alpha_1 \varphi_1(x_{1,2},\dots, x_{N,2}) & + & \cdots &+ \alpha_m\varphi_m(x_{1,2},\dots, x_{N,2}) - y_2\\ &\vdots\\ r_i &=&\alpha_0 \varphi_0(x_{1,i},\dots, x_{N,i}) & + \alpha_1 \varphi_1(x_{1,i},\dots, x_{N,i}) & + & \cdots &+ \alpha_m\varphi_m(x_{1,i},\dots, x_{N,i}) - y_i\\ &\vdots\\ r_n &=&\alpha_0 \varphi_0(x_{1,n},\dots, x_{N,n}) & + \alpha_1 \varphi_1(x_{1,n},\dots, x_{N,n}) & + & \cdots &+ \alpha_m\varphi_m(x_{1,n},\dots, x_{N,n}) - y_n\\ \end{matrix}$

oder in Matrixschreibweise

r = Aα − y,

wobei der Vektor r die r_i zusammenfasst, die Matrix A die Basisfunktionswerte $A_{ij}:=\varphi_j(x_{1,i},\dots, x_{N,i})$, der Parametervektor α die Parameter α_j und der Vektor y die Beobachtungen y_i.

Das Minimierungsproblem kann dann in der Form

$\min_{\alpha} \sum_{i=1}^n r_i^2 = \min_\alpha\|f(\alpha)-y\|^2_2 =\min_{\alpha}\|A\alpha-y\|_2^2$.

geschrieben werden.

Lösung des Minimierungsproblems

Das Minimierungsproblem ergibt sich wie im allgemeinen linearen Fall gezeigt als

$\min_\alpha\|A\alpha-y\|_2^2 = \min_\alpha(A\alpha -y)^T(A\alpha-y) = \min_\alpha (\alpha^TA^TA\alpha - 2y^TA\alpha + y^Ty).$

Dieses Problem hat immer eine Lösung. Hat die Matrix A vollen Rang, so ist sie sogar eindeutig. Die partiellen Ableitungen bezüglich der α_j und Nullsetzen derselben zum Bestimmen eines Extremums ergeben ein lineares System von Normalgleichungen (auch Normalengleichungen)

A^TAα = A^Ty,

das die Lösung des Minimierungsproblems liefert und im Allgemeinen numerisch gelöst werden muss. Die Matrix A^TA ist positiv definit, so dass es sich beim gefundenen Extremum um ein Minimum handelt.^[7] Damit kann das Lösen des Minimierungsproblems der linearen Modellfunktionen auf das Lösen eines Gleichungssystems reduziert werden. Im einfachen Fall einer Ausgleichsgeraden kann dessen Lösung, wie gezeigt wurde, sogar direkt als einfache Formel angegeben werden.

Alternativ lassen sich die Normalgleichungen in der Darstellung

$A^TA\alpha-A^Ty= \begin{pmatrix} \left\langle \varphi_0,\varphi_0\right\rangle & \left\langle \varphi_0,\varphi_1\right\rangle & \cdots & \left\langle \varphi_0,\varphi_m\right\rangle \\ \left\langle \varphi_1,\varphi_0\right\rangle & \left\langle \varphi_1,\varphi_1\right\rangle & \cdots & \left\langle \varphi_1,\varphi_m\right\rangle \\ \vdots & \vdots & \ddots & \vdots \\ \left\langle \varphi_m,\varphi_0\right\rangle & \left\langle \varphi_m,\varphi_1\right\rangle & \cdots & \left\langle \varphi_m,\varphi_m\right\rangle \\ \end{pmatrix} \begin{pmatrix} \alpha_{0} \\ \alpha_{1} \\ \vdots \\ \alpha_{m} \end{pmatrix} - \begin{pmatrix} \left\langle y,\varphi_0\right\rangle \\ \left\langle y,\varphi_1\right\rangle \\ \vdots \\ \left\langle y,\varphi_m\right\rangle \\ \end{pmatrix} = 0.$

ausschreiben, wobei $\left\langle \cdot,\cdot\right\rangle$ das Skalarprodukt symbolisiert. Die Basisfunktionen φ_i sind als Vektoren $\vec{\varphi_i} =$$(\varphi_i(x_{1,1},\dots,x_{N,1}),$ $\varphi_i(x_{1,2},\dots,x_{N,2}),$, $\ldots, \varphi_i(x_{1,n},\dots,x_{N,n}))$ zu lesen mit den n diskreten Stützstellen am Ort der Beobachtungen $y = \vec y = (y_1, y_2, \ldots,y_n)$.

Ferner lässt sich das Minimierungsproblem mit einer Singulärwertzerlegung gut analysieren. Diese motivierte auch den Ausdruck der Pseudoinversen, einer Verallgemeinerung der normalen Inversen einer Matrix. Diese liefert dann eine Sichtweise auf nichtquadratische lineare Gleichungssysteme, die einen nicht stochastisch, sondern algebraisch motivierten Lösungsbegriff erlaubt.

Numerische Behandlung der Lösung

Zur numerischen Lösung des Problems gibt es zwei Wege. Zum Einen können die Normalgleichungen

A^TAα = A^Ty

gelöst werden, die eindeutig lösbar sind, falls die Matrix A vollen Rang hat. Ferner hat die Systemmatrix A^TA die Eigenschaft, positiv definit zu sein, ihre Eigenwerte sind also alle positiv. Zusammen mit der Symmetrie von A^TA kann dies beim Einsatz von numerischen Verfahren zur Lösung ausgenutzt werden: beispielsweise mit der Cholesky-Zerlegung oder dem CG-Verfahren. Da beide Methoden von der Kondition der Matrix stark beeinflusst werden, ist dies manchmal keine empfehlenswerte Herangehensweise: Ist schon A schlecht konditioniert, so ist A^TA quadratisch schlecht konditioniert. Dies führt dazu, dass Rundungsfehler so weit verstärkt werden können, dass sie das Ergebnis unbrauchbar machen.

Zum anderen liefert das ursprüngliche Minimierungsproblem eine stabilere Alternative, da es bei kleinem Wert des Minimums eine Kondition in der Größenordnung der Kondition von A, bei großen Werten des Quadrats der Kondition von A hat. Um die Lösung zu berechnen wird eine QR-Zerlegung verwendet, die mit Householdertransformationen oder Givens-Rotationen erzeugt wird. Grundidee ist, dass orthogonale Transformationen die euklidische Norm eines Vektors nicht verändern. Damit ist

$\|A\alpha-y\|_2 = \|Q(A\alpha-y)\|_2$

für jede orthogonale Matrix Q. Zur Lösung des Problems kann also eine QR-Zerlegung von A berechnet werden, wobei man die rechte Seite direkt mittransformiert. Dies führt auf eine Form

$\|R\alpha-Q^T y\|_2$

mit $R=\begin{pmatrix} \tilde{R} \\ 0 \end{pmatrix},$ wobei $\tilde{R} \in \mathbb{R}^{m\times m}$ eine rechte obere Dreiecksmatrix ist. Die Lösung des Problems ergibt sich somit durch die Lösung des Gleichungssystems

$\tilde{R}\begin{pmatrix} \alpha_1 \\ \vdots \\ \alpha_m\end{pmatrix}=\begin{pmatrix} (Q^T y)_1 \\ \vdots \\ (Q^T y)_m \end{pmatrix}.$

Die Norm des Minimums ergibt sich dann aus den restlichen Komponenten der transformierten rechten Seite (Qy)_{m + 1},...,(Qy)_n, da die dazugehörigen Gleichungen aufgrund der Nullzeilen in R nie erfüllt werden können.

In der statistischen Regressionsanalyse spricht man bei mehreren gegebenen Variablen x_j von multipler Regression. Der Ansatz ist auch als OLS (ordinary least squares) bekannt, im Gegensatz zu GLS (generalised least squares), dem verallgemeinerten Regressionsmodell bei Residuen, die von der Verteilungsannahme wie Unkorreliertheit und Homoskedastie abweichen. Dagegen liegen bei multivariater Regression für jede Beobachtung i (i = 1,...,n) r viele y-Werte vor, so dass statt eines Vektors eine $n \times r$-Matrix Y vorliegt. Die linearen Regressionsmodelle sind in der Statistik wahrscheinlichkeitstheoretisch intensiv erforscht worden. Besonders in der Ökonometrie werden beispielsweise komplexe rekursiv definierte lineare Strukturgleichungen analysiert, um volkswirtschaftliche Systeme zu modellieren.

Probleme mit Nebenbedingungen

Häufig sind Zusatzinformationen an die Parameter bekannt, die durch Nebenbedingungen formuliert werden, die dann in Gleichungs- oder Ungleichungsform vorliegen. Gleichungen tauchen beispielsweise auf, wenn bestimmte Datenpunkte interpoliert werden sollen. Ungleichungen tauchen häufiger auf, in der Regel in der Form von Intervallen für einzelne Parameter. Im Einführungsbeispiel wurde die Federkonstante erwähnt, diese ist immer größer Null und kann für den konkret betrachteten Fall immer nach oben abgeschätzt werden.

Im Gleichungsfall können diese bei einem sinnvoll gestellten Problem genutzt werden, um das ursprüngliche Minimierungsproblem in eines niedrigerer Dimension umzuformen, dessen Lösung die Nebenbedingungen automatisch erfüllt.

Schwieriger ist der Ungleichungsfall. Hier ergibt sich bei linearen Ungleichungen das Problem

$\min_\alpha \|\vec{f}-\vec{y}\|_2$ mit $l \leq C\alpha\leq u$, $C \in \mathbb{R}^{n\times n},$

wobei die Ungleichungen komponentenweise gemeint sind. Dieses Problem ist als konvexes Optimierungsproblem eindeutig lösbar und kann beispielsweise mit Methoden zur Lösung solcher angegangen werden.

Quadratische Ungleichungen ergeben sich beispielsweise bei der Nutzung einer Tychonow-Regularisierung zur Lösung von Integralgleichungen. Die Lösbarkeit ist hier nicht immer gegeben. Die numerische Lösung kann beispielsweise mit speziellen QR-Zerlegungen erfolgen.

Nichtlineare Modellfunktionen

Mit dem Aufkommen leistungsfähiger Rechner gewinnt insbesondere die nichtlineare Regression an Bedeutung. Hierbei gehen die Parameter nichtlinear in die Funktion ein. Nichtlineare Modellierung ermöglicht im Prinzip die Anpassung von Daten an jede Gleichung der Form y = f(α). Da diese Gleichungen Kurven definieren, werden die Begriffe nichtlineare Regression und „curve fitting“ zumeist synonym gebraucht.

Manche nichtlineare Probleme lassen sich durch geeignete Substitution in lineare überführen und sich dann wie oben lösen. Ein multiplikatives Modell von der Form

$y = \alpha_0 \cdot x^{\alpha_1}$

lässt sich beispielsweise durch Logarithmieren in ein additives System überführen. Dessen Parameter können dann berechnet werden. Dieser Ansatz findet unter Anderem in der Wachstumstheorie Anwendung.

Im Allgemeinen ergibt sich bei nichtlinearen Modellfunktionen ein Problem der Form

$\min_\alpha\|f(\alpha)-y\|_2,$

mit einer nichtlinearen Funktion f. Partielle Differentiation ergibt dann ein System von Normalgleichungen, das nicht mehr analytisch gelöst werden kann. Eine numerische Lösung kann hier iterativ mit dem Gauß-Newton-Verfahren erfolgen. Jenes hat allerdings das Problem, dass die Konvergenz des Verfahrens nicht gesichert ist.

Aktuelle Programme arbeiten häufig mit einer Variante, dem Levenberg-Marquardt-Algorithmus. Bei diesem Verfahren ist zwar die Konvergenz ebenfalls nicht gesichert, jedoch wird durch eine Regularisierung die Monotonie der Näherungsfolge garantiert. Zudem ist das Verfahren bei größerer Abweichung der Schätzwerte toleranter als die Ursprungsmethode. Beide Verfahren sind mit dem Newton-Verfahren verwandt und konvergieren meist quadratisch, in jedem Schritt verdoppelt sich also die Zahl der korrekten Nachkommastellen.

Wenn die Differentiation auf Grund der Komplexität der Zielfunktion zu aufwändig ist, stehen eine Reihe anderer Verfahren als Ausweichlösung zu Verfügung, die keine Ableitungen benötigen, siehe bei Methoden der lokalen nichtlinearen Optimierung.

Beispiel aus der Enzymkinetik einer nicht linearisierbaren Modellfunktion

Ein Beispiel für Regressionsmodelle, die voll nichtlinear sind, ist die Enzymkinetik. Hier ist zu fordern, dass nur y (Reaktionsgeschwindigkeit) und nicht α (Substratkonzentration) einem Fehler unterliegt und damit α als Variable genutzt werden kann. Die Lineweaver-Burk-Beziehung ist zwar eine algebraisch korrekte Umformung der Michaelis-Menten-Gleichung v = Vmax × [S] / (Km + [S]), ihre Anwendung liefert aber nur korrekte Ergebnisse, wenn die Messwerte fehlerfrei sind. Dies ergibt sich aus der Tatsache, dass sich die Realität nur mit einer erweiterten Michaelis-Menten-Beziehung

$\nu_i = \frac{V_\max\left[S_i\right]}{K_m+\left[S_i\right]}(1+e_i)\ \boldsymbol{\nu}_i$

mit e_i als Fehlerparameter, beschreiben lässt. Diese Gleichung lässt sich nicht mehr linearisieren, also muss hier die Lösung iterativ ermittelt werden.

Anforderungen an die Daten

Die Normalverteilungsannahme für die abhängige Variable y ist nicht zwingend notwendig. Es sollen lediglich keine Ausreißer in den Messwerten vorliegen, da diese das Schätzergebnis verzerren. Außerdem ist Multikollinearität zwischen den zu schätzenden Parametern ungünstig, da diese numerische Probleme verursacht. Im übrigen können auch Regressoren, die weit von den anderen entfernt liegen, die Ergebnisse der Ausgleichsrechnung stark beeinflussen. Man spricht hier von Werten mit großer Hebelkraft (High Leverage Value).

Multikollinearität

Multikollinearität entsteht, wenn die Messreihen zweier gegebener Variablen x_i und x_j sehr hoch korreliert sind, also fast linear abhängig sind. Im linearen Fall bedeutet dies, dass die Determinante der Normalgleichungsmatrix A^TA sehr klein und die Norm der Inversen umgekehrt sehr groß, die Kondition von A^TA ist also stark beeinträchtigt. Die Normalgleichungen sind dann numerisch schwer zu lösen. Die Lösungswerte können unplausibel groß werden und bereits kleine Änderungen in den Beobachtungen bewirken große Änderungen in den Schätzwerten.

Ausreißer

Ausreißer von y:
Der Wert zieht die Gerade nach oben

Als Ausreißer sind Datenwerte definiert, die „nicht in eine Messreihe passen“. Diese Werte beeinflussen die Berechnung der Parameter stark und verfälschen das Ergebnis. Um dies zu vermeiden, müssen die Daten auf fehlerhafte Beobachtungen untersucht werden. Die entdeckten Ausreißer können beispielsweise aus der Messreihe ausgeschieden werden oder es sind alternative ausreißerresistente Berechnungsverfahren wie gewichtete Regression oder das Drei-Gruppen-Verfahren anzuwenden.

Im ersten Fall wird nach der ersten Berechnung der Schätzwerte durch statistische Tests geprüft, ob Ausreißer in einzelnen Messwerten vorliegen. Diese Messwerte werden dann ausgeschieden und die Schätzwerte erneut berechnet. Dieses Verfahren eignet sich dann, wenn nur wenige Ausreißer vorliegen.

Bei der gewichteten Regression werden die abhängigen Variablen y in Abhängigkeit von ihren Residuen gewichtet. Ausreißer, d. h. Beobachtungen mit großen Residuen, erhalten ein geringes Gewicht, das je nach Größe des Residuums abgestuft sein kann. Beim Algorithmus nach Mosteller und Tukey (1977), der als „biweighting“ bezeichnet wird, werden unproblematische Werte mit 1 und Ausreißer mit 0 gewichtet, was die Unterdrückung des Ausreißers bedingt. Bei der gewichteten Regression sind in der Regel mehrere Iterationsschritte erforderlich, bis sich die Menge der erkannten Ausreißer nicht mehr ändert.

Verallgemeinerte Kleinste-Quadrate-Modelle

Hauptartikel: Regressionsanalyse.

Weicht man die starken Anforderungen im Verfahren an die Beobachtungsfehler auf, erhält man so genannte verallgemeinerte Kleinste-Quadrate-Ansätze. Wichtige Spezialfälle haben dann wieder eigene Namen, etwa die gewichteten kleinsten Quadrate, bei denen die Fehler zwar weiter als unkorreliert angenommen werden, aber nicht mehr von gleicher Varianz. Dies führt auf ein Problem der Form

$\|D(A\alpha-y)\|_2,$

wobei D eine Diagonalmatrix ist. Variieren die Varianzen stark, so haben die entsprechenden Normalgleichungen eine sehr große Kondition, weswegen das Problem direkt gelöst werden sollte.

Nimmt man noch weiter an, dass die Fehler in den Messdaten auch in der Modellfunktion berücksichtigt werden sollten, ergeben sich die „totalen kleinsten Quadrate“ in der Form

$\min_{E, r}\|(E, r)\|_F, (A+E)\alpha = b+r,$

wobei E der Fehler im Modell und r der Fehler in den Daten ist^[8][9].

Schließlich gibt es noch die Möglichkeit, keine Normalverteilung zugrunde zu legen. Dies entspricht beispielsweise der Minimierung nicht in der euklidischen Norm, sondern der 1-Norm. Solche Modelle sind Thema der Regressionsanalyse.

Literatur

• Åke Björck: Numerical Methods for Least Squares Problems. SIAM, Philadelphia 1996, ISBN 0-89871-360-9.
• Walter Großmann: Grundzüge der Ausgleichsrechnung. Springer Verlag, Berlin Heidelberg New York 1969 (3. erw. Aufl.), ISBN 3-540-04495-7.
• Richard J. Hanson, Charles L. Lawson: Solving least squares problems. SIAM, Philadelphia 1995, ISBN 0-89871-356-0.
• Frederick Mosteller, John W. Tukey: Data Analysis and Regression – a second course in statistics. Addison-Wesley, Reading MA 1977, ISBN 0-201-04854-X.
• Gerhard Opfer: Numerische Mathematik für Anfänger. Eine Einführung für Mathematiker, Ingenieure und Informatiker. Vieweg, Braunschweig 2002 (4. Aufl.), ISBN 3-528-37265-6.
• Peter Schönfeld: Methoden der Ökonometrie. 2 Bd. Vahlen, Berlin-Frankfurt 1969–1971.
• Eberhard Zeidler (Hrsg.): Taschenbuch der Mathematik. Begründet v. I.N. Bronstein, K.A. Semendjajew. Teubner, Stuttgart-Leipzig-Wiesbaden 2003, ISBN 3-8171-2005-2.
• T. Strutz: Data Fitting and Uncertainty (A practical introduction to weighted least squares and beyond). Vieweg+Teubner, ISBN 978-3-8348-1022-9.

Weblinks

 Wikibooks: Einführung in die Regressionsrechnung – Lern- und Lehrmaterialien

Einzelnachweise

1. ↑ Vgl. S. 436 von Moritz Cantor: Gauß: Karl Friedrich G.. In: Allgemeine Deutsche Biographie (ADB). Band 8, Duncker & Humblot, Leipzig 1878, S. 430–445.
2. ↑ Adrien-Marie Legendre: Nouvelles méthodes pour la détermination des orbites des comètes. Paris 1805, S. 72-80 (Anhang): Sur la Méthode des moindres quarrés.
3. ↑ Carl Friedrich Gauß: Theoria Motus Corporum Coelestium in sectionibus conicis solem ambientium. Göttingen 1809; Carl Haase (Übers.): Theorie der Bewegung der Himmelskörper, welche in Kegelschnitten die Sonne umlaufen. Hannover 1865.
4. ↑ http://www-history.mcs.st-andrews.ac.uk/HistTopics/Matrices_and_determinants.html
5. ↑ Carl Friedrich Gauß: Theoria combinationis observationum erroribus minimis obnoxiae. 2 Tle. Göttingen 1821-23. (Commentationes Societatis Regiae Scientiarum Gottingensis recentiores, classis mathematicae, Bd. 5.); Supplementum Theoria combinationis observationum erroribus minimis obnoxiae. Göttingen 1826/28. (Commentationes Societatis Regiae Scientiarum Gottingensis recentiores, classis mathematicae, Bd. 6.); Anton Börsch; Paul Simon (Hgg.): Abhandlungen zur Methode der kleinsten Quadrate von Carl Friedrich Gauss. In deutscher Sprache. Berlin 1887.
6. ↑ Pete Stewart, 21. Juni 1991: http://www.netlib.org/na-digest-html/91/v91n26.html#4
7. ↑ Hans R. Schwarz, Norbert Köckler: Numerische Mathematik. Teuber, 2009 (7.,überarb. Aufl.). doi:10.1007/978-3-8348-9282-9 – ISBN 978-3-8348-9282-9, Seite 141, Kapitel 3.6 (Gauß-Approximation), Satz 3.23.
8. ↑ Sabine Van Huffel, Joos Vandewalle: The Total Least Squares Problem: Computational Aspects and Analysis, SIAM Publications, Philadelphia PA, 1991. ISBN 978-0-898712-75-9
9. ↑ Martin Plesinger: The Total Least Squares Problem and Reduction of Data in AX ≈ B, Dissertation, TU Liberec und ICS Prague, 2008.

Dieser Artikel wurde in die Liste der lesenswerten Artikel aufgenommen.