When is Newton-Krylov not an appropriate solver?

Question

Recently I have been comparing different non-linear solvers from scipy and was particularly impressed with the Newton-Krylov example in the Scipy Cookbook in which they solve a second order differential equation equation with non-linear reaction term in about 20 lines of code.

I modified the example code to solve the non-linear Poisson equation (also called the Poisson-Boltzmann equation, see page 17 in these notes) for semiconductor heterostructures, which has the form,

$$ \frac{d^2\phi}{dx^2} - k(x) \left(p(x,\phi) - n(x,\phi) + N^{+}(x)\right) = 0 $$

(This is the residual function that is passed to the solver.)

This is an electrostatics problem where $n(x,\phi)$ and $p(x,\phi)$ are nonlinear functions for the form $n_i(x) e^{-(E_i(x,\phi) - E_f)}$. The details here are not important, but the point is that the non-linear function vary exponentially with $\phi$ so the residual function can vary over a huge range ($10^{-6} - 10^{16})$ with a slight change in $\phi$.

I numerically solve this equation with scipy's Newton-Krylov, but it would never converge (in fact it would always report an error with the calculating the Jacobian). I switched from a Newton-Krylov solver to fsolve (which is based on MINPACK hybrd) and it worked first time!

Are there general reasons why Newton-Krylov is not a good fit to certain problems? Do the input equations need to be conditioned somehow?

Maybe more information is needed to comment, but why do you think fsolve worked in this case?

Answer 1

There are two issues that you are likely to be encountering.

Ill-conditioning

First, the problem is ill-conditioned, but if you only provide a residual, Newton-Krylov is throwing away half your significant digits by finite differencing the residual to get the action of the Jacobian:

$$ J[x] y \approx \frac{F(x+\epsilon y) - F(x)}{\epsilon}$$

If you provide an analytic Jacobian, you get to keep all the digits (e.g., 16 in double precision). Krylov methods also rely on inner products, so if your Jacobian is ill-conditioned to the tune of $10^{16}$, it is effectively singular and Krylov can stagnate or return erroneous solutions. This can also prevent direct solvers from converging. Sometimes you can use multigrid methods to homogenize into a coarse grid with tractable conditioning. When a problem cannot be formulated with better conditioning, it can be worth working in quad precision. (This is supported by PETSc, for example.)

Note that the same issues apply to quasi-Newton methods, though without finite differencing. All scalable methods for problems with non-compact operators (e.g., differential equations) must use Jacobian information for preconditioning.

It is likely that fsolve either did not use Jacobian information or that it used a dogleg method or a shift to make progress with a "gradient descent" method, despite an essentially singular Jacobian (i.e., finite differencing would have a lot of "noise" from finite precision arithmetic). This is not scalable and fsolve likely gets slower as you increase the problem size.

Globalization

If the linear problems are solved accurately, we can rule out problems relating to the linear problem (Krylov) and focus on those due to nonlinearity. Local minima and nonsmooth features slow convergence or cause stagnation. Poisson-Boltzmann is a smooth model so if you start with a good enough initial guess, Newton will converge quadratically. Most globalization strategies involve some sort of continuation to produce a high-quality initial guess for the final iterations. Examples include grid continuation (e.g., Full Multigrid), parameter continuation, and pseudotransient continuation. The latter is generally-applicable to steady-state problems, and offers some global convergence theory, see Coffey, Kelley, and Keyes (2003). A search turns up this paper, which may be useful to you: Shestakov, Milovich, and Noy (2002) Solution of the nonlinear Poisson-Boltzmann equation using pseudo-transient continuation and the finite element method. Pseudotransient continuation is closely related to the Levenberg-Marquardt algorithm.

Further reading

• Why is Newton's method not converging?
• Why is my iterative linear solver not converging?