2
$\begingroup$

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.

Improve this question
$\endgroup$
7
  • $\begingroup$ If it is a linear operator it should return a vector of dimensions (ndim, 1). $\endgroup$
    – nicoguaro ♦
    Mar 1 at 19:26
  • $\begingroup$ Hi, thanks for your answer. I however can't really make sense of it. Would you mind further explaining what you mean ? $\endgroup$
    – G. Fougeron
    Mar 2 at 13:22
  • $\begingroup$ 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. $\endgroup$
    – nicoguaro ♦
    Mar 2 at 16:08
  • $\begingroup$ Yes, I agree. But isn't it what my example does ? $\endgroup$
    – G. Fougeron
    Mar 2 at 17:42
  • $\begingroup$ I am not sure. I could not grasp the whole thing in your code. $\endgroup$
    – nicoguaro ♦
    Mar 2 at 17:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.