I asked this question a few days ago on stackoverflow, but I figure scicomp.stackexchange is probably a better place. Sorry for the double post.
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 ?
EDIT : A partial answer to this question can be found here.