scipy.optimize.root as a nonlinear solver

[gdmcbain]

#118 referred to the translation of one of the FEniCS tutorial examples into scikit-fem https://github.com/gdmcbain/fenics-tuto-in-skfem/tree/master/05_poisson_nonlinear. As discussed there, it's a nonlinear Poisson equation like docs/examples/ex10.py.

In ex10, Newton iteration is coded from scratch. I wondered whether it would be possible to use a canned implementation such as scipy.optimize.root. This was done.

Actually in the scikit-fem translation, the optional Jacobian is not passed to root; it solves in a reasonable time (4 seconds) without for this very small problem. It shouldn't be that hard to add the Jacobian (by hand as suggested in #27 or perhaps using SymPy #118).

Would ex10 be better for using a canned nonlinear solver? I tend to prefer to encapsulate algorithms like Newton iteration than recode them each time.

SciPy is already a dependency.

[kinnala]

Assumed that scipy implements it in a reasonable fashion (with regards to the sparsity etc.), I don't mind that the example is modified.

[gdmcbain]

Yes. I'm not immediately sure how to check that, but will look into it. Perhaps repeatedly refine with each method until the machine falls over? Or I could read the source.

This touches on a big area which hasn't been addressed much yet: scaling up. Most or all of the examples to date only demonstrate formal correctness, whereas very often the approaches that work best for small problems have to be swapped to scale up; e.g. switching from direct linear algebraic solvers to iterative and the need for preconditioning. I wonder whether there are good testing frameworks for this kind of thing, that would record not only correctness but speed and the consumption of memory. As these wouldn't have to run on-line, they might not introduce an on-going dependency. In general I do think that it is reasonable to offer multiple ways of solving the same problem, one method being simpler and clearer if only suitable for small problems, another that scales better but is more complicated. Of course iterative methods are usually slower than direct where the latter are applicable. And small problems do arise in practice; e.g. most one- and many two-dimensional problems can be treated as small. I shall give this general matter some thought.

[gdmcbain]

More specifically, another thing is that scipy.optimize.root is more an interface than a definite routine. SciPy ships with ten different implements, selected by the method keyword argument, and each of those take various options. I suspect these behave very differently with regard to speed and usage of memory. In particular, in 05_poisson_nonlinear, I didn't even bother passing the Jacobian—it wasn't necessary, the scikit-fem solution is obtained in seconds anyway. For most problems, I imagine passing the Jacobian would be worth the additional analytical and programming effort; that has already been done for ex10.

[gdmcbain]

On further reflection, I think that the use of a canned numerical continuation algorithm, e.g.

scikit-fem/docs/examples/ex23.py

Lines 85 to 87 in c8ad891

	pacopy.euler_newton(
	problem, u0, lmbda0, callback, max_steps=500, newton_tol=1.0e-10
	)

is better than a simple Newton iteration. The former will typically contain a robust implementation of the latter and in real examples:

• Numerical continuation in a parameter governing the nonlinearity is often required to obtain an initial guess within the radius of convergence of the Newton iteration.
• The solution of the nonlinear problem across a range of values of the parameter is very often of interest.

For example: the Reynolds number for the backward facing step (ex24).

[gdmcbain]

Taking another look at this, there is an unforeseen difficulty with scipy.optimize.rootas described in ‘Which SciPy nonlinear solver when Jacobian is analytically known and sparse?’:

I am surprised that I can't find a method in scipy.optimize that takes a Jacobian in sparse matrix format.

Indeed, passing a function jac that returns a scipy.sparse.csr_matrix raises

*** TypeError: fsolve: there is a mismatch between the input and output shape of the 'fprime' argument 'derivative'.Shape should be (1089, 1089) but it is (1,).

which is because root calls numpy.atleast_1d on the Jacobian csr_matrix returned by jac

https://github.com/scipy/scipy/blob/41f26a297baf7bf2ce7b480f58bd96bebd6b6071/scipy/optimize/minpack.py#L26

which mucks things up by turning that into an array of one object which is the Jacobian.

array([<1089x1089 sparse matrix of type '<class 'numpy.float64'>'
with 7361 stored elements in Compressed Sparse Row format>], dtype=object)

Perhaps this was forewarned in the stipulation

that scipy implements it in a reasonable fashion (with regards to the sparsity etc.)

[gdmcbain]

Following the (offline) discussion of the 'action' of a bilinear form (Kirby & Logg, 'Finite element variational forms') and then the tantalizing last paragraph of Hannukainen & Jantunen about the then planned future work on 'matrix free implementation of the finite element assembly' (did they ever get to that?), I was reminded of ‘Jacobian-free Newton–Krylov methods’ and took another look at this using root with method='krylov'. It works quite nicely on ex10: just omit the BilinearForm and replace the for-loop with

def residual(u: np.ndarray) -> np.ndarray:
    res = asm(rhs, basis, w=basis.interpolate(u))
    res[D] = 0.
    return res

x = root(residual, x, method='krylov', options={'disp': True}).x

Complete file: ex10_jfnk.py.txt