数値勾配を計算するための「標準的な」方法はありますか？

Question

私は、C++で滑らかな関数の数値勾配を計算しようとしています。パラメータ値はゼロから非常に大きな値（おそらく1e10から1e20まで）で変化する可能性があります。

私はテストベンチとして関数f（x、y）= 10 * x^3 + y^3を使用しましたが、 xまたはyが大きすぎると、正しい勾配を得ることができないことがわかりました。ここで

はgraidientを計算するために私のコードです：

#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; 
} 

私は勾配を計算するために上記のコードを使用している場合は、勾配がゼロになります！

テストベンチ関数10 * x^3 + y^3は単なるデモですが、私が解決する必要がある実際の問題は実際にはブラックボックス関数です。

したがって、数値勾配を計算するための「標準的な」方法はありますか？

Answer 1

勾配を計算する方法は計算です。 （i、j）はxおよびy方向の単位ベクトルである

g(x, y) = Df/Dx i + Df/Dy j 

、それぞれ：

グラデーションベクトルです。おおよそのデリバティブへ

一つの方法は、一次違いです：

Df/Dx ~ (f(x2, y)-f(x1, y))/(x2-x1) 

と

何をやっているようには見えません

Df/Dy ~ (f(x, y2)-f(x, y1))/(y2-y1) 

。

あなたが閉じた形の式を有する：あなたは（X、Y）の値をプラグインし、任意の時点で正確に勾配を計算することができる

g(x, y) = 30*x^2 i + 3*y^2 j 

を。それをあなたの違いと比較し、あなたの近似がどれほどうまくいくかを見てください。

数字で実装する方法は、あなたの責任です。 （10^19）^ 3 = 10^57、そうですか？

マシン上のダブルのサイズはどのくらいですか？ 64ビットIEEE倍精度浮動小数点数ですか？

Answer 2

使用

dx = (1+abs(x))*eps, dfdx = (f(x+dx,y) - f(x,y))/dx 
dy = (1+abs(y))*eps, dfdy = (f(x,y+dy) - f(x,y))/dy 

大きな引数の意味のステップサイズを取得します。

片側差分式の場合はeps = 1e-8、中央差分差の場合はeps = 1e-5を使用してください。

差分の商と誤差が非常に小さい派生品の自動微分（autodiff.orgを参照）を調べてください。

Answer 3

必要な精度を考慮する必要があります。一見

、|y| = 5.49756e14とepsi = 1e-4ので、あなたが考慮すべきyとy+epsiため（つまりも仮数として知られている電話番号の桁数を符号化するために使用されるビットの数である）仮精度の少なくとも⌈log2(5.49756e14)-log2(1e-4)⌉ = 63ビットを必要とします異なる。

倍精度浮動小数点フォーマットは、（8バイトと仮定して）53ビットの有効となる精度を持っています。したがって、,y+epsiおよびy-epsiが等しいため、現在、f1,f2およびf3は全く同じです。

ここでは、制限：y = 1e20と、関数の結果、10x^3 + y^3について考えてみましょう。今はxを無視してみましょう。f = y^3としましょう。今度は、f(y)とf(y+epsi)が異なることになる精度を計算することができます：f(y) = 1e60とf(epsi) = 1e-12。これにより、最小有効数字精度は⌈log2(1e60)-log2(1e-12)⌉ = 240ビットになります。あなたはそれが16バイトであると仮定すると、long doubleタイプを使用していたとしても

、あなたの結果は異なっていないでしょう：f1、f2とf3はまだ等しくなり、たとえyとy+epsiはないでしょう。

xを考慮すると、fの最大値は11e60（x = y = 1e20）になります。したがって、精度の上限は⌈log2(11e60)-log2(1e-12)⌉ = 243ビット、または少なくとも31バイトです。

問題を解決する1つの方法は、別のタイプ、多分固定小数点として使用されるbignumを使用することです。

もう1つの方法は、問題を再考し、それを別の方法で処理することです。最終的に、あなたが望むのはf1 - f2です。 f(y+epsi)を分解しようとすることができます。繰り返しますが、xを無視すると、f(y+epsi) = (y+epsi)^3 = y^3 + 3*y^2*epsi + 3*y*epsi^2 + epsi^3です。だからf(y+epsi) - f(y) = 3*y^2*epsi + 3*y*epsi^2 + epsi^3。

Answer 4

最初は、より正確な（テイラー開発のもう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、相対値は動作しないことができるとあなたが戻っ定数に落下する必要がある場合ということに注意してください。しかし、何を使うべきかは何も教えてくれません！

Answer 5

以下のプログラムを使用して、微分の誤差の挙動を調べることができます。変化するステップサイズを使用して片側微分と中心差分微分を計算します。ここではxとy〜10^10を使用していますが、これは使用していたものよりも小さくなりますが、同じ点を説明する必要があります。

#include <iostream> 
#include <cmath> 
#include <cassert> 
using namespace std; 
double f(double x, double y) { 
    return 10 * pow(x, 3) + pow(y, 3); 
} 

double f_x(double x, double y) { 
    return 3 * 10 * pow(x,2); 
} 

double f_y(double x, double y) { 
    return 3 * pow(y,2); 
} 

int main() 
{ 
    // double x = -5897182590.8347721; 
    // double y = 269857217.0017581; 
    double x = 1.13041e+10; 
    double y = -5.49756e+10; 
    //double x = 10.1; 
    //double y = -5.2; 

    double epsi = 1e8; 
    for(int i=0; i<60; ++i) { 
    double dfx_n = (f(x+epsi,y) - f(x,y))/epsi; 
    double dfx_cd = (f(x+epsi,y) - f(x-epsi,y))/(2*epsi); 
    double dfx = f_x(x,y); 
    cout<<epsi<<" "<<fabs(dfx-dfx_n)<<" "<<fabs(dfx - dfx_cd)<<std::endl; 
    epsi/=1.5; 
    } 
    return 0; 
} 

出力は1両面差はおよそ100.0のステップ長で私たちについて1.37034e+13の最適なエラーを取得することを示しています。この誤差が大きく見えるが、相対誤差としては3.5746632302764072e-09と比較して

（正確な値は3.833e+21あるので）両面の差が約のステップサイズで約1.89493e+10の最適誤差を取得することに注意してください45109.3。これは3桁の大きさがあります（はるかに大きな刻み幅）。

ステップサイズはどのように調整できますか？ Yves Daostsのコメントには、

h=x_c sqrt(eps)（片面）と、h=x_c cbrt(eps)（両面）のリンクがあります。

いずれにしても、x〜10^10での正確な精度のための必要なステップサイズが100.0であれば、x〜10^20で必要なステップサイズも10^10になります。だから問題は単純にあなたのステップサイズがの方法が小さすぎることです。

これは、上記コードの開始ステップサイズを大きくし、x/y値を元の値にリセットすることで確認できます。

そして予想誘導体O(1e39)で、約O(1e31)の最良片面エラーが5.9e10のステップ長の近くに発生し、約O(1e29)の最良両面エラーが6.1e13のステップ長の近くに発生します。

Answer 6

数値微分が調整されている（小さな誤差が結果を大きく変える可能性があることを意味する）ので、Cauchy's integral formulaを使用することを検討する必要があります。このようにして、積分でn次導関数を計算することができます。これにより、精度と安定性を考慮することで問題は少なくなります。