Die Regression (Ausgleichsrechnung, Methode der kleinsten Fehlerquadrate, ordinary least squares) ist ein Standardwerkzeug in sehr vielen Disziplinen und hat viele praktische Anwendungen (Datenanalyse, Statistik, Prognoseerstellung, Machine Learning, Beschreibung von Zusammenhängen, Parameterschätzung für Systemidentifikation etc.).
Aus mathematischer Sicht ist die Regression ein quadratisches (und deshalb relativ einfaches) Optimierungsproblem, das immer lösbar ist und unter einfachen Voraussetzungen eine eindeutige, globale Lösung hat. Geometrisch betrachtet ist die Regression die orthogonale Projektion eines Vektors auf eine lineare Hülle von Vektoren.
Geschichte: Die Methode der kleinsten Fehlerquadrate wurde von Carl Friedrich Gauß (1777-1855) und Adrien-Marie Legendre (1752-1833) unabhängig voneinander entwickelt. Legendre veröffentlichte seine Methode 1805 in seinem Buch “Nouvelles méthodes pour la détermination des orbites des comètes”. Gauß veröffentlichte seine Methode 1809 in seinem Buch “Theoria motus corporum coelestium in sectionibus conicis solem ambientium”. Link: Wikipedia
Beispiel Ausgleichsgerade: Wir beginnen mit dem klassischen Beispiel des Fits einer Geraden durch \(n\) Datenpunkte.
import numpy as np
import matplotlib.pyplot as plt
n = 13 # number of data points
t = np.linspace(1, 10, num = n) # time points measurements
noise = 0.4*np.random.randn(n) # noise in measurement values
y = 2 + 0.5*t + noise # measurement values: line + noise
plt.figure(figsize=(5, 3))
plt.plot(t, y, 'o')
plt.xlabel('t')
plt.ylabel('y')
plt.grid(True)
Ziel: Fit einer Geraden \(y(t) = d + kt\) durch die Datenpunkte \((t_i,y_i)\), sodass der “Fehler”, der noch zu definieren ist, klein ist. In anderen Worten: Finde die “optimalen” Werte für \(d\) und \(k\) aus den Datenpunkten.
Vorgehensweise: Wir formulieren das lineare Gleichungssystem \[ \begin{aligned} d + kt_1 & = y_1 \\ d + kt_2 & = y_2 \\ \vdots & = \vdots \\ d + kt_n & = y_n \end{aligned} \] für die zwei unbekannten Größen \(d\) und \(k\) als Matrixgleichung \(Ax = b\): \[ \begin{pmatrix} 1 & t_1 \\ 1 & t_2 \\ \vdots & \vdots \\ 1 & t_n \end{pmatrix} \begin{pmatrix} d \\ k \end{pmatrix} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix}. \] Das lineare Gleichungssystem ist typischer Weise nicht lösbar, da meist keine Gerade durch alle Datenpunkte legbar ist. Stattdessen soll der “Fehler” minimiert werden. Wir definieren als \(i\)-ten Einzelfehler \(e_i\) die \(y\)-Differenz zwischen dem Fit \(\hat{y}_i := d + kt_i\) und dem zugehörigen Datenwert \(y_i\), d. h. \(e_i := \hat{y}_i - y_i = d + kt_i - y_i\). Der Fehlervektor \(e\) besteht aus allen \(n\) Einzelfehlern. In Matrixnotation lässt sich der Fehlervektor aus dem Fitvektor \(\hat{y} = Ax\) und dem Datenvektor \(y\), den wir allgemein als \(b\) schreiben, durch \(e = \hat{y} - y = Ax - b\) berechnen. Wir verlangen, dass der Fehlervektor minimale Länge haben soll. Das heißt \(\lVert e \rVert = \sqrt{e_1^2 + e_2^2 + \ldots + e_n^2}\) soll minimal sein. Das ist äquivalent zur Forderung, dass \(\lVert e \rVert^2 = e_1^2 + e_2^2 + \ldots + e_n^2\) minimal sein soll, weil Quadrieren über den nicht-negativen Zahlen eine streng monotone Funktion ist. In Matrixnotation schreibt man \(\lVert e \rVert = \lVert Ax - b \rVert\). Diese Vorgehensweise heißt “Methode der kleinsten Fehlerquadrate”, andere Bezeichnungen sind “Regression” und “least squares”.
Formulierung des allgemeinen Optimierungsproblems: Das obige Optimierungsproblem lässt sich in Matrixform schreiben als \[ \min\; \lVert Ax - b \rVert. \] Umgekehrt wird jedes Optimierungsproblem, das sich als \(\min\; \lVert Ax - b \rVert\) schreiben lässt, auch als Regressionsproblem bezeichnet.
Hinweis: Die Probleme \(\min\; \lVert Ax - b \rVert\) und \(\min\; \lVert Ax - b \rVert^2\) sind äquivalent.
Lösung in Matrixform: Falls der Rang der Matrix \(A\) maximal ist, dann hat das Regressionsproblem die eindeutige Lösung \[ \hat{x} = (A^T A)^{-1}A^T b. \] Falls der Rang von \(A\) nicht maximal ist, dann gibt es unendlich viele Lösungen \(\hat{x}\), die alle denselben optimalen Fit \(A\hat{x}\) bilden. Der (bzw. ein) optimale(r) Vektor \(\hat{x}\) ist auch Lösung des Gleichungssystems \[ A^T A \hat{x} = A^T b. \] Die Gleichungen dieses Gleichungssystems heißen Normalgleichungen. Der optimale Fit ergibt sich als \(\hat{y} = A\hat{x}\).
Lösungsvarianten in Python:
np.linalg.solvenp.linalg.lstsqcol_of_ones = np.ones(n)
A = np.stack((col_of_ones, t), axis=1)
b = y.reshape(13,1)
# print(f"{A = }")
# print(f"{b = }")
# via Formel:
x_hat_1 = np.linalg.inv(A.T @ A) @ A.T @ b
print(f"x_hat_1 = \n{x_hat_1}")
# via Normalgleichungen
x_hat_2 = np.linalg.solve(A.T @ A, A.T @ b)
print(f"x_hat_2 = \n{x_hat_2}")
# via lstsq
x_hat_3 = np.linalg.lstsq(A, b, rcond=None)[0]
print(f"x_hat_3 = \n{x_hat_3}")
# Fit A*x_hat:
y_hat = A @ x_hat_1 # same as y_hat = x_hat_1[0] + x_hat_1[1]*t
plt.figure(figsize=(5,3))
plt.plot(t, y, 'o', label='data')
plt.plot(t, y_hat, 'o-', label='fit')
plt.xlabel('t')
plt.ylabel('y')
plt.legend(loc='best')
plt.grid(True)x_hat_1 =
[[2.7576247 ]
[0.39732196]]
x_hat_2 =
[[2.7576247 ]
[0.39732196]]
x_hat_3 =
[[2.7576247 ]
[0.39732196]]
Orthogonale Projektion: Wir betrachten das Matrixprodukt \(Ax\) als Linearkombination der \(m\) Spalten \(a_i\) von \(A\) mit den Komponenten von \(x\): \[ Ax = x_1 a_1 + x_2 a_2 + \ldots + x_m a_m \] Diese Linearkombination der Vektoren \(a_i\) im \(\mathbb{R}^n\) soll durch optimale Wahl der \(x_i\) dem Vektor \(b\in\mathbb{R}^n\) möglichst nahe kommen, sodass der Differenzvektor \(e = Ax-b\) eine minimale Länge hat.
Für den Fall \(n=3\) und \(m=2\) lässt sich die Situation wie in Abbildung Abbildung 1 darstellen:
Die Abbildung suggeriert, dass die optimale Linearkombination \(\hat{y} = A\hat{x}\) der orthogonalen Projektion von \(b\) auf die von den Spalten \(a_i\) aufgespannte Ebene entspricht.
Wir zeigen nun, dass diese Intuition in allen Dimensionen \(n\) und \(m\) stimmt, sofern der Rang der Matrix \(A\) maximal ist und deshalb \(A^T A\) invertierbar ist. Dazu überprüfen wir, ob der Fit \(A\hat{x}\) und sein Fehler \(e = A\hat{x} - b\) rechtwinklig sind, indem wir das innere Produkt der beiden berechnen: \[ \begin{aligned} (A\hat{x})^Te & = \hat{x}^T A^T(A\hat{x}-b) \\ & = \hat{x}^T A^T (A (A^T A)^{-1} A^T b - b) \\ & = \hat{x}^T A^T A (A^T A)^{-1} A^T b - \hat{x}^T A^T b \\ & = \hat{x}^T A^T b - \hat{x}^T A^T b = 0. \end{aligned} \] Der optimale Fit \(\hat{y} = A\hat{x}\) ist also orthogonal zum Fehlervektor \(e = A\hat{x} - b\). Es gilt allgemein, dass die Länge des Fehlers \(e\) genau dann minimal ist, wenn der Fehler zu der von den Spaltenvektoren der Matrix \(A\) aufgespannten Hyperebene orthogonal ist. Die orthogonale Projektion von \(b\) auf diese Hyperebene ist der optimale Fit \(\hat{y}= A\hat{x} = A (A^T A)^{-1} A^T b\), und die Matrix \(A (A^T A)^{-1} A^T\) die zugehörige Projektionsmatrix.
Bemerkungen:
Für ein LGS \(Ax = b\) liefert die Regressionslösung \(\hat{x}\) eine Lösung des LGS, falls es mindestens eine Lösung hat. Falls das LGS keine Lösung hat, liefert die Regressionslösung \(\hat{x}\) die Least-Squares-Approximation!
np.linalg.lstsq(A, b, rcond=None)[0] verwendet werden:
sp.linalg.null_space(A) liefert dafür eine Basis des Nullraums.np.linalg.solve(A, b) die zugehörige eindeutige Lösung. Auch np.linalg.lstsq und np.linalg.inv(A)@b könnten verwendet werden, aber np.linalg.solve ist effizienter.Regressionen werden sehr oft für Prognosen verwendet: Die Spalten der \(n\times m\) Matrix \(A\) bilden dabei die jeweils \(n\) historischen Werte der \(m\) erklärenden Größen. Nach dem Fit an die \(n\) historischen Werte \(b\) der zu erklärenden Größe stehen die optimalen Koeffizienten im Vektor \(\hat{x}\) zur Verfügung. Damit können neue Werte der zu erklärenden Größe berechnet werden, indem die zugehörigen neuen Werte der \(m\) erklärenden Größen mit den Koeffizienten des Vektors \(\hat{x}\) multipliziert werden.
t_new = np.array([11, 12, 13])
col_of_ones_new = np.ones(3)
A_new = np.stack((col_of_ones_new, t_new), axis=1)
y_new = A_new @ x_hat_1 # same as y_new = x_hat_1[0] + x_hat_1[1]*t_new
plt.figure(figsize=(5,3))
plt.plot(t, y, 'o', label='data')
plt.plot(t_new, y_new, 'o-', label='forecast')
plt.xlabel('t')
plt.ylabel('y')
plt.legend(loc='best')
plt.grid(True)
Verallgemeinerung der Ausgleichsgerade als Polynom 1. Ordnung auf z. B. Polynome 3. Ordnung mit gesuchten Koeffizienten \(c_0\), \(c_1\), \(c_2\) und \(c_3\): \[ \begin{aligned} c_0 + c_1 t_1 + c_2 t_1^2 + c_3 t_1^3 & = y_1 \\ c_0 + c_1 t_2 + c_2 t_2^2 + c_3 t_2^3 & = y_2 \\ \vdots & = \vdots \\ c_0 + c_1 t_n + c_2 t_n^2 + c_3 t_n^3 & = y_n \end{aligned} \] als Matrixgleichung \(Ax=b\): \[ \begin{pmatrix} 1 & t_1 & t_1^2 & t_1^3 \\ 1 & t_2 & t_2^2 & t_2^3 \\ \vdots & \vdots \\ 1 & t_n & t_n^2 & t_n^3 \end{pmatrix} \begin{pmatrix} c_0 \\ c_1 \\ c_2 \\c_3 \end{pmatrix} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} \] Die optimalen Koeffizienten sind Lösung des Regressionsproblems \(\min\; \lVert Ax - b \rVert\).
Beispiel: Das exponentielle Wachstumsmodell \[ n(t) = \alpha^{t - t_0} \] mit unbekannten Parametern \(\alpha\) und \(t_0\) kann durch Logarithmieren zu einem linearen Modell \[ \begin{aligned} \log(n(t)) & = (t - t_0) \log(\alpha) \\ \log(n(t)) & = -t_0 \log(\alpha) + \log(\alpha)t \end{aligned} \] transformiert werden. Das resultierende lineare Modell kann als Ausgleichsgerade mittels Regression behandelt werden. Dabei entspricht \(d=-t_0 \log(\alpha)\) und \(k=\log(\alpha)\).
Jede quadratische, skalare Zielfunktion \(g(x)\) mit \(x\in\mathbb{R}^n\) kann (bis auf einen für die Optimierung irrelevanten additiven Term) in der Regressionsform \(\lVert Ax - b \rVert^2\) geschrieben und somit mittels Regression behandelt werden.
Beispiel: \[ \begin{aligned} g(x_1, x_2) =& 13 x_1^2 - 30 x_1 x_2 + 18 x_2^2 + 14 x_1 - 12 x_2 + 10 \\ =& \left\lVert \begin{pmatrix} 2 & -3 \\ -3 & 3 \\ \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} - \begin{pmatrix} 1 \\ 3 \end{pmatrix} \right\rVert^2 \end{aligned} \]
Die zentrale Python-Funktion ist numpy.linalg.lstsq. Sie liefert einen Vektor zurück, der als ersten Eintrag \(\hat{x}\) enthält. Die weiteren Einträge sind für uns nicht relevant.
Für die Beispiele und Aufgaben verwenden wir folgende Python-Bibliotheken, siehe Kapitel Python Tutorial:
import numpy as np
import scipy as sp
import matplotlib.pyplot as pltNicht-quadratisches, inhomogenes LGS mit unendlich vielen Lösungen. Partikuläre Lösung mit lstsq:
A = np.array([[ 1,-3, 1],
[-2, 0, 5]])
b = np.array([[ 0],
[-7]])
Ab = np.hstack((A, b))
print("Rang von A =", np.linalg.matrix_rank(A))
print("Rang von Ab =", np.linalg.matrix_rank(Ab))
print("Anzahl der Variablen =", np.shape(A)[1])
print("Daher: unendlich viele Lösungen.")
# Eine Lösung erhält man mit lstsq:
x = np.linalg.lstsq(A, b, rcond=None)[0]
print("Eine Lösung ist x = \n{}".format(x))
# Check, dass x eine Lösung ist
print("Ax - b =\n", A@x - b)
# Alle Lösungen erhält man aus x + beliebige Linearkombination der Spalten von N:
r = np.linalg.matrix_rank(A)
N = sp.linalg.null_space(A)
# Check, dass die Spalten von N mit A multipliziert Null ergeben:
print("AN =\n", A@N)Rang von A = 2
Rang von Ab = 2
Anzahl der Variablen = 3
Daher: unendlich viele Lösungen.
Eine Lösung ist x =
[[ 0.56451613]
[-0.20322581]
[-1.17419355]]
Ax - b =
[[ 4.44089210e-16]
[-1.77635684e-15]]
AN =
[[3.33066907e-16]
[4.44089210e-16]]Nicht-quadratisches, inhomogenes LGS ohne Lösung. lstsq liefert eine Least-Squares-Approximation:
A = np.array([[ 1,-3],
[-2, 0],
[ 5, 1]])
b = np.array([[ 0],
[-7],
[ 1]])
Ab = np.hstack((A, b))
print("Rang von A =", np.linalg.matrix_rank(A))
print("Rang von Ab =", np.linalg.matrix_rank(Ab))
print("Daher: keine Lösung.")
# Achtung : Das mit lstsq berechnete x ist keine Lösung!
x = np.linalg.lstsq(A, b, rcond=None)[0]
print("lstsq liefert x = \n{}".format(x))
print("Ax - b =\n", A@x - b)Rang von A = 2
Rang von Ab = 3
Daher: keine Lösung.
lstsq liefert x =
[[ 0.63513514]
[-0.02702703]]
Ax - b =
[[0.71621622]
[5.72972973]
[2.14864865]]Das nicht-quadratische LGS mit eindeutiger Lösung. Diese kann mit lstsq berechnet werden:
A = np.array([[ 1,-3],
[-2, 0],
[ 5, 1]])
b = np.array([[ 4],
[-2],
[ 4]])
Ab = np.hstack((A, b))
rank_A = np.linalg.matrix_rank(A)
rank_Ab = np.linalg.matrix_rank(Ab)
(m, n) = np.shape(A)
print("Rang von A =", rank_A)
print("Rang von Ab =", rank_Ab)
print("Anzahl der Variablen =", n)
print("Daher: genau eine Lösung.")
x = np.linalg.lstsq(A, b, rcond=None)[0]
print("eindeutige Lösung ist x = \n{}".format(x))
# Check, dass x eine Lösung ist
print("Ax - b =\n", A@x - b)Rang von A = 2
Rang von Ab = 2
Anzahl der Variablen = 2
Daher: genau eine Lösung.
eindeutige Lösung ist x =
[[ 1.]
[-1.]]
Ax - b =
[[-3.10862447e-15]
[ 0.00000000e+00]
[ 8.88178420e-16]]Laden Sie mit dem Befehl die Datei Global_Temperatures_2024.csv in Python.
D = np.genfromtxt('daten/Global_Temperatures_2024.csv', delimiter=',')Die Werte stammen von der NASA-Homepage GISS Surface Temperature Analysis. Sie repräsentieren jährliche Mittel der Abweichungen der Oberflächentemperaturen vom Mittel der Jahre 1951 bis 1980. Die erste Spalte enthält die Jahre, die zweite die jährlichen Mittelwerte, die dritte die 5-jährlichen Mittelwerte.
Lösung:
D = np.genfromtxt('daten/Global_Temperatures_2024.csv', delimiter=',')
t = D[:,0]
T = D[:,1]
n = len(t)
plt.figure(figsize=(5, 3))
plt.plot(t, T, '.-', label='historische Daten')
plt.xlabel('Jahr')
plt.ylabel('Temperaturabweichung')
plt.legend(numpoints=1, loc='best')
plt.grid(True)
Linearer (Gerade) und quadratischer Fit via Least-Squares:
A = np.column_stack( (np.ones((n, 1)), t) )
x_linear = np.linalg.lstsq(A, T, rcond=None)[0]
y_fit_linear = A@x_linear
# quadratisch:
A = np.column_stack( (np.ones((n, 1)), t, t**2) )
x_quadr = np.linalg.lstsq(A, T, rcond=None)[0]
y_fit_quadr = A@x_quadr
plt.figure(figsize=(5, 3))
plt.plot(t, T, '.-', label='historische Daten')
plt.plot(t, y_fit_linear, '-r', label='linearer Fit')
plt.plot(t, y_fit_quadr , '-g', label='quadratischer Fit')
plt.xlabel('Jahr')
plt.ylabel('Temperaturabweichung')
plt.legend(loc='best', numpoints=1)
plt.grid(True)
Prognosen:
t_ = np.arange(2018, 2047)
n_ = len(t_)
# linear:
A_ = np.column_stack( (np.ones((n_, 1)), t_) )
y_prog_linear = A_ @ x_linear
# quadratisch:
A_ = np.column_stack( (np.ones((n_, 1)), t_, t_**2) )
y_prog_quadr = A_ @ x_quadr
plt.figure(figsize=(5, 3))
plt.plot(t, T, '.-', label='historische Daten')
plt.plot(t_, y_prog_linear, '-r', label='linear')
plt.plot(t_, y_prog_quadr , '-g', label='quadratisch')
plt.xlabel('Jahr')
plt.ylabel('Temperaturabweichung')
plt.legend(numpoints=1, loc = 'best')
plt.grid(True)
Das Breast Cancer Wisconsin (Diagnostic) Data Set enthält Daten von 569 Patientinnen, die Sie aus der Datei wdbc.csv mit dem folgenden Code laden.
data = np.genfromtxt('daten/wdbc.csv', delimiter=',')
X = data[:, 2:]
b = data[:, 1]Die 30 Spalten der Matrix \(X\) enthalten die Werte von 30 Features, die Aufschluss über die Gut- oder Bösartigkeit des Krebs jeder Patientin (=Zeile) geben können. Die Einträge des Vektors \(b\) sind 1 für Patientinnen mit bösartigem Krebs und 0 für Patientinnen mit gutartigem Krebs.
n = len(b)
col_of_ones = np.ones((n, 1))
A = np.hstack((col_of_ones, X))
x_hat = np.linalg.lstsq(A, b, rcond=None)[0]
b_hat_raw = A@x_hat # raw predicted values
b_hat = np.round(b_hat_raw) # round to 0 or 1
plt.figure(figsize=(8, 1))
plt.plot(b, 'o') # ture values
plt.plot(b_hat, '+') # predicted values
plt.ylim(-0.5, 1.5)
plt.grid(True)
# accuray (=the fraction of correctly classified samples)
accuracy = sum(b_hat == b)/n
print(f"{accuracy = :.2}")accuracy = 0.96
In der Praxis wird aus einem Trainingsdatensatz \(A_{\text{train}}\) und \(b_{\text{train}}\), der nur einen Teil der Patientinnendaten enthält, ein Modell \(A_{\text{train}} \hat{x} \sim b_{\text{train}}\) gefittet (=gelernt). Dieses Modell wird verwendet, um aus den restlichen Daten \(A_{\text{test}}\) Prognosen \(A_{\text{test}} \hat{x}\) von \(b_{\text{test}}\) zu erstellen, die dann mit den echten Werten \(b_{\text{test}}\) verglichen werden. Wieso tut man das?
Ergebnisse:
Laden Sie mit folgenden Python-Befehlen die Daten der Datei Polyfit.csv in ein array D, dessen Spalten als Daten-Vektoren \(t\) und \(y\) weiterverwendet werden:
D = np.genfromtxt('daten/Polyfit.csv', delimiter=';')
y = D[:, 1]
t = D[:, 0]lstsq), um ein Polynom vom Grad 4 durch die Punktwolke zu fitten. Erzeugen Sie ein Stamm-Blatt-Diagramm der Koeffizienten. Warum sind die ungeraden Koeffizienten null?numpy.polyfit und numpy.poly1d um die gleichen Resultate zu erzielen.Mehr Information zum Beispiel unter Wikipedia: Mooresches Gesetz.
D = np.genfromtxt('daten/Moores_Law.csv', delimiter=';')
t = D[:, 0]
n = D[:, 1]und plotten Sie \(n\) gegen \(t\) mit den matplotlib Befehlen plot und semilogy. 2. Erklären Sie, wie die Methode der kleinsten Fehlerquadrate (Regression) verwendet werden kann, um \(\alpha\) und \(t_0\) zu bestimmen, so dass \[
n_i \approx \alpha^{t_i -t_0} \text{ for } i= 1, . . . , 13
\] gilt. 3. Lösen Sie das Least-Squares-Problem (d.h. die Regressionsaufgabe) in Python mit dem Befehl lstsq. Vergleichen Sie Ihre Resultate mit dem Mooreschen Gesetz, das behauptet, dass sich die Anzahl der Transistoren in Mikroprozessoren alle 2 Jahre verdoppelt.
Die Müllmenge in Millionen Tonnen pro Tag des Landes Lower Slobbovia von 1960 bis 1995 ist in folgender Tabelle wiedergegeben.
| Jahr \(t\) | 1960 | 1965 | 1970 | 1975 | 1980 | 1985 | 1990 | 1995 |
|---|---|---|---|---|---|---|---|---|
| Menge \(y\) | 86.0 | 99.8 | 135.8 | 155.0 | 192.6 | 243.1 | 316.3 | 469.5 |
Sie haben in einer Messreihe die Spannung \(U\) und den Strom \(I\) an einem Ohmschen Widerstand gemessen und die folgenden Daten erhalten:
| \(I\) (A) | 0.11 | 0.15 | 0.27 | 0.40 | 0.50 | 0.54 | 0.63 | 0.89 | 0.99 | 1.10 |
|---|---|---|---|---|---|---|---|---|---|---|
| \(U\) (V) | 1.10 | 2.08 | 2.94 | 4.39 | 4.87 | 5.99 | 7.70 | 8.30 | 8.50 | 10.03 |