# Implementation of the Jacobian-free Newton method

In my calculation (of a simple heat equation, for testing) using the Newton method, I tried to replace the full Jacobian matrix with an approximation vector, i.e. replacing $$J$$ in $$J(u)\delta u=-F(u)$$ with $$J\vec{v}\approx\frac{F(\vec{u}+\varepsilon \vec{v}) - F(\vec{u})}{\varepsilon}$$ which should be possible, after I solve the system using a GMRES-solver. Here $$\varepsilon$$ is calculated using $$\varepsilon=\frac{1}{n\vert\vert\vec{v}\vert\vert_2}\sum_{i=1}^n\left(b\vert u_i\vert+b\right)\text{, }b=10^{-6}$$ $$F(u)$$ is the function itself, i.e. $$F(u)=\nabla^2u+f$$ Now I encountered two problems:

• Usually in the first iteration $$\vert\vert\vec{v}\vert\vert_2=0$$, resulting in infinity for $$\varepsilon$$. My current solution is to set $$\varepsilon=10^{-6}$$ in that case, but is that correct? I could not find any solution for that in the literature
• When considering the value for $$F(u)$$ at each step, it decreases (as it should) when using the full version of $$J$$, but increases when using the approximation. Is there a way how can I narrow down the possible problems here?

How can I solve those problems (or is there literature about them which I could not find yet)?

• GMRES does not need the matrix itself, only a response in the form of a matrix-vector product. This matrix-vector product can be approximated using a finite difference resulting in a vector (also called the jacobian-free newton-krylov method). Dec 12, 2018 at 9:53
• You have a linear PDE. Why would you compute the application of the matrix via finite differences when you know what the exact value is? Dec 12, 2018 at 21:57
• @WolfgangBangerth: Because I use that equation as test case for the solver. My final application contains a non-linear system of equations. Before applying the method to the more complex stuff, I would like to calculate something with a known solution. Dec 13, 2018 at 7:08

Not sure where you get your equation for $$\epsilon$$, but ultimately your approximation for the Jacobian matvec operation is a finite difference approximation to the directed derivative of $$F(\cdot)$$. This means that you will want as accurate of an approximation to this directed derivative as possible while being sure to avoid numerical round-off issues.

For this situation, I have personally used the following method referred to as the Complex-Step Derivative Approximation:

\begin{align} \nabla F(\boldsymbol{x}) \cdot \boldsymbol{v} = \frac{\text{Im}(F(\boldsymbol{x} + i \epsilon \boldsymbol{v}))}{\epsilon} + O(\epsilon^2) \end{align}

where $$\text{Im}(\cdot)$$ maps a complex vector in $$\mathbb{C}^n$$ to $$\mathbb{R}^n$$ by just taking the imaginary part of each component. Note that if $$\boldsymbol{v}$$ is not a unit vector, you should compute the above quantity with the direction of $$\boldsymbol{v}$$ and then multiply the result by $$\left\lVert \boldsymbol{v}\right\rVert$$. The obvious thing to note here is that the above assumes you can modify your code to allow for feeding in complex numbers into $$F(\cdot)$$. If you can do this, this approximation has great numeric precision since it does not have round-off issues related to subtractions, like typical Finite Difference approximations. You can readily choose $$\epsilon = 10^{-8}$$ and get close to machine precision.

• Unfortunately my code does not allow the usage of imaginary values directly. Of course I can implement that, but I would like to avoid doing that... Dec 12, 2018 at 16:02
• @spektr It's called "Complex-Step Derivative Approximation". You might want to name the method if you are introducing advanced techniques. Dec 12, 2018 at 16:31
• @P.G. I didn’t know the name, I read the paper about it years ago. I’ll update my response. Dec 12, 2018 at 20:52
• +1$\,$ I'm a huge advocate of the Complex Step Approximation, but there are some subtle caveats. Since you will most likely be using a language with built-in support for complex arithmetic, you must be careful about using its built-in functions. For example, in Julia X' is simply the transpose when working with real variables, but automagically becomes the hermitian conjugate when X is promoted to a complex variable. This will ruin the CSA approximation. Analytic functions such as exp(X), sin(X), X^2, etc are fine, but another common built-in to avoid is norm(X)
– greg
Feb 15 at 15:28
• However, norm(X) can be replaced by something like sqrt(sum(X.*X)) and X' by permutedims(X)
– greg
Feb 15 at 15:28

I finally found the solution to that problem: In order to be able to reuse as much code as possible, I had one single function for calculating the right hand side of my problem, i.e. $$-F(u)$$, and for calculating the residual value. This function returned $$-F(u)$$ when calling, but I forgot to remove the minus-sign when using the same function for the approximation, thus resulting in $$J\vec{v}\approx\frac{-F(\vec{u}+\varepsilon \vec{v}) + F(\vec{u})}{\varepsilon}$$ After multiplying everything with $$-1$$, the approach works as expected.