Files
eth-summaries/semester3/numcs/parts/05_curve-fitting/02_non-linear/02_gauss-newton.tex

80 lines
3.0 KiB
TeX

\subsubsection{Gauss-Newton Verfahren}
Direkt das Newton-Verfahren auf ein Problem anzuwenden kann unmöglich oder schwer praktikabel sein.
Die Idee des Gauss-Newton Verfahrens ist es, die komplizierte Funktion $F(x)$ lokal durch eine lineare Funktion approximiert, also:
\begin{align*}
F(x) \approx F(y) + DF(y) (x - y) = F(y) + DF(y)x - DF(y)y
\end{align*}
Falls man $A := DF(y)$ und $b = DF(y)y - F(y)$ definiert, so erhält man ein lineares Ausgleichsproblem:
\rmvspace
\begin{align*}
\argmin{x \in \R^n} \frac{1}{2} ||F(x)||^2_2 \approx \argmin{x \in \R^n} \frac{1}{2} ||F(y) + DF(y) x||^2_2 = \argmin{x \in \R^n} \frac{1}{2} ||Ax - b||^2_2
\end{align*}
\drmvspace
wobei $y$ eine Näherung der Lösung $x$ ist.
Die Iterationsvorschrift ist gegeben durch:
\rmvspace
\begin{align*}
x^{(k + 1)} = x^{(k)} - s \smallhspace \text{ mit } s := \argmin{z \in \R^n} ||F(x^{(k)}) - DF(x^{(k)})z||^2_2
\end{align*}
\begin{code}{python}
import numpy as np
def gauss_newton(start_vec, Func, Jacobian, tolerance):
# Start vector has to be chosen intelligently
s = np.linalg.lstsq(Jacobian(start_vec), Func(start_vec))[0]
start_vec = start_vec - s
# now we perform the iteration
while np.linalg.norm(s) > tolerance * np.linalg.norm(start_vec):
# every time we update x by subtracting s, found with the least square method
s = np.linalg.lstsq(Jacobian(start_vec), Func(start_vec))[0]
start_vec = start_vec - s
return start_vec
\end{code}
Der Vorteil ist, dass die zweite Ableitung nicht benötigt wird, jedoch ist die Konvergenzordnung niedrieger ($p \leq 2$)
\newpage
\setLabelNumber{all}{3}
\inlineex Wir haben zwei Modellfunktionen, $F_1(t) = a_1 + b_1 e^{-c_1 t}$ and $F_2(t) = a_2 - b_2 e^{-c_2 t}$. ($F_1$ ist ein Heizvorgang, $F_2$ ist ein Abkühlvorgang).
Untenstehender code berechnet die Lösung des nichtlinearen Ausgleichsproblems
\begin{code}{python}
import numpy as np
t = np.arange(0, 30, 5); n = len(t)
curr_heating = np.array([24.34, 18.93, 17.09, 16.27, 15.97, 15.91])
curr_cooling = np.array([9.66, 18.8, 22.36, 24.07, 24.59, 24.91])
# define the functions that have to be minimized
F_1 = lambda a: a[0] + a[1] * np.exp(-a[2] * t) - curr_heating
F_2 = lambda a: a[0] - a[1] * np.exp(-a[2] * t) - curr_cooling
# define the corresponding Jacobi matrices
def J_1(a):
mat = np.zeros((n, 3))
for k in range(n):
mat[k, 0] = 1.0
mat[k, 1] = np.exp(-t[k] * a[2])
mat[k, 2] = -t[k] * a[1] * np.exp(-t[k] * a[2])
return mat
def J_2(a):
mat = np.zeros((n, 3))
for k in range(n):
mat[k, 0] = 1.0
mat[k, 1] = -np.exp(-t[k] * a[2])
mat[k, 2] = t[k] * a[1] * np.exp(-t[k] * a[2])
return mat
# guess starting vector
x_1 = np.array([10.0, 5.0, 0.0])
x_2 = np.array([30.0, 10.0, 0.0])
# use the Gauss-Newton algorithm declared above
a_1 = gauss_newton(x_1, F_1, J_1, tolerance=10e-6)
a_2 = gauss_newton(x_2, F_2, J_2, tolerance=10e-6)
print("Heating ", a_1)
print("Cooling ", a_2)
\end{code}