# 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). – arc_lupus Dec 12 '18 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? – Wolfgang Bangerth Dec 12 '18 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. – arc_lupus Dec 13 '18 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... – arc_lupus Dec 12 '18 at 16:02
• @arc_lupus well at least if nothing else works, the above probably will if you are willing to implement it, so it would be a back up approach. Sticking with the normal finite difference you're doing, where do you get the definition for $\epsilon$ that you were using? Have you tried just using $\epsilon = 10^{-6}$ for every iteration? – spektr Dec 12 '18 at 16:22
• @spektr It's called "Complex-Step Derivative Approximation". You might want to name the method if you are introducing advanced techniques. – P. G. Dec 12 '18 at 16:31
• @spektr: I got it from "Review: Jacobian-free Newton–Krylov methods: a survey of approaches and applications", in "Journal of Computational Physics 193 (2004)". My problem with the implementation is that my whole program does not support complex numbers. If it should support complex numbers, I have to rewrite it from base (which is the reason I would like to avoid that method)... – arc_lupus Dec 12 '18 at 20:48
• @P.G. I didn’t know the name, I read the paper about it years ago. I’ll update my response. – spektr Dec 12 '18 at 20:52

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.