Can an approximated Jacobian with finite differences cause instability in the Newton method?

Question

I have implemented a backward-Euler solver in python 3 (using numpy). For my own convenience and as an exercise, I also wrote a small function that computes a finite difference approximation of the gradient so that I don't always have to determine the Jacobian analytically (if it is even possible!).

Using the descriptions provided in Ascher and Petzold 1998, I wrote this function which determines the gradient at a given point x:

def jacobian(f,x,d=4):
    '''computes the gradient (Jacobian) at a point for a multivariate function.

    f: function for which the gradient is to be computed
    x: position vector of the point for which the gradient is to be computed
    d: parameter to determine perturbation value eps, where eps = 10^(-d).
        See Ascher und Petzold 1998 p.54'''

    x = x.astype(np.float64,copy=False)
    n = np.size(x)
    t = 1 # Placeholder for the time step
    jac = np.zeros([n,n])
    eps = 10**(-d)
    for j in np.arange(0,n):
        yhat = x.copy()
        ytilde = x.copy()
        yhat[j] = yhat[j]+eps
        ytilde[j] = ytilde[j]-eps
        jac[:,j] = 1/(2*eps)*(f(t,yhat)-f(t,ytilde))
    return jac

I tested this function by taking a multivariate function for the pendulum and comparing the symbolic Jacobian to the numerically determined gradient for a range of points. I was pleased with the results of the test, the error was around 1e-10. When I solved the ODE for the pendulum using the approximated Jacobian, it also worked very well; I couldn't detect any difference between the two.

Then I tried testing it with the following PDE (Fisher's equation in 1D):

$$ \partial_t \mathbf{u} = \partial_x(k\partial_x\mathbf{u}) + \lambda(\mathbf{u}(C-\mathbf{u}))$$

using a finite difference discretization.

Now Newton's method blows up in the first timestep:

/home/sfbosch/Fisher-Equation.py:40: RuntimeWarning: overflow encountered in multiply
  du = (k/(h**2))*np.dot(K,u) + lmbda*(u*(C-u))
./newton.py:31: RuntimeWarning: invalid value encountered in subtract
  jac[:,j] = 1/(2*eps)*(f(t,yhut)-f(t,yschlange))
Traceback (most recent call last):
  File "/home/sfbosch/Fisher-Equation.py", line 104, in <module>
    fisher1d(ts,dt,h,L,k,C,lmbda)
  File "/home/sfbosch/Fisher-Equation.py", line 64, in fisher1d
    t,xl = euler.implizit(fisherode,ts,u0,dt)
  File "./euler.py", line 47, in implizit
    yi = nt.newton(g,y,maxiter,tol,Jg)
  File "./newton.py", line 54, in newton
    dx = la.solve(A,b)
  File "/usr/lib64/python3.3/site-packages/scipy/linalg/basic.py", line 73, in solve
    a1, b1 = map(np.asarray_chkfinite,(a,b))
  File "/usr/lib64/python3.3/site-packages/numpy/lib/function_base.py", line 613, in asarray_chkfinite
    "array must not contain infs or NaNs")
ValueError: array must not contain infs or NaNs

This happens for a variety of eps values, but strangely, only when the PDE spatial step size and time step size are set so that the Courant–Friedrichs–Lewy condition is not met. Otherwise it works. (This is the behaviour you'd expect if solving with forward Euler!)

For completeness, here is the function for the Newton Method:

def newton(f,x0,maxiter=160,tol=1e-4,jac=jacobian):
    '''Newton's Method.

    f: function to be evaluated
    x0: initial value for the iteration
    maxiter: maximum number of iterations (default 160)
    tol: error tolerance (default 1e-4)
    jac: the gradient function (Jacobian) where jac(fun,x)'''

    x = x0
    err = tol + 1
    k = 0
    t = 1 # Placeholder for the time step
    while err > tol and k < maxiter:
        A = jac(f,x)
        b = -f(t,x)
        dx = la.solve(A,b)
        x = x + dx
        k = k + 1
        err = np.linalg.norm(dx)
    if k >= maxiter:
        print("Maxiter reached. Result may be inaccurate.")
        print("k = %d" % k)
    return x

(The function la.solve is scipy.linalg.solve.)

I am confident that my backward Euler implementation is in order, because I have tested it using a function for the Jacobian and get stable results.

I can see in the debugger that newton() manages 35 iterations before the error occurs. This number remains the same for every eps I have tried.

An additional observation: when I compute the gradient with FDA and a function using the initial condition as an input and compare the two while varying the size of epsilon, the error grows as epsilon shrinks. I would expect it to be large at first, then get smaller, then larger again as epsilon shrinks. So an error in my implementation of the Jacobian is a reasonable assumption, but if so, it's so subtle that I am unable to see it. EDIT: I modified jacobian() to use forward instead of central differences, and now I observe the expected development of the error. However, newton() still fails to converge. Observing dx in the Newton iteration, I see that it only grows, there is not even a fluctuation: it nearly doubles (factor 1.9) with each step, with the factor getting progressively larger.

Ascher and Petzold do mention that difference approximations for the Jacobian do not always work well. Can an approximated Jacobian with finite differences cause instability in Newton's method? Or is the cause somewhere else? How else might I approach this problem?

Answer 1

More of a long comment than anything else:

Using the descriptions provided in Ascher and Petzold 1998, I wrote this function which determines the gradient at a given point x:

Look at the code for SUNDIALS' difference quotient approximation to get a better idea of what you should do in an implementation. Ascher and Petzold is a good book to get started, but SUNDIALS is actually used in production work, and thus has been better tested. (Also, SUNDIALS is related to DASPK, which Petzold worked on.)

Ascher and Petzold do mention that difference approximations for the Jacobian do not always work well. Can an approximated Jacobian with finite differences cause instability in Newton's method?

Empirically, approximate Jacobians can cause convergence failures in Newton's method. I don't know that I'd characterize them as "instability"; in some cases, it's just not possible to achieve the desired error tolerances in the termination criteria. In other cases, it could manifest as an instability. I'm almost certain there's a more quantitative result on this phenomenon in Higham's numerical methods book, or Hairer and Wanner's discussion of W-methods.

Or is the cause somewhere else? How else might I approach this problem?

It depends where you think the error might be. If you're extremely confident in your implementation of backward Euler, I wouldn't start there. Experience has made me paranoid in my implementations of numerical methods, so if it were me, I'd start by coding up a few really basic test problems (a couple nonstiff and stiff linear problems, the heat equation with a centered finite difference approximation, things like that) and I'd use the method of manufactured solutions to assure myself that I know what the solution will be, and what I should be comparing against.

However, you've already done some of that:

I am confident that my backward Euler implementation is in order, because I have tested it using a function for the Jacobian and get stable results.

That would be the next thing I would test: use an analytical Jacobian. After that, you might also look at the extremal eigenvalues of your finite difference Jacobian on the off chance you're in the unstable region of backward Euler. Looking at extremal eigenvalues of your analytical Jacobian as a basis for comparison might give you some insight. Assuming those all check out, the problem is probably in the Newton solve.