Nonlinear systems of equations

Automatic differentiation

Dual numbers

We all know the field of complex numbers $\mathbb{C}$: they extend the real numbers $\mathbb{R}$ based on the introduction ot $i$ with $i^2=-1$.

Dual numbers are defined by extending the real numbers by formally adding an number $\epsilon$ with ${\epsilon }^2=0$:

$$D=\{a+b\epsilon |a,b\in \mathbb{R}\}=\{\left(\begin{array}{cc}a& b\\ 0& a\end{array}\right)|a,b\in \mathbb{R}\}\subset {\mathbb{R}}^{2\times 2}$$

They form a ring, not a field.

• Evaluating polynomials on dual numbers: Let $p(x)=\sum _{i=0}^n{p}_i{x}^i$. Then

$$\begin{array}{rl}p(a+b\epsilon )& =\sum _{i=0}^n{p}_i{a}^i+\sum _{i=1}^{n}i{p}_i{a}^{i-1}b\epsilon \\ & =p(a)+b{p}^{\prime }(a)\epsilon \end{array}$$

• This can be generalized to any analytical function. $\Rightarrow$ automatic evaluation of function and derivative at once

• $\Rightarrow$ forward mode automatic differentiation

• Multivariate dual numbers: generalization for partial derivatives

Dual numbers in Julia

Constructing a dual number:

Accessing its components:

Comparison with known derivative:

Polynomial expressions:

Standard functions:

Function composition:

Conclusion: if we apply dual numbers in the right way, we can do calculations with derivatives of complicated nonlinear expressions without the need to write code to calculate derivatives.

The forwardiff package provides these facilities.

Let us plot some complicated function:

Solving nonlinear systems of equations

Let $A_1\dots A_n$ be functions depending on $n$ unknowns $u_1\dots u_n$. Solve the system of nonlinear equations:

$$A(u)=\left(\begin{array}{c}{A}_1(u_1\dots u_n)\\ A_2(u_1\dots u_n)\\ ⋮\\ A_n(u_1\dots u_n)\end{array}\right)=\left(\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ f_n\end{array}\right)=f$$

$A(u)$ can be seen as a nonlinar operator $A:D\to {\mathbb{R}}^n$ where $D\subset {\mathbb{R}}^n$ is its domain of definition.

There is no analogon to Gaussian elimination, so we need to solve iteratively.

Fixpoint iteration scheme:

Assume $A(u)=M(u)u$ where for each $u$, $M(u):{\mathbb{R}}^n\to {\mathbb{R}}^n$ is a linear operator.

Then we can define the iteration scheme: choose an initial value $u_0$ and at each iteration step, solve

$$M(u^i)u^{i+1}=f$$

Terminate if

$$||A(u^i)-f||<\epsilon (\text{residual based})$$

or

$$||u_{i+1}-u_i||<\epsilon (\text{update based}).$$

• Large domain of convergence

• Convergence may be slow

• Smooth coefficients not necessary

Example problem

Newton iteration scheme

The fixed point iteration scheme assumes a particular structure of the nonlinear system. Can we do better ?

Let $A^{\prime }(u)$ be the Jacobi matrix of first partial derivatives of $A$ at point $u$:

$$A^{\prime }(u)=(a_{kl})$$

'with

$$a_{kl}=\frac{\partial }{\partial u_l}{A}_k(u_1\dots u_n)$$

The one calculates in the $i$-th iteration step:

$$u_{i+1}=u_i-(A^{\prime }(u_i){)}^{-1}(A(u_i)-f)$$

One can split this a folows:

• Calculate residual: $r_i=A(u_i)-f$

• Solve linear system for update: $A^{\prime }(u_i)h_i=r_i$

• Update solution: $u_{i+1}=u_i-h_i$

General properties are:

• Potenially small domain of convergence - one needs a good initial value

• Possibly slow initial convergence

• Quadratic convergence close to the solution

Linear and quadratic convergence

Let $e_i=u_i-\hat{u}$.

• Linear convergence: observed for e.g. linear systems: Asymptically constant error contraction rate

$$\frac{||e_{i+1}||}{||e_i||}\sim \rho <1$$

• Quadratic convergence: $\mathrm{\exists }{i}_0>0$ such that $\mathrm{\forall }i>i_0$, $\frac{||e_{i+1}||}{||e_i|{|}^2}\le M<1.$

  • As $||e_i||$ decreases, the contraction rate decreases:

$$\begin{array}{rl}\frac{\frac{||e_{i+1}||}{||e_i||}}{\frac{||e_i||}{||e_{i-1}||}}& =\frac{||e_{i+1}||}{\frac{||e_i|{|}^2}{||e_{i-1}||}}\le ||e_{i-1}||M\end{array}$$

• In practice, we can watch $||r_i||$ or $||h_i||$

Automatic differentiation for Newton's method

This is the situation where we could apply automatic differentiation for vector functions of vectors.

Create a result buffer for $n=2$

Calculate function and derivative at once:

A Newton solver with automatic differentiation

Let us take a more complicated example:

Here, we observe that we have to use lots of iteration steps and see a rather erratic behaviour of the residual. After $\approx$ 55 steps we arrive in the quadratic convergence region where convergence is fast.

Damped Newton iteration

There are may ways to improve the convergence behaviour and/or to increase the convergence radius in such a case. The simplest ones are:

• find a good estimate of the initial value

• damping: do not use the full update, but damp it by some factor which we increase during the iteration process

• linesearch: automatic detection of a damping factor

The example shows: damping indeed helps to improve the convergece behaviour. However, if we keep the damping parameter less than 1, we loose the quadratic convergence behavior.

Parameter embedding

Another option is the use of parameter embedding for parameter dependent problems.

• Problem: solve $A(u_{\lambda },\lambda )=f$ for $\lambda =1$.

• Assume $A(u_0,0)$ can be easily solved.

• Choose step size $\delta$

1. Solve $A(u_0,0)=f$

2. Set $\lambda =0$

3. Solve $A(u_{\lambda +\delta },\lambda +\delta )=f$ with initial value $u_{\lambda }$

4. Set $\lambda =\lambda +\delta$

5. If $\lambda <1$ repeat with 3.

• If $\delta$ is small enough, we can ensure that $u_{\lambda }$ is a good initial value for $u_{\lambda +\delta }$.

• Possibility to adapt $\delta$ depending on Newton convergence

• Parameter embedding + damping + update based convergence control go a long way to solve even strongly nonlinear problems!

• As we will see later, a similar apporach can be used for time dependent problems.