在 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' 方法的雅可比有限差分近似印象非常深刻。