% vim: set sts=4 et :
\documentclass[11 pt,xcolor={dvipsnames,svgnames,x11names,table}]{beamer}
\usetheme[
bullet=circle, % Other option: square
bigpagenumber, % circled page number on lower right
topline=true, % colored bar at the top of the frame
shadow=false, % Shading for beamer blocks
watermark=BG_lower, % png file for the watermark
]{Flip}
\usepackage{bm}
\usepackage{tikz,pgfplots}
\usepackage{graphicx}
% SOME COMMANDS THAT I FIND HANDY
% \renewcommand{\tilde}{\widetilde} % dinky tildes look silly, dosn't work with fontspec
\newcommand{\comment}[1]{\textcolor{comment}{\footnotesize{#1}\normalsize}} % comment mild
\newcommand{\Comment}[1]{\textcolor{Comment}{\footnotesize{#1}\normalsize}} % comment bold
\newcommand{\COMMENT}[1]{\textcolor{COMMENT}{\footnotesize{#1}\normalsize}} % comment crazy bold
\newcommand{\Alert}[1]{\textcolor{Alert}{#1}} % louder alert
\newcommand{\ALERT}[1]{\textcolor{ALERT}{#1}} % loudest alert
%% "\alert" is already a beamer pre-defined
\newcommand*{\Scale}[2][4]{\scalebox{#1}{$#2$}}%
\graphicspath{{images/}} % Put all images in this directory. Avoids clutter.
% suppress frame numbering for backup slides
% you always need the appendix for this!
\newcommand{\backupbegin}{
\newcounter{framenumberappendix}
\setcounter{framenumberappendix}{\value{framenumber}}
}
\newcommand{\backupend}{
\addtocounter{framenumberappendix}{-\value{framenumber}}
\addtocounter{framenumber}{\value{framenumberappendix}}
}
\begin{document}
\tikzstyle{every picture}+=[remember picture]
{
\setbeamertemplate{sidebar right}{\llap{\includegraphics[width=\paperwidth,height=\paperheight]{bubble2}}}
\begin{frame}[c]%{\phantom{title page}}
\begin{center}
\begin{center}
\begin{columns}
\begin{column}{0.9\textwidth}
\flushright%\fontspec{Trebuchet MS}
\bfseries \Huge {Root Finding}
\end{column}
\begin{column}{0.2\textwidth}
%\includegraphics[width=\textwidth]{SHiP-2}
\end{column}
\end{columns}
\end{center}
\quad
\vspace{3em}
\begin{columns}
\begin{column}{0.6\textwidth}
\flushright \vspace{-1.8em} {%\fontspec{Trebuchet MS}
\large Marcin Chrzaszcz, Danny van Dyk\\\vspace{-0.1em}\small \href{mailto:mchrzasz@cern.ch}{mchrzasz@cern.ch}, \href{mailto:dany.van.dyk@gmail.com}{danny.van.dyk@gmail.com}}
\end{column}
\begin{column}{0.4\textwidth}
\includegraphics[height=1.3cm]{uzh-transp}
\end{column}
\end{columns}
\vspace{1em}
\vspace{0.5em}
\textcolor{normal text.fg!50!Comment}{Numerical Methods, \\ 26. September, 2016}
\end{center}
\end{frame}
}
\begin{frame}{Plan for today}
\begin{itemize}
\item \alert{General problem of finding the root of a function}\\
$\xi$ is a root of $f$ iff $f(\xi) = 0$.
\vfill
\item \alert{Begin with $D=1$}\\
What properties of $f(x)$ can we use to find one, any, or all roots
of $f$? What are the requirements on $f(x)$?
\vfill
\item \alert{How about $D>1$?}\\
Can we generalize root finding from $D=1$ do arbitrary $D$?
\end{itemize}
\end{frame}
\begin{frame}{General iterative procedure}
Let's assume that
\begin{itemize}
\item $f(x)$ is a our function of interest,
\item $\xi$ is the only root of $f$,
\item we have a point $x_0$ close to $\xi$
\end{itemize}
We now want to find a sequence $\lbrace x_0, x_1, \dots \rbrace$ that
\begin{enumerate}
\item converges toward $\xi$: $\lim_{k\to \infty} x_k = \xi$,
\item is iterative: $x_{k + 1} = \Phi[f](x_k)$
\end{enumerate}
We can attempt to Taylor expand $f$ around $x_0$ to order $N$, in order
to obtain the generator $\Phi[f]$ for the iteration:
\begin{block}{Expansion around $x_0$}
\begin{equation*}
f(\xi) = 0 = \sum_{n=0}^N \frac{(\xi - x_0)^n}{n!} f^{(n)}(x_0)
\end{equation*}
where $f^{(n)}(x_0)$ is the $n$th derivative of $f$ at the position $x_0$.
\end{block}
\end{frame}
\begin{frame}{Newton-Raphson method \hfill ($N=1$)}
\begin{block}{Expansion for $N=1$}
\begin{equation}
\tag{*}
f(\xi) = 0 = f(x_0) + (\xi - x_0) \cdot f'(x_0) + \mathcal{O}((\xi - x_0)^2)
\end{equation}
\end{block}
\vfill
Algorithm:
\begin{enumerate}
\item start with index $k=0$
\item solve equation (*), assuming a vanshing approximation error:
\begin{gather*}
x_{k + 1} \leftarrow x_k - \frac{f(x_k)}{f'(x_k)} \approx \xi + \mathcal{O}((\xi - x_k)^2)\\
\Phi[f](x) \equiv x - \frac{f(x)}{f'(x)}
\end{gather*}
\item if $f(x_{k + 1}) \leq t$, where $t$ is an a-prior threshold, then stop; otherwise
jump back to step \#2.
\end{enumerate}
\end{frame}
\begin{frame}{Illustration of Newton-Raphson}
\begin{columns}
\begin{column}[T]{.5\textwidth}
\begin{overlayarea}{\textwidth}{\textheight}
$\bm{k=0}$:
\only<2->{%
\begin{equation*}
x_0 := 1
\end{equation*}
}
\only<3->{%
$\bm{k=1}$:
\begin{align*}
x_0 & = 1 &
f(x_0) & = +0.54 \\
& &
f'(x_0) & = -0.84 \\
\end{align*}
}
\only<4->{%
\begin{equation*}
{\color{red} x_1 \leftarrow 1.64}
\end{equation*}
}
\only<6->{%
$\bm{k=2}$:
\begin{align*}
x_1 & = 1.64 &
f(x_0) & = -0.07 \\
& &
f'(x_0) & = -0.997 \\
\end{align*}
}
\only<7->{%
\begin{equation*}
{\color{red} x_2 \leftarrow 1.57}
\end{equation*}
}
\end{overlayarea}
\end{column}
\begin{column}[T]{.5\textwidth}
\hspace{-1cm}
\resizebox{1.1\textwidth}{!}{
\begin{tikzpicture}
\begin{axis}[%
axis x line=center,
axis y line=center,
ymin=-0.2,ymax=0.8,
xmin=0.0,xmax=1.8,
title=$\cos(x)$
]
\addplot[thick,black,domain=0:3.14] {
cos(deg(x))
};
\only<2-3>{
\addplot[thick,blue,only marks,mark=+,mark size=5pt] coordinates {
( 1, 0.54030 )
};
}
\only<3-4>{
\addplot+[thick,blue,domain=0:3.14,mark=none] {
0.54030 - 0.84147 * (x - 1)
};
}
\only<4-5>{
\addplot[thick,red,only marks,mark=o,mark size=5pt] coordinates {
( 1.6421, 0.0 )
};
\draw[thick,red] (axis cs:1.6421, 0) -- (axis cs:1.6421, -0.0712);
}
\only<5-6>{
\addplot[thick,blue,only marks,mark=+,mark size=5pt] coordinates {
( 1.6421, -0.0712 )
};
}
\only<6-7>{
\addplot+[thick,blue,domain=0:3.14,mark=none] {
-0.0712 - 0.99746 * (x - 1.6421)
};
}
\only<7->{
\addplot[thick,red,only marks,mark=o,mark size=5pt] coordinates {
( 1.5707, 7.8e-5 )
};
}
\end{axis}
\end{tikzpicture}
}
\vfill
\begin{overlayarea}{\textwidth}{.3\textheight}
\only<8->{
\begin{tabular}{c|cc}
$k$ & $x_k$ & $f(x_k)$\\
\hline
$0$ & $1.00$ & $+0.54$\\
$1$ & $1.64$ & $-0.07$\\
$2$ & $1.57$ & $+0.00$\\
\hline
$\xi$&$1.5708$ & $0$
\end{tabular}
}
\end{overlayarea}
\end{column}
\end{columns}
\end{frame}
\begin{frame}{Modified Newton-Raphson method \hfill ($N=2$)}
\begin{block}{Expansion for $N=2$}
\vspace{-2\medskipamount}
\begin{equation}
\tag{*}
f(\xi) = 0 = f(x_0) + (\xi - x_0) \cdot f'(x_0) + \frac{1}{2} (\xi - x_0)^2 f''(x_0) + \mathcal{O}((\xi - x_0)^3)
\end{equation}
\end{block}
\vfill
Algorithm:
\begin{enumerate}
\item start with index $k=0$
\item solve equation (*), assuming a vanshing approximation error:
\begin{align*}
a_\pm
& \leftarrow x_k - \frac{f'(x_k) \pm \sqrt{[f'(x_k)]^2 - 2 f(x_k) f''(x_k)}}{f''(x_k)}\\
& \approx \xi + \mathcal{O}((\xi - x_k)^3)
\end{align*}
\item if $|f(a_+)| < |f(a_-)|$, then $x_{k + 1} \leftarrow a_+$; otherwise $x_{k + 1} \leftarrow a_-$
\item if $f(x_{k + 1}) \leq t$, where $t$ is an a-priori threshold, then stop; otherwise
jump back to step \#2.
\end{enumerate}
\end{frame}
\begin{frame}{Rate of convergence}
Does this iterative procedure converge? If yes, how fast? How can we quantify the rate of convergence?
\begin{itemize}
\item we have \alert{local} convergence of order $p \geq 1$, if for all $x \in U(\xi)$
\begin{equation*}
|| \Phi[f](x) - \xi|| \leq C \cdot || x - \xi||^p\,,\quad\text{with}\quad C \geq 0\,.
\end{equation*}
Note: if $p = 1$ then we must have $C < 1$.
\item we have \alert{global} convergence if $U(\xi) = \mathbb{R}$
\item for $D=1$ we can calculate $p$ if $\Phi[f](x)$ is differentiable to sufficient degree:
\begin{align*}
\Phi(x) - \xi = \Phi(x) - \Phi(\xi) = \frac{(x - \xi)^p}{p!} + o(||x - \xi||^p)
\end{align*}
\end{itemize}
\end{frame}
\begin{frame}{Rate of convergence for Newton-Raphson}
We have $\Phi(x) \equiv \Phi[f](x)$:
\begin{itemize}
\item $\Phi(\xi) = \xi$ by construction
\item $\Phi'(\xi) = \frac{f(\xi) \cdot f''(\xi)}{[f'(\xi)]^2} = 0$ \alert{(*)} by construction
\item $\Phi''(\xi) = \frac{f''(\xi)}{f'(\xi)}$
\end{itemize}
We can therefore write
\begin{equation}
\Phi(x) - \Phi(\xi) = \frac{(x - \xi)^2}{2!} \Phi''(\xi) + o(||x - \xi||^2)\,.
\end{equation}
The Newton-Raphson method converges therefore \emph{at least quadratically}, i.e.: it is a
second-order method.
\vfill
\alert{(*): only if $f'(\xi) \neq 0$}, which is equivalent to $\xi$ is a simple root of $f$.
\end{frame}
\begin{frame}{Pathological Example}
\begin{columns}
\begin{column}[T]{.5\textwidth}
\begin{overlayarea}{\textwidth}{\textheight}
$\bm{k=0}$:
\only<2->{%
\begin{equation*}
x_0 := 0.5
\end{equation*}
}
\only<3->{%
$\bm{k=1}$:
\begin{align*}
x_0 & = 0.5 &
f(x_0) & = +0.707 \\
& &
f'(x_0) & = -0.707 \\
\end{align*}
}
\only<4->{%
\begin{equation*}
{\color{red} x_1 \leftarrow -0.5}
\end{equation*}
}
\only<6->{%
$\bm{k=2}$:
\begin{align*}
x_1 & = -0.5 &
f(x_0) & = +0.707 \\
& &
f'(x_0) & = -0.707 \\
\end{align*}
}
\only<7->{%
\begin{equation*}
{\color{red} x_2 \leftarrow +0.5}
\end{equation*}
}
\end{overlayarea}
\end{column}
\begin{column}[T]{.5\textwidth}
\hspace{-1cm}
\resizebox{1.1\textwidth}{!}{
\begin{tikzpicture}
\begin{axis}[%
axis x line=center,
axis y line=center,
ymin=-0.1,ymax=0.9,
xmin=-0.7,xmax=0.7,
title=$\sqrt{|x|}$,
samples=160
]
\addplot[thick,black,domain=-1:+0] {
sqrt(abs(x))
};
\addplot[thick,black,domain=-0:+1] {
sqrt(abs(x))
};
\only<2-3>{
\addplot[thick,blue,only marks,mark=+,mark size=5pt] coordinates {
( +0.5, 0.707 )
};
}
\only<3-4>{
\addplot+[thick,blue,domain=-1:1,mark=none] {
0.707 + 0.707 * (x - 0.5)
};
}
\only<4-5>{
\addplot[thick,red,only marks,mark=o,mark size=5pt] coordinates {
( -0.5, 0.707 )
};
\draw[thick,red] (axis cs:-0.5, 0) -- (axis cs:-0.5, +0.707);
}
\only<5-6>{
\addplot[thick,blue,only marks,mark=+,mark size=5pt] coordinates {
( -0.5, 0.707 )
};
}
\only<6-7>{
\addplot+[thick,blue,domain=-1:1,mark=none] {
0.707 - 0.707 * (x + 0.5)
};
}
\only<7->{
\addplot[thick,red,only marks,mark=o,mark size=5pt] coordinates {
( 0.5, 0.707 )
};
\draw[thick,red] (axis cs:+0.5, 0) -- (axis cs:+0.5, +0.707);
}
\end{axis}
\end{tikzpicture}
}
\vfill
\begin{overlayarea}{\textwidth}{.3\textheight}
~\\\bigskip
endless loop
\end{overlayarea}
\end{column}
\end{columns}
\end{frame}
\begin{frame}{Pitfalls of Newton-Raphson}
We had previously discussed that Newton-Raphson is guaranteed to converge
in an environment $U(\xi)$ of the root $\xi$.\\ \medskip
However, it may fail to converge or converge only very slowly\dots
\begin{itemize}
\item if $U(\xi)$ is a small interval and the initial point lies outside of $U$,
\item if the algorithm encounters a stationary point $\chi$ of $f$ (i.e. $f'(\chi) = 0$),
\item if $\xi$ is a $m$-multiple root of $f$ (i.e. $f(x) \sim (x - \xi)^m \dots$)
\end{itemize}
The last point can be accomodated for: If you know the mulitplicity $m$ of the root,
modify the Newton-Raphson step to read:
\begin{equation*}
x_{k + 1} \leftarrow x_{k} - m \frac{f(x_k)}{f'(x_k)}
\end{equation*}
Basically, this rescales the derivative from $f'(x) \to f'(x) / m$.
\end{frame}
\begin{frame}{Honorable mention: Horner Scheme}
If $f(x) \equiv p_n(x)$ is a polynomial of degree $n$ in $x$, there might be multiple roots of $f$, i.e.
roots $\xi_{n}$ with $f'(\xi_n) = 0$.\\
\vfill
The Horner Scheme allows to efficiently calculate all the (multiple) roots
of the polynomial $p_n(x)$.\\
\vfill
However, this is of limited use: The polynomials dominantly arise in numerical computations
as the characteric polynomial $\chi_M$ of a matrix $M$. In these cases, the roots of $\chi_M$
correspond to eigenvalue of $M$. However, the are better/more stables ways to \alert{numerically} compute
all eigenvalues of $M$, instead of finding all root of $\chi_M$.
\end{frame}
\begin{frame}{Intermission: Derivatives}
Both the original and the modified Newton-Raphson methods involve taking
the derivative $f'$ of our target function $f$ at arbitrary points. Of
course, if we can analytically calculate the derivative, there is no problem.\\
\vfill
What about functions that we can only evaluate numerically? (Think: Experiments,
numerical calculations)\\
\vfill
Enter the field of \emph{numerical differentiation}!
\end{frame}
\begin{frame}{Intermission: Derivatives}
How can we obtain the derivative $f'(x)$ numerically? Let's start from an expansion of
$f(x)$ close to the point of interest $x_0$:
\begin{tikzpicture}
\node at (0,0) {
\begin{minipage}{\textwidth}
\begin{align*}
f(x_0 + 2 h) & = f(x_0) + f'(x_0) \cdot 2 h + \frac{f''(x_0)}{2} \cdot 4 h^2 + \dots\\
f(x_0 + \phantom2 h) & = f(x_0) + f'(x_0) \cdot \phantom2 h + \frac{f''(x_0)}{2} \cdot \phantom4 h^2 + \dots\\
f(x_0 {\color{white}+ 2h}) & = f(x_0)\\
f(x_0 - \phantom2 h) & = f(x_0) - f'(x_0) \cdot \phantom2 h + \frac{f''(x_0)}{2} \cdot \phantom4 h^2 + \dots\\
f(x_0 - 2 h) & = f(x_0) - f'(x_0) \cdot 2 h + \frac{f''(x_0)}{2} \cdot 4 h^2 + \dots
\end{align*}
\end{minipage}
};
%\draw[step=0.50,black,thin] (-4.5,-2.0) grid (4.5,2.0);
\draw[green!75!black,thick,rounded corners=5pt] (-4.5,+0.15) rectangle (+4.5,+1.05);
\only<1>{
\draw[red, thick,rounded corners=5pt] (-4.5,-0.5) rectangle (+4.5,+0.1);
}
\only<2>{
\draw[red, thick,rounded corners=5pt] (-4.5,-1.5) rectangle (+4.5,-0.5);
}
\end{tikzpicture}
\begin{overlayarea}{\textwidth}{.4\textheight}
\only<1>{
\begin{block}{Variant \#1: Forward difference quotient}
\begin{equation*}
\frac{{\color{green!75!black}f(x_0 + h)} - {\color{red!90!black}f(x_0)}}{h} = f'(x_0) + \frac{f''(x_0)}{2} \cdot h = f'(x_0) + o(h)
\end{equation*}
\end{block}
}
\only<2>{
\begin{block}{Variant \#2: Central difference quotient}
\vspace{-1.5\medskipamount}
\begin{equation*}
\frac{{\color{green!75!black}f(x_0 + h)} - {\color{red!90!black}f(x_0 - h)}}{2h} = f'(x_0) + \frac{f'''(x_0)}{3} h^2 = f'(x_0) + o(h^2)
\end{equation*}
\end{block}
}
\end{overlayarea}
\end{frame}
\begin{frame}{Intermission: How to choose $h$?}
\begin{columns}
\begin{column}{.5\textwidth}
\includegraphics[width=\textwidth]{eps.png}\\
{\tiny \hfill source: Wikipedia}
\begin{itemize}
\item two sources of numerical errors compete:
\begin{itemize}
\item formula errors, and
\item round-off (finite precision) errors
\end{itemize}
\end{itemize}
\end{column}
\begin{column}{.5\textwidth}
\begin{itemize}
\item formula error can be reduced through ``intelligent formulae''
\item round-off error can be reduced through (schematically):
\begin{align*}
x_\text{approx}
& \leftarrow x + h\\
h_\text{approx}
& \leftarrow x_\text{approx} - x\\
f' & \leftarrow \frac{f(x_\text{approx}) - f(x)}{h_\text{approx}}
\end{align*}
\item optimal choice of $h$ will be close to
\begin{equation*}
h \sim \sqrt{\varepsilon} \cdot x
\end{equation*}
with $\varepsilon$: working precision
\end{itemize}
\end{column}
\end{columns}
\end{frame}
\begin{frame}{Intermission: Higher derivatives}
What about $f''(x_0)$, $f^{(3)}(x_0)$ or even higher derivatives?\\
\medskip
\begin{block}{Central difference quotient for $f^{(n)}$}
\begin{align*}
\sum_{k=0}^n (-1)^k \binom{n}{k} \frac{f(x_0 + (n - 2 k) h)}{(2 h)^n}
& = f^{(n)}(x_0) + o(h^2)
\end{align*}
\end{block}
\end{frame}
\newcommand{\del}{\partial}
\begin{frame}{Intermission: Multiple partial derivatives}
What about (multiple) partial derivatives $\del_x \del_y f(x, y)$?\\
\medskip
\begin{block}{Central difference ``stencil'' for $\del_x \del_y f(x, y)$}
Combine what we learned previously:
\begin{multline*}
\frac{\only<1>{a}\only<2>{+} f(+ h_x, + h_y) \only<1>{+b}\only<2>{-} f(+ h_x, - h_y) \only<1>{+c}\only<2>{-} f(- h_x, + h_y) \only<1>{+d}\only<2>{+} f(- h_x, - h_y)}{4 h_x h_y}\\
= \only<1>{\dots}\only<2>{\del x\del_y f(x,y)|_{x=0,y=0} + o(h_x^2) + o(h_y^2)}
\end{multline*}
\end{block}
\vfill
\only<2>{End of intermission.}
\end{frame}
\begin{frame}{Newton-Raphson for $D > 1$?}
What happens if we generalize to $D > 1$?\\
\medskip
Let's consider the case $D=2$ first, for a polynomial of degree $1$
\begin{equation*}
f(x_1, x_2) = (1 - x) \cdot (2 - y)
\end{equation*}
\medskip
$f$ has infinitely many roots:
\begin{itemize}
\item $x = 1$ and $y$ arbitrary
\item $y = 2$ and $x$ arbitrary
\end{itemize}
In order to pin-point only one point in $D=2$, we need to solve $2$
equations:
\begin{align*}
f_1(x_1, x_2) & = 0, &
f_2(x_1, x_2) & = 0
\end{align*}
In $D$ dimensions, we will need to solve $D$ simultaneous equations!
\end{frame}
\begin{frame}{Generalization to $D > 1$}
Substitute:
\begin{align*}
f(x) \in \mathbb{R} \quad & \to \quad f(x) \in \mathbb{R}^D\\
x \in \mathbb{R} \quad & \to \quad x \in \mathbb{R}^D\\
f'(x) \in \mathbb{R} \quad & \to \quad [D_f]_{ij} \equiv \frac{\del f_i(x)}{\del x_j} \in \mathbb{R}^{D\times D}
\end{align*}
Then the Newton-Raphson steps reads
\begin{equation*}
x^{(k+1)} \leftarrow x^{(k)} - \left[D_f(x^{(k)})\right]^{-1} f(x^{(k)})\,.
\end{equation*}
Problems:
\begin{itemize}
\item the functional matrix $D_f(x^{(k)})$ might be singular, or numerically close to being
singular
\item inverting $D_f$ is generally computationally expensive!
\end{itemize}
In the next two lectures, Marcin will discuss how to solve systems of linear equations.
These methods will also introduce you to methods to numerically compute the above $x^{(k + 1)}$
from $D_f(x^{(k)})$.
\end{frame}
\backupbegin
\begin{frame}\frametitle{Backup}
\end{frame}
\backupend
\end{document}
Danny van Dyk