Ausgleichsverfahren - Numerische Mathematik (i282)

Ausgleichsverfahren erzeugen ein Modell für den Zusammenhang zwischen gemessenen Größen. Obwohl nur endlich viele punktuelle Messungen vorliegen, beschreibt das Modell den Zusammenhang lückenlos. Messfehler und Widersprüche in den Messwerten werden bei der Erzeugung des Modells ausgeglichen.

Das wichtigeste Ausgleichsverfahren ist die Methode der kleinsten Quadrate (siehe unten). Ausgleichsverfahren kommen immer dann zum Einsatz, wenn die darstellung eines funktionalen Zusammenhangs durch eine endliche Anzahl von Punkte nicht genügt, z.B. bei

der Auswertung des funktionalen Zusammenhangs an beliebigen Zwischenstellen,
der Berechnung von Flächeninhalten unter dem Funktionsgraph (numerische Integration, siehe späteres Kapitel),
der Visualisierung von funktionalen Zusammenhängen.

Im Kontext des maschinellen Lernens werden Ausgleichsverfahren als Regressionsverfahren bezeichnet.

7.1Zu lösendes Problem¶

Der Zusammenhang zwischen zwei Größen $x$ und $y$ soll aus $n$ Messungen $(x_1,y_1),\ldots,(x_n,y_n)$ ermittelt werden. Dabei sind die Messstellen $x_k$ vorgegeben oder praktisch fehlerfrei gemessen. Die von den $x_k$ abhängenden Werte $y_k$ können Messfehler aufweisen.

Im Allgemeinen wird man deutlich mehr Messwerte zur Verfügung haben als (bei fehlerfreier Messung) eigentlich benötigt.

Beispiele:

$x_k$ sind Zeitpunkte zu denen Temperaturen $y_k$ gemessen werden.
$x_k$ sind Positionen entlang einer Strecke an denen Höhen $y_k$ gemessen werden.

import matplotlib.pyplot as plt
import numpy as np

x_meas = np.array([1, 2, 5, 7, 9])
y_meas = np.array([2.3, 6.4, 5.1, 12.8, 11.9])

fig, ax = plt.subplots()

ax.plot(x_meas, y_meas, 'ob', label='Messungen $(x_k,y_k)$')

ax.set_xlim(0, 10)
ax.set_xticks(x_meas)
ax.set_xticklabels(['$x_1$', '$x_2$', '$x_3$', '$x_4$', '$x_5$'])
ax.set_xlabel('Messstellen')

ax.set_ylim(0, 15)
ax.set_yticks(y_meas)
ax.set_yticklabels(['$y_1$', '$y_2$', '$y_3$', '$y_4$', '$y_5$'])
ax.set_ylabel('Messwerte')

ax.grid()
ax.legend(loc='upper left')

plt.show()

Aus dem Kontext der Messaufgabe ist der grundlegende funktionale Zusammenhang $f:\bbR\to\bbR$ zwischen den Größen $x$ und $y$ bekannt. Es müssen nur noch einige Parameter bestimmt werden, sodass $f$ die Messwerte gut repräsentiert. Beachte: Die Funktion $f$ kann auch nur auf einem Intervall von Interesse sein, was am Verfahren aber nichts ändert.

Beispiele:

Die Temperatur steigt linear mit der Zeit (Gerade, 2 Parameter).
Das Höhenprofil hat keine starken Gefälle oder Steigungen und kann durch eine stückweise lineare Funktion mit 20 Stützstellen hinreichend genau beschrieben werden (20 Parameter).

fig, ax = plt.subplots()

ax.plot(x_meas, y_meas, 'ob', label='Messungen $(x_k,y_k)$')

a = 1.171875
b = 2.075
ax.plot([0, 10], [a * 0 + b, a * 10 + b], '-r', label='$f$')

ax.set_xlim(0, 10)
ax.set_xticks(x_meas)
ax.set_xticklabels(['$x_1$', '$x_2$', '$x_3$', '$x_4$', '$x_5$'])
ax.set_xlabel('Messstellen')

ax.set_ylim(0, 15)
ax.set_yticks(y_meas)
ax.set_yticklabels(['$y_1$', '$y_2$', '$y_3$', '$y_4$', '$y_5$'])
ax.set_ylabel('Messwerte')

ax.grid()
ax.legend(loc='upper left')

plt.show()

7.2Wie Ansatzfunktion wählen?¶

Das hängt vom Kontext ab! Oft ergibt sich die Ansatzfunktion aus physikalischen Gesetzen oder aus der verfolgten Zielstellung heraus.

Übliche Ansätze für $f$ :

Gerade $f(x)={\color{red}a}\,x+{\color{red}b}$
Parabel $f(x)={\color{red}a}\,x^2+{\color{red}b}\,x+{\color{red}c}$
Polynom höheren Grades
stückweise linear Funktion mit $m$ gleichmäßig verteilten Stützstellen und Werten ${\color{red}a_1},\ldots,{\color{red}a_m}$ an diesen Stellen
Summen von verschobenen Hütchenfunktionen:
$f(x)=\sum_{l=1}^m{\color{red}a_l}\,h_l(x),$
wobei $h_1,\ldots,h_m$ verschobene Kopien einer Funktion $h$ sind

def hat(x):
    return 0.5 * np.maximum(0, 2 - np.abs(x))

fig, ax = plt.subplots()

# Hütchen plotten
x_grid = np.array([2, 4, 6, 8])
for k, x in enumerate(x_grid):
    ax.plot(x_grid, hat(x_grid - x), label=f'$h_{int(k+1)}$', lw=1)

# Koeffizienten vor verschobenen Hütchen
a = [1.5, 0.5, 2, 2.5]

# stückweise lineare Funktion als Summe von verschobenen Hütchen bauen
y = np.zeros(len(a))
for k in range(len(a)):
    y += a[k] * hat(x_grid - x_grid[k])

# Funktion plotten (mit Formel in Legende)
label = '+'.join([f'{a[k]}\\,h_{k+1}(x)' for k in range(len(a))])
ax.plot(x_grid, y, '-r', label=f'$f(x)={label}$', lw=3)

ax.set_xlim(x_grid[0], x_grid[-1])
ax.set_xticks(x_grid)
ax.set_xlabel('x')

ax.set_ylim(0, 4)
ax.set_yticks(a)
ax.set_yticklabels(['$a_1$', '$a_2$', '$a_3$', '$a_4$'])
ax.set_ylabel('y')

ax.grid()
ax.legend(loc='upper left')

plt.show()

7.3Was heißt „gut“?¶

Es gibt viele verschiedene richtige Antworten auf diese Frage!

Grundsätzlich möchte man die Differenzen $f(x_k)-y_k$ zwischen $f$ an den Messstellen und den zugehörigen Messwerten für alle $k$ klein bekommen.

fig, ax = plt.subplots()

a = 1.171875
b = 2.075

# Differenzen
for x, y in zip(x_meas, y_meas):
    ax.plot([x, x], [a * x + b, y], '-y', lw=3)
# letzte Differenz doppelt, damit ein (und nur ein!) Legendeneintrag erzeugt wird
ax.plot([x, x], [a * x + b, y], '-y', label='$|f(x_k)-y_k|$', lw=3)

# Messpunkte
ax.plot(x_meas, y_meas, 'ob', label='Messungen $(x_k,y_k)$')

# Gerade und Punkte darauf
ax.plot(x_meas, a * x_meas + b, 'or', label='Punkte $(x_k,f(x_k))$')
ax.plot([0, 10], [a * 0 + b, a * 10 + b], '-r', label='$f$')

ax.set_xlim(0, 10)
ax.set_xticks(x_meas)
ax.set_xticklabels(['$x_1$', '$x_2$', '$x_3$', '$x_4$', '$x_5$'])
ax.set_xlabel('Messstellen')

ax.set_ylim(0, 15)
ax.set_yticks(y_meas)
ax.set_yticklabels(['$y_1$', '$y_2$', '$y_3$', '$y_4$', '$y_5$'])
ax.set_ylabel('Messwerte')

ax.grid()
ax.legend(loc='upper left')

plt.show()

7.3.1Variante 1 (Summe aller Fehler minimieren)¶

Minimiere

\sum_{k=1}^n|{\color{red}f}(x_k)-y_k|

bezüglich der Parameter in $f$ .

Ist $f$ beispielsweise eine Gerade, so müssen wir

\sum_{k=1}^n|{\color{red}a}\,x_k+{\color{red}b}-y_k|

bezüglich $a$ und $b$ minimieren.

Probleme:

Funktion nicht differenzierbar (aber es gibt numerische Verfahren dafür, siehe Subgradient method)
Lösung oft nicht eindeutig
kleine Änderungen in den Messpunkten können zu großen Änderungen in der Lösung führen (schlechte Kondition)

Für Details siehe Least absolute deviations.

7.3.2Variante 2 (maximalen Fehler minimieren)¶

Minimiere

\max_{k=1,\ldots,n}|{\color{red}f}(x_k)-y_k|

bezüglich der Parameter in $f$ .

Probleme:

Funktion nicht differenzierbar (und es gibt keine wirklich guten numerischen Verfahren dafür)

7.3.3Variante 3 (Summe aller Fehlerquadrate minimieren)¶

„Rund-Lutschen“ von Variante 1 liefert Minimierung von

\sum_{k=1}^n\bigl({\color{red}f}(x_k)-y_k\bigr)^2

bezüglich der Parameter in $f$ .

Das ist differenzierbar und hat praktisch immer einen eindeutigen Minimierer!

Für die Eindeutigkeit werden mehr Messstellen benötigt als $f$ Parameter hat.

Diese Variante zum Lösen des Ausgleichsproblems ist als Methode der kleinesten Quadrate bekannt. Das Verfahren wurde von Carl Friedrich Gauß für Anwendungen in der Astronomie entwickelt und von ihm auch bei der Vermessung des Königreichs Hannover angewendet. Heute ist die Methode der kleinsten Quadrate in praktisch allen Wissenschaftszweigen im Einsatz, in denen Erkenntnisse aus Messungen gewonnen werden.

7.4Warum Quadrate?¶

Warum nicht dritte oder vierte Potenzen?

Einfache Gründe:

Die Ableitungen werden schön einfach.
$p=2$ ist die kleinste Potenz, sodass $t\mapsto t^p$ zweimal differenzierbar ist (hinreichende Optimalitätsbedingung!).
Je höher die Potenz, desto weniger Einfluss haben kleine Messfehler auf die Zielfunktion. Die gesuchte Funktion $f$ orientiert sich bei großen Potenzen im wesentlichen an Ausreißern in den Messungen.

Tiefer liegende Gründe:

Betrachtet man die Abweichungen nicht einzeln an jeder Stelle $x_k$ , sondern interpretiert sie als Abweichung zwischen den beiden Vektoren $[f(x_1),\ldots,f(x_n)]^{\mathrm{T}}$ und $[y_1,\ldots,y_n]^{\mathrm{T}}$ , so ist die Zielfunktion mit Potenz 2 gerade der übliche Abstand zwischen Vektoren im $\mathbb{R}^n$ .
Man kann zeigen, dass die Kleineste-Quadrate-Lösung in gewissem Sinn die bestmögliche Lösung ist, wenn die Messfehler gewisse statistische Eigenschaften haben (siehe Satz von Gauß-Markow für Details).

7.5Lösung des Minimierungsproblems¶

Lösen das Kleinste-Quadrate-Problem hier nur beispielhaft für eine Gerade $f(x)=a\,x+b$ . Für andere Ansätze $f$ erfolgen die Rechenschritte völlig analog, solange $f$ linear von den Parametern abhängt (siehe unten für den nichtlinearen Fall).

Die allgemeine Lösung wird anschließend ohne Rechenweg angeben.

7.5.1Ausgleichsgerade¶

Gegeben sind Punkte $(x_k,y_k)$ für $k=1,\ldots,n$ . Gesucht ist der Minimierer von

g({\color{red}a},{\color{red}b}):=\sum_{k=1}^n\bigl({\color{red}a}\,x_k+{\color{red}b}-y_k\bigr)^2

bezüglich $a,b\in\mathbb{R}$ .

Gradient von $g$ berechnen:

{\begin{align} \nabla g({\color{red}a},{\color{red}b}) &=\begin{bmatrix} \displaystyle 2\,\sum_{k=1}^n x_k\,\bigl({\color{red}a}\,x_k+{\color{red}b}-y_k\bigr)\\ \displaystyle 2\,\sum_{k=1}^n \bigl({\color{red}a}\,x_k+{\color{red}b}-y_k\bigr) \end{bmatrix}\\ &=\begin{bmatrix} \displaystyle 2\,{\color{red}a}\,\sum_{k=1}^n x_k^2+2\,{\color{red}b}\,\sum_{k=1}^n x_k-2\,\sum_{k=1}^n x_k\,y_k\\ \displaystyle 2\,{\color{red}a}\,\sum_{k=1}^n x_k+2\,{\color{red}b}\,\sum_{k=1}^n 1-2\,\sum_{k=1}^n y_k \end{bmatrix}\\ &=\begin{bmatrix} 2\,c_1\,{\color{red}a}+2\,c_2\,{\color{red}b}-2\,c_3\\ 2\,c_2\,{\color{red}a}+2\,n\,{\color{red}b}-2\,c_4 \end{bmatrix}, \end{align}}

(1)

wobei

c_1:=\sum_{k=1}^n x_k^2,

(2)

c_2:=\sum_{k=1}^n x_k,

(3)

c_3:=\sum_{k=1}^n x_k\,y_k,

(4)

c_4:=\sum_{k=1}^n y_k

(5)

konkrete Zahlen sind.

Der Gradient wird genau dann zum Nullvektor, wenn das lineare Gleichungssystem

\begin{bmatrix} c_1&c_2\\ c_2&n \end{bmatrix} \begin{bmatrix} {\color{red}a}\\ {\color{red}b} \end{bmatrix} =\begin{bmatrix} c_3\\ c_4 \end{bmatrix}

erfüllt ist.

Man kann zeigen, dass die Lösung des Gleichungssystems nicht nur ein stationärer Punkt, sondern ein globales Minimum von $g$ ist.

7.5.2Allgemeine Lösung¶

Hat die gesuchte Ausgleichsfunktion die Form

f(x)=\sum_{l=1}^m{\color{red}a_l}\,h_l(x),

wobei $h_1,\ldots,h_m$ beliebige Funktionen sind (nicht notwendig verschobene Hütchen), so soll die Funktion

g({\color{red}a_1},\ldots,{\color{red}a_m}):=\sum_{k=1}^n\left(\sum_{l=1}^m{\color{red}a_l}\,h_l(x_k)-y_k\right)^2

minimiert werden. Mit

\underline{a}:=\begin{bmatrix}a_1\\\vdots\\a_m\end{bmatrix},

(6)

H:=\begin{bmatrix}h_1(x_1)&\cdots&h_m(x_1)\\\vdots&\ddots&\vdots\\h_1(x_n)&\cdots&h_m(x_n)\end{bmatrix},

(7)

\underline{y}:=\begin{bmatrix}y_1\\\vdots\\y_n\end{bmatrix}

(8)

kann dies als

g({\color{red}\underline{a}})=\bigl\|H\,{\color{red}\underline{a}}-\underline{y}\bigr\|_2^2

geschrieben werden.

Für den Gradient erhalten wir

\nabla g({\color{red}\underline{a}})=2\,H^{\mathrm{T}}\,\bigl(H\,{\color{red}\underline{a}}-\underline{y}\bigr),

sodass die stationären Punkte Lösungen des linearen Gleichungssystems

\boxed{H^{\mathrm{T}}\,H\,{\color{red}\underline{a}}=H^{\mathrm{T}}\,\underline{y}}

sind. Dieses $(m\times m)$ -Gleichungssystem hat immer eine Lösung, die in allen praktisch relevanten Fällen ( $n\geq m$ ) auch eindeutig ist. Die Systemmatrix ist immer symmetrisch und in allen praktisch relevanten Fällen auch positiv definit. Über die Kondition des Gleichungssystems und damit über die Kondition des Kleinste-Quadrate-Verfahrens sind keine allgemeinen Aussagen möglich. Die Kondition wird maßgeblich vom Zusammenspiel der Messstellen mit den Ansatzfunktionen bestimmt. Bei üblichen Ansatzfunktionen und gleichmäßig verteilten Messstellen kann die Kondition als gut angenommen werden.

Beispiel

Wenden die allgemeine Formel für das Finden einer Ausgleichsgerade

f(x)=a_1\,x+a_2

an.

Haben $\quad m=2$ , $\quad h_1(x)=x$ , $\quad h_2(x)=1\quad$ und

\underline{a}:=\begin{bmatrix}a_1\\a_2\end{bmatrix},\qquad H=\begin{bmatrix} x_1&1\\ \vdots&\vdots\\ x_n&1 \end{bmatrix},\qquad \underline{y}:=\begin{bmatrix}y_1\\\vdots\\y_n\end{bmatrix}.

Systemmatrix und rechte Seite sind also

H^{\mathrm{T}}\,H=\begin{bmatrix} \displaystyle \sum_{k=1}^n x_k^2&\displaystyle \sum_{k=1}^n x_k\\ \displaystyle \sum_{k=1}^n x_k&n \end{bmatrix},\qquad H^{\mathrm{T}}\,\underline{y}=\begin{bmatrix} \displaystyle \sum_{k=1}^n x_k\,y_k\\ \displaystyle \sum_{k=1}^n y_k \end{bmatrix}.

7.6Hürden im praktischen Einsatz¶

7.6.1Verteilung der Messstellen¶

Die Messstellen sollten möglichst gleichmäßig im Definitionsbereich der gesuchten Funktion $f$ verteilt sein. In Bereichen mit überdurchschnittlich vielen Messstellen wird sich die Funktion stärker den Messpunkten anpassen als in Bereichen mit nur wenigen Messstellen.

Umgekehrt kann man diesen Effekt nutzen, um durch Vervielfachung von Messpunkten diesen Messpunkten mehr Einfluss auf die Lösung zu geben. Das kann sinnvoll sein, wenn gleichzeitig Messpunkte mit kleinem Messfehler und Messpunkte mit großen Messfehler vorliegen. Eine Vervielfachung der genaueren Punkte gibt diesen mehr Einfluss auf die Kleinste-Quadrate-Lösung.

7.6.2Nichtlineare Ansatzfunktionen¶

Das Lösen eines Kleinste-Quadrate-Problems als lineares Gleichungssystem funktioniert nur, wenn die Ansatzfunktion $f$ linear bezüglich ihrer Parameter ist.

Manche nichtlineare Ansatzfunktionen, wie z.B.

f(x)={\color{red}a}\,\mathrm{e}^{{\color{red}b}\,x},

können durch eine geeignete Transformation der Messwerte zu linearen Ansatzfunktionen umgebaut werden (IDVID 790).

Meistens, wie z.B. bei

f(x)=\sin({\color{red}a}\,x+{\color{red}b}),

müssen Lösungsverfahren für nichtlineare Optimierungsprobleme eingesetzt werden. Beispielsweise liefert $\nabla g(a_1,\ldots,a_m)=0$ für die Zielfunktion $g$ des Kleinste-Quadrate-Problems ein nichtlineares Gleichungssystem, welches mittels Newton-Verfahren gelöst werden kann.

7 Ausgleichsverfahren

7.1Zu lösendes Problem¶

7.2Wie Ansatzfunktion wählen?¶

7.3Was heißt „gut“?¶

7.3.1Variante 1 (Summe aller Fehler minimieren)¶

7.3.2Variante 2 (maximalen Fehler minimieren)¶

7.3.3Variante 3 (Summe aller Fehlerquadrate minimieren)¶

7.4Warum Quadrate?¶

7.5Lösung des Minimierungsproblems¶

7.5.1Ausgleichsgerade¶

7.5.2Allgemeine Lösung¶

7.6Hürden im praktischen Einsatz¶

7.6.1Verteilung der Messstellen¶

7.6.2Nichtlineare Ansatzfunktionen¶