在 scipy.root 中使用 LinearOperator 作为雅可比行列式

如何解决在 scipy.root 中使用 LinearOperator 作为雅可比行列式

我想使用 scipy.root 求解非线性方程组。出于性能原因，我想使用 LinearOperator 提供系统的 jacobian。但是，我无法让它工作。这是一个使用 Rosenbrock 函数梯度的最小示例，我首先将雅可比矩阵（即 Rosenbrock 函数的 Hessian）定义为 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)

执行后，我收到以下错误：

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,).

我在这里遗漏了什么？

解决方法

部分答案：

我能够将我的“精确”雅可比输入到 scipy.optimize.nonlin.nonlin_solve 中。这真的感觉很hacky。长话短说，我定义了一个继承自 scipy.optimize.nonlin.Jacobian 的类，我在其中定义了“update”和“solve”方法，以便求解器使用我的确切 jacobian。

我希望性能结果因问题而异。让我详细说明我对“几乎”强制函数（即，如果我花时间移除 4 维对称生成器，问题将是强制的）的 ~10k 维临界点求解的经验，其中有许多局部最小值（和因此大概有很多很多关键点）。

长话短说，这给出了远非最佳结果的可怕结果，但在更少的优化周期内实现了局部收敛。每个优化周期的成本（就我手头的个人问题而言）远远高于“标准”krylov lgmres，因此最终即使接近最佳状态，我也不能说这值得麻烦。

老实说，我对 scipy.optimize.root 的 'krylov' 方法的雅可比有限差分近似印象非常深刻。