# Providing Jacobian as LinearOperator in scipy.optimize.root

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.

• If it is a linear operator it should return a vector of dimensions (ndim, 1). Mar 1 at 19:26
• Hi, thanks for your answer. I however can't really make sense of it. Would you mind further explaining what you mean ? Mar 2 at 13:22
• A linear operator does not return a matrix but a function that when you evaluate it gives you the multiplication of the matrix times a vector. Mar 2 at 16:08
• Yes, I agree. But isn't it what my example does ? Mar 2 at 17:42
• I am not sure. I could not grasp the whole thing in your code. Mar 2 at 17:49