SciPy 的有限差分自适应步长方法起源于哪里？

Question

在 SciPy 的

KrylovJacobian

类中，有这个方法：

def _update_diff_step(self):
    mx = abs(self.x0).max()
    mf = abs(self.f0).max()
    self.omega = self.rdiff * max(1, mx) / max(1, mf)

这与：

相同

这修改了有限差分法使用的步长，但是，我找不到这个表达式的起源，或者它为什么起作用。

有人知道这个方法的来源或背后的原因吗？

Answer 1

阅读 SciPy 引用的期刊文章*，我找不到任何与 SciPy 正在做的事情完全相同的 omega 选择。然而，有几个案例是相似的。

有人知道这个方法的来源或背后的原因吗？

阅读D.A.诺尔和 D.E.凯斯，J.Comp。物理。 193, 357 (2004)。 DOI:10.1016/j.jcp.2003.08.010，SciPy 引用的文章之一，我找到了 SciPy 选择的高级理由。

如上所示，雅可比向量乘积近似基于泰勒级数展开。在这里，我们讨论选择方程中扰动参数 ε 的各种选项。 (10)，[编者注：这是 SciPy 称为 omega 的变量，除以 v 的范数。] 给定 u 和 v，它显然对缩放敏感。如果 ε 太大，则导数的近似性很差，如果太小，有限差分的结果会受到浮点舍入误差的污染。用于单个参数的标量有限差分的最佳 ε 可以作为这两个可量化权衡的平衡进行精确优化。然而，矢量有限差分 ε 的选择既是一门艺术，也是一门科学。

因此这里需要平衡两个问题：

• 斜率的变化：如果步长很大，那么函数的雅可比可能会在

  f(x)

  和

  f(x + step)

  之间变化。
• 舍入：如果x很大，那么步长也必须很大，否则计算

  x + step

  时会出现舍入误差。

理想情况下，为了解决第一个问题，您应该查看函数的二阶导数。然而，我们不知道函数的一阶或二阶导数。这就是找到步长的重点。我认为它将

f(x)

 的大小视为下一个最好的事情：如果 

f(x)

 很大，那么要么用户在启动解算器时对 x 进行了非常错误的猜测，要么这是一个区域函数快速变化的函数。

舍入的处理方式类似，如果x很大，那么步长也会很大。这有点类似于同一篇论文中的方程（11）。

在这个方程中，n 代表问题的维数，v 代表我们试图找到雅可比向量乘积的点，u 代表乘积的方向，b 代表近似于平方的任意参数机器 epsilon 的根。

我们可以通过代数操作来找到它与 SciPy 公式之间的异同。

# difference between SciPy's omega and Knoll's epsilon
epsilon = omega / norm(v)
epsilon = 1/(n*norm(v)) * (sum(b * abs(u[i]) for i in range(n)) + b)

# Combine the two equations
omega / norm(v) = 1/(n*norm(v)) * (sum(b * abs(u[i]) for i in range(n)) + b)

# Multiply both sides by norm of v
omega = 1/n * (sum(b * abs(u[i]) for i in range(n)) + b)

# Factor out b
omega = b/n * (sum(abs(u[i]) for i in range(n)) + 1)

# Move n into parens
omega = b * (sum(abs(u[i]) for i in range(n))/n + 1/n)

# Recognize sum(...)/n as mean
omega = b * (mean(abs(u)) + 1/n)

这有点类似于 SciPy 正在做的事情，除了：

它使用平均值而不是最大值。
通过添加
• 1/n

   而不是使用 

  max(1, ...)

  ，避免采取零大小的步骤。

它不会调整
• f(x)

  。

避免无限或零步长

有两种步长必须不惜一切代价避免：零和无穷大。如果一步为零，则会导致被零除。无穷大的步长不会告诉我们任何有关局部雅可比行列式的信息。

这是一个问题，因为 x 和 f(x) 都可以为零。

max(1, ...)

 步骤很可能被放置在那里以避免这种情况。

结论

我找不到采用相同方法的现有期刊论文。我怀疑这个方程只是一种经过实验验证并在实践中有效的方法。

*注：我只读过 D.A. 的论文。诺尔和 D.E.凯斯、A.H. 贝克、E.R. 杰瑟普和 T. 曼托菲尔。该页面上的第一个参考是在选择 omega 之后添加的，所以我没有阅读它。