% Except otherwise noted,
% this work is © 2009 Dmitri Fedorov, under CC-BY-SA license, version 3+.

\chapter{Nonlinear equations}

\section{Introduction}

{\em Non-linear equations} (or {\em
root-finding}\index{root-finding}) is a problem of finding a set of
$n$ variables $\ve{x}=\{x_1,\dots,x_n\}$ which satisfy a system of $n$
(generally) non-linear equations
	\begin{equation}\label{root:eq:f}
f_{i}(x_{1},...,x_{n})=0 \;\Big|_{i=1,\dots,n} \:.
	\end{equation}
In matrix notation the system is written as
	\begin{equation}
\ve{f}(\ve{x})=0 \;,
	\end{equation}
where $\ve{f}\doteq\{f_1,\dots,f_n\}$.

In one-dimensional cases it is generally possible to plot the function
in the region of interest and see whether the graph crosses the $x$-axis.
One can then be sure the root exists and even figure out its approximate
position to start one's root-finding algorithm with.  In multi-dimensions
one generally does not know if the root exists at all, until it is found.

The root-finding algorithms generally proceed by iteration, starting
from some approximate solution and making consecutive steps (hopefully)
in the direction of the suspected root until some convergence criterion
is satisfied.  The procedure is generally not even guaranteed to converge
unless starting from a point close enough to the sought root.

We shall only consider the multi-dimensional case here since i)
the multi-dimensional root-finding is more difficult, and ii) the
multi-dimensional routines can also be used in a one-dimensional case.

\section{Newton's method}
Newton's method (also reffered to as Newton-Raphson method, after Isaac
Newton and Joseph Raphson) is a root-finding algorithm that uses the
first term of the Taylor series of the functions $f_i$ to linearise the
system (\ref{root:eq:f}) in the vicinity of a suspected root. It is one
of the oldest and best known methods and is a basis of a number of more
refined methods.

Suppose that the point $\ve{x}=\{x_1,\dots,x_n\}$ is close to the root.
The Newton's algorithm tries to find the step $\Delta\ve{x}$ which would
move the point towards the root, such that
	\begin{equation}\label{root:eq:fxdx}
f_i(\ve{x}+\Delta\ve{x})=0 \;\Big|_{i=1,\dots,n}\;.
	\end{equation}
The first order Taylor expansion of (\ref{root:eq:fxdx}) gives a system
of linear equations,
	\begin{equation}
f_i(\ve{x})+\sum_{k=1}^n \frac{\partial f_i}{\partial x_k}\Delta x_k =
0 \;\Big|_{i=1,\dots,n}\;,
	\end{equation}
or, in the matrix form,
	\begin{equation}\label{eq:axf}
\ma{J}\Delta\ve{x} = -\ve{f}(\ve{x}),
	\end{equation}
where $\ve{f}(\ve{x})\doteq\{f_0(\ve{x}),\dots,f_n(\ve{x})\}$ and $\ma{J}$
is the matrix of partial derivatives\footnote{in practice if derivatives
are not available analytically one uses finite differences
	$$
\frac{\partial f_i}{\partial x_k}\approx
\frac{f_i(x_1,\dots,x_{k-1},x_k+\delta x,x_{k+1},\dots,x_n)
-f_i(x_1,\dots,x_k,\dots,x_n)} {\delta x}
	$$
with $\delta x\ll s$ where $s$ is the typical scale of the problem
at hand.},
	\begin{equation}\label{root:eq:ns}
J_{ik}\doteq\frac{\partial f_i}{\partial x_k}\:,
	\end{equation}
called the {\em Jacobian matrix}\index{Jacobian matrix}.

The solution $\Delta\ve{x}$ to the linear system
(\ref{eq:axf})---called the Newton's step---gives the approximate
direction and the step-size towards the solution.

The Newton's method converges quadratically if sufficiently close
to the solution. Otherwise the full Newton's step $\Delta\ve{x}$
might actually diverge from the solution. Therefore in practice a more
conservative step $\lambda\Delta\ve{x}$ with $\lambda<1$ is usually
taken. The strategy of finding the optimal $\lambda$ is referred to as
{\em line search}\index{line search}.

It is typically not worth the effort to find $\lambda$ which minimizes
$\|\ve{f}(\ve{x}+\lambda\Delta\ve{x})\|$ exactly, since
$\Delta\ve{x}$ is only an approximate direction towards
the root.  Instead an inexact but quick minimization strategy is usually
used, like the {\em backtracking line search}\index{backtracking} where
one first attempts the full step, $\lambda=1$, and then backtracks,
$\lambda\leftarrow\lambda/2$, until the condition
	\begin{equation}\label{root:eq:c}
\|\ve{f}(\ve{x}+\lambda\Delta\ve{x})\|
	<\left(1-\frac{\lambda}{2}\right)
\|\ve{f}(\ve{x})\|
	\end{equation}
is satisfied. If the condition is not satisfied for sufficiently small
$\lambda_\mathrm{min}$ the step is taken with $\lambda_\mathrm{min}$
simply to step away from this diffictul place and try again.  A typical
algrorithm of the Newton's method with backtracking line search and
condition (\ref{root:eq:c}) is shown in Table~\ref{root:t:n}.

	\begin{table}[htb]
\caption{Newton's root-finding algorithm with backtracking.}
	\label{root:t:n}
	\hrule
	\begin{code}
repeat \\
\hspace*{1em} solve $\ma{J}\Delta\ve{x}=-\ve{f}(\ve{x})$ for $\Delta\ve{x}$ \\
\hspace*{1em} $\lambda\leftarrow 1$ \\
\hspace*{1em} while
	$\big($ $\|\ve{f}(\ve{x}+\lambda\Delta\ve{x})\|
	>\left(1-\frac{\lambda}{2}\right)\|\ve{f}(\ve{x})\|$
	and $\lambda>{1\over 64}$ $\big)$ do $\lambda \leftarrow \lambda/2$ \\
\hspace*{1em} $\ve{x}\leftarrow\ve{x}+\lambda\Delta\ve{x}$ \\
until converged (e.g. $\|\ve{f}(\ve{x})\|<\mathrm{tolerance}$)
	\end{code}
\hrule
	\end{table}

	\begin{table}
\caption{C implementation of Newton's algorithm with backtracking.}
\lstinputlisting{c/roots/newton.c} 
	\end{table}

A somewhat more refined backtracking linesearch is based on
an approximate minimization of the function
	\begin{equation}\label{}
g(\lambda) \doteq \frac12
\| \ve{f}(\ve{x}+\lambda\Delta\ve{x})\|^2 \:
	\end{equation}
using interpolation.  The values $g(0)=\frac12
\|\ve{f}(\ve{x})\|^2$ and $g'(0)=-\|\ve{f}(\ve{x})\|^2$
are already known (check this). If the previous step with
certain $\lambda_\mathrm{trial}$ was rejected, we also have
$g(\lambda_\mathrm{trial})$.  These three quantities allow to build a
quadratic approximation,
	\begin{equation}\label{}
g(\lambda) \approx g(0) + g'(0)\lambda + c\lambda^2 \:,
	\end{equation}
where
	\begin{equation}\label{}
c = \frac{g(\lambda_\mathrm{trial})-g(0)-g'(0)\lambda_\mathrm{trial}}
{\lambda_\mathrm{trial}^2} \:.
	\end{equation}
The minimum of this approximation (determined by the condition
$g'(\lambda)=0$),
	\begin{equation}\label{}
\lambda_\mathrm{next} = - \frac{g'(0)}{2c} \;,
	\end{equation}
becomes the next trial step-size.

The procedure is repeated recursively until either
condition~(\ref{root:eq:c}) is satisfied or the step becomes too small
(in which case it is taken unconditionally in order to simply get away
from the difficult place).

\section{Broyden's quasi-Newton method}
The Newton's method requires calculation of the Jacobian matrix at every
iteration. This is generally an expensive operation.  Quasi-Newton
methods avoid calculation of the Jacobian matrix at the new point
$\ve{x}+\Delta\ve{x}$, instead trying to use certain approximations,
typically rank-1 updates.

Broyden's algorithm estimates the Jacobian $\ma{J}+\Delta\ma{J}$ at the
point $\ve{x}+\Delta\ve{x}$ using the finite-difference approximation,
	\begin{equation}\label{root:eq:b1}
(\ma{J}+\Delta\ma{J})\Delta \ve{x}=\Delta \ve{f} \:,
	\end{equation}
where $\Delta\ve{f}\doteq\ve{f}(\ve{x}+\Delta\ve{x})- \ve{f}(\ve{x})$
and $\ma{J}$ is the Jacobian at the point $\ve{x}$.

The matrix equation (\ref{root:eq:b1}) is under-determined in more than
one dimension as it contains only $n$ equations to determine $n^2$ matrix
elements of $\Delta\ma{J}$. Broyden suggested to choose $\Delta\ma{J}$
as a rank-1 update, linear in $\Delta\ve{x}$,
	\begin{equation}\label{root:eq:b2}
\Delta\ma{J}=\ve{c}\,\Delta \ve{x}^T\:,
	\end{equation}
where the unknown vector $\ve{c}$ can be found by substituting
(\ref{root:eq:b2}) into (\ref{root:eq:b1}), which gives
	\begin{equation}
\Delta\ma{J}=\frac{\Delta\ve{f}-\ma{J}\Delta\ve{x}}{\|\Delta \ve{x}\|^2}
\Delta\ve{x}^T\:.
	\end{equation}

In practice if one wanders too far from the point where $\ma{J}$ was first
calculated the accuracy of the updates may decrease significantly. In
such case one might need to recalculate $\ma{J}$ anew.  For example,
two successive steps with $\lambda_\mathrm{min}$ might be interpreted as
a sign of accuracy loss in $\ma{J}$ and subsequently trigger its recalculation.
