在 scipy.root 中使用 LinearOperator 作为 Jacobian

Question

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

我在这里错过了什么？

Answer 1

部分答案：

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

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

长话短说，这给出了远非最佳的糟糕结果，但在更少的优化周期内实现了局部收敛。每个优化周期的成本（对于我手头的个人问题）远远大于“标准”krylov lgmres，所以最终即使接近最佳，我也不能说这是值得的。

说实话，scipy.optimize.root的'krylov'方法的雅可比有限差分逼近让我印象深刻。

在 scipy.root 中使用 LinearOperator 作为 Jacobian

Using a LinearOperator as Jacobian in scipy.root

scipy-optimize