是否有任何“标准”方式来计算数字梯度？

Question

我正在尝试计算C ++中平滑函数的数值梯度。并且参数值可以从零到很大的数字（可能是1e10到1e20？）变化。

我使用函数f（x，y）= 10 * x ^ 3 + y ^ 3作为测试平台，但是我发现如果x或y太大，我将无法获得正确的渐变。

这是我的代码来计算应付款：

#include <iostream>
#include <cmath>
#include <cassert>
using namespace std;
double f(double x, double y)
{
    // black box expensive function
    return 10 * pow(x, 3) + pow(y, 3);
}
int main()
{
    // double x = -5897182590.8347721;
    // double y = 269857217.0017581;
    double x = 1.13041e+19;
    double y = -5.49756e+14;
    const double epsi = 1e-4;

    double f1 = f(x, y);
    double f2 = f(x, y+epsi);
    double f3 = f(x, y-epsi);
    cout << f1 << endl;
    cout << f2 << endl;
    cout << f3 << endl;
    cout << f1 - f2 << endl; // 0
    cout << f2 - f3 << endl; // 0
    return 0;
}

如果我使用上面的代码来计算梯度，则梯度将为零！

testbench函数10 * x ^ 3 + y ^ 3只是一个演示，我需要解决的真正问题实际上是黑盒函数。

那么，有没有“标准”的方法来计算数值梯度呢？

Answer 1

首先，您应该使用中央差分方案，该方案更加准确（通过取消泰勒发展的另一个任期）。

(f(x + h) - f(x - h)) / 2h

而不是

(f(x + h) - f(x)) / h

然后，选择h至关重要，并且使用固定常数是您最糟糕的事情。因为对于x，h将太大，以至于逼近公式不再起作用；对于x，h将变得太小，从而导致严重的截断误差。

一个更好的选择是采取的相对值，h = x??其中?是机器的ε-（1个ULP），其给出了良好的折衷。

(f(x(1 + ??)) - f(x(1 - ??))) / 2x??

注意，当时x = 0，相对值不起作用，您需要退回给常数。但是，什么也没有告诉您使用哪一个！