如何让scipy newton_krylov使用其他导数近似方法？

查阅scipy官方文档后得知，newton_krylov默认采用前向差分近似导数，但未找到直接切换或自定义导数近似方法的途径。尝试自定义前向差分实现雅可比向量积，封装成LinearOperator传入inner_M参数测试，结果自定义实现的误差比默认大一个数量级，对工具用法存在疑惑。

测试代码 #

初始化代码 #

import numpy as np
from scipy.optimize import newton_krylov
from scipy.sparse.linalg import LinearOperator

# 定义目标函数
def uniform_problem(x, A, b):
    return b - A@x

size = 12

A = np.random.uniform(-1, 1, size=(size, size))
b = np.random.uniform(-1, 1, size=(size, ))

xr = np.random.uniform(-1, 1, size=(size, ))  # 真实根
x0 = np.random.uniform(-1, 1, size=(size, ))  # 初始猜测值

# 构造F(x)，使得F(xr)=0
F = lambda x: uniform_problem(x, A, b) - uniform_problem(xr, A, b)

# 参数设置
max_iter = 10
tol = 1e-3
h = 1e-4

repeats = 5000

自定义前向差分实现 #

# 使用自定义前向差分实现雅可比向量积
def get_jacobian_vector_product_fdf(F, x, v, h=1e-5):
  step = h * v
  return (F(x + step) - F(x)) / h

error1 = 0
for i in range(repeats):
    x = x0.copy()

    lambdaJv = lambda v: get_jacobian_vector_product_fdf(F, x, v, h)
    linear_operator = LinearOperator((size, size), matvec=lambdaJv)

    solution1 = newton_krylov(F, x, method="gmres", inner_maxiter=max_iter, iter=max_iter, callback=None, f_tol=tol, rdiff=h, inner_M=linear_operator)
    error1 += np.linalg.norm(F(solution1))
error1 /= repeats
print(error1)  # 约1.659173186802721

默认方法实现 #

# 使用默认方法
error2 = 0
for i in range(repeats):
    x = x0.copy()
    
    solution2 = newton_krylov(F, x, method="gmres", inner_maxiter=max_iter, iter=max_iter, callback=None, f_tol=tol, rdiff=h)
    error2 += np.linalg.norm(F(solution2))
error2 /= repeats
print(error2)  # 约0.024629534404425796

print(error1/error2)  # 误差相差一个数量级

问题分析与解决 #

核心问题1：对inner_M参数的误解 #

inner_M不是用来替换雅可比向量积的计算逻辑，而是作为Krylov子空间方法（此处为GMRES）的预条件子，作用是预处理线性系统以加速迭代收敛。当你传入自定义的雅可比向量积算子到inner_M时，newton_krylov依然会用默认前向差分计算雅可比向量积，同时叠加你的算子作为预条件子，这相当于给GMRES加了一个完全不合适的预处理，反而干扰了收敛过程。

核心问题2：自定义算子未随迭代更新x #

你的lambdaJv捕获的是循环中初始的x0.copy()值，但newton_krylov在迭代过程中会不断更新当前的x值，而你的自定义算子始终基于初始点计算雅可比向量积，完全不符合牛顿法“用当前点的导数近似进行迭代”的要求，这是误差巨大的直接原因。

正确的处理方式 #

1. 放弃用inner_M替换雅可比向量积：该参数的设计目的是预条件，而非替换导数计算逻辑。
2. 若需自定义导数近似，手动实现牛顿-Krylov框架：newton_krylov并未暴露替换内部导数计算的接口，若要自定义，建议基于scipy.sparse.linalg.gmres手动实现牛顿迭代逻辑，这样可以完全控制雅可比向量积的计算方式。
3. 修正自定义算子的x更新问题（仅作演示，无法直接适配newton_krylov）：

   class DynamicJVP:
       def __init__(self, F, h):
           self.F = F
           self.h = h
           self.current_x = None
       
       def set_x(self, x):
           self.current_x = x.copy()
       
       def matvec(self, v):
           if self.current_x is None:
               raise ValueError("Current x not set")
           step = self.h * v
           return (self.F(self.current_x + step) - self.F(self.current_x)) / self.h
   

内容的提问来源于stack exchange，提问作者Ultrinik