I want to solve a system of nonlinear equations using scipy.root. For performance reason, I want to provide the jacobian of the system using a LinearOperator. However, I cannot get it to work. Here is a minimal example using the gradient of the Rosenbrock function, where I first define the Jacobian (i.e. the Hessian of the Rosenbrock function) as a LinearOperator.
import numpy as np
import scipy.optimize as opt
import scipy.sparse as sp
ndim = 10
def rosen_hess_LO(x):
return sp.linalg.LinearOperator((ndim,ndim) ,matvec = (lambda dx,xl=x : opt.rosen_hess_prod(xl,dx)))
opt_result = opt.root(fun=opt.rosen_der,x0=np.zeros((ndim),float),jac=rosen_hess_LO)
Upon execution, I get the following error :
TypeError: fsolve: there is a mismatch between the input and output shape of the 'fprime' argument 'rosen_hess_LO'.Shape should be (10, 10) but it is (1,).
What am I missing here ?
Partial answer :
I was able to input my "exact" jacobian into scipy.optimize.nonlin.nonlin_solve . This really felt hacky. Long story short, I defined a class inheriting from scipy.optimize.nonlin.Jacobian, where I defined "update" and "solve" method so that my exact jacobian would be used by the solver.
I expect performance results to greatly vary from problem to problem. Let me detail my experience for a ~10k dimensional critial point solve of an "almost" coercive function (i.e. the problem would be coercive if I had taken the time to remove a 4-dimensional symmetry generator), with many many local minima (and thus presumably many many critical points).
Long story short, this gave terrible results far from the optimum, but local convergence was achieved in fewer optimization cycles. The cost of each of those optimization cycle was (for my personal problem at hand) far greater than the "standard" krylov lgmres, so in the end even close to the optimum, I cannot really say it was worth the trouble.
To be honest, I am very impressed with the Jacobian finite difference approximation of the 'krylov' method of scipy.optimize.root.