eth-summaries/semester3/numcs/parts/03_zeros/08_quasi-newton.tex

\newsectionNoPB
\subsection{Quasi-Newton-Verfahren}
Falls $DF(x)$ zu teuer ist oder nicht zur Verfügung steht, können wir im Eindimensionalen das Sekantenverfahren verwenden.

Im höherdimensionalen Raum ist dies jedoch nicht direkt möglich und wir erhalten die Broyden-Quasi-Newton Methode:
\rmvspace
\begin{align*}
    J_{k + 1} := J_k + \frac{F(x^{(k + 1)}) (\Delta x^{(k)})^\top}{||\Delta x^{(k)}||_2^2}
\end{align*}

\drmvspace
Dabei ist $J_0$ z.B. durch $DF(x^{(0)})$ definiert.

\fancyremark{Broyden-Update} Das Broyden-Update ergibt bezüglich der $||\cdot||_2$-Norm die minimale korrektur der Jakobi-Matrix $J_k$ an, so dass die Sekantenbedingung erfüllt ist.
Die Implementierung erzielt man folgendermassen mit der \bi{Sherman-Morrison-Woodbury} Update-Formel:
\begin{align*}
    J_{k + 1}^{-1} = \left( I - \frac{J_k^{-1} F(x^{(k + 1)}) (\Delta x^{(k)})^\top}{||\Delta x^{(k)}||_2^2 + (\Delta x^{(k)})^\top J_k^{-1} F(x^{(k + 1)})} \right) J^{-1}_k
\end{align*}

Das Broyden-Quasi-Newton-Verfahren konvergiert langsamer als das Newton-Verfahren, aber schneller als das vereinfachte Newton-Verfahren. (\texttt{sp} ist \texttt{Scipy} und \texttt{np} logischerweise \texttt{Numpy} im untenstehenden code)

\begin{code}{python}
def fast_broyden(x0: np.ndarray, F, J, tol=1e-12, maxIter=20):
    x   = x0.copy()
    lup = lu_factor(J)
    
    s  = lu_solve(lup, F(x))
    sn = np.dot(s, s)
    x -= s
    
    # Book keeping, for Broyden Update
    dx     = np.zeros((maxIter, len(x)))
    dxn    = np.zeros(maxIter)
    dx[0]  = s
    dxn[0] = sn
    k = 1
    
    while sn > tol and k < maxIter:
        w = lu_solve(lup, F(x))     # Simplified Newton Update
        
        # Apply Broyden correction (Shermann-Morrison-Woodbury formula)
        for r in range(1, k):
            w += dx[r] * np.dot(dx[r-1], w) / dxn[r-1]
        z = np.dot(s, w)
        s = (1 + z/(sn-z)) * w
        x -= s      # Apply the iteration
        
        # Book keeping again
        sn = np.dot(s, s)
        dx[k] = s
        dxn[k] = sn
        k += 1
    
    return x, k
\end{code}