導関数の計算が面倒な場合の数値近似法

最急降下法やLevenberg-Marquardt法で求めたい式の導関数をイチイチ計算したり実装するのが面倒な場合に数値計算で近似する方法です。
以下の3種類くらいがよく使われるらしく、計算量と精度のバランスを考えて決めればよいようです。

• 2点近似
• 3点近似（中心差分）
• 5点近似

近似法

2点近似

一番簡単な近似で計算も一番早い。
その分精度はそこそこ。
数式で表すと

f'(x_0) \approx \frac{f(x_0+h)-f(x_0)}{h}

ただし$h$は十分小さい値です。
C言語で実装すると、

#include <float.h>

extern double f( double x );

double df( double x )
{
    const double h = FLT_EPSILON;
    double y1 = f( x + h );
    double y2 = f( x );

    return ( y1 - y2 ) / h;
}

$h$は正でも負でも計算はできるようですが結果が変わります。

FLT_EPSILONについて

FLT_EPSILONはfloat.hに定義されている値です。
微小な値で、機械イプシロンという「1より大きい最小の数」と1との差の値として定義されてます。
FLT_EPSILON自体はfloat用の機械イプシロンですがxの値が小さい場合丸め誤差が発生しそうなのでdouble用ではなくとりあえずこちらを使用してます。

定義としてはdouble用もありDBL_EPSILONが同じfloat.hに定義されてます。

3点近似（中心差分）

2点近似だと導関数の値がゼロになる点$x_0$でもゼロになりにくいのでもう少し改良して点$x_0$の前後の点を使う方法です。
計算量は2点近似よりちょっと増える分精度は向上。

数式で表すと

f'(x_0) \approx \frac{f(x_0+h)-f(x_0-h)}{2h}

C言語で実装すると、

#include <float.h>

extern double f( double x );

double df( double x )
{
    const double h = FLT_EPSILON;
    double y1 = f( x + h );
    double y2 = f( x - h );

    return ( y1 - y2 ) / ( 2 * h );
}

3点近似といいながら計算では2点しか出てこないんですね。
平均値の定理を考えるとこちらの方がxの本当の値に近い可能性が高そうです。

$h$は正負どちらでも結果は一緒になります。

5点近似

f'(x_0) \approx \frac{f(x_0-2h)-8f(x_0-h)+8f(x_0+h)-f(x_0+2h)}{12h}

一気にわけわからない式になりました。
Taylor展開とかLagrange補間多項式とかそんなのを駆使して導くようです。

C言語で実装すると、

#include <float.h>

extern double f( double x );

double df( double x )
{
    const double h = FLT_EPSILON;
    double y1 = f( x + h );
    double y2 = f( x - h );
    double y3 = f( x + 2 * h );
    double y4 = f( x - 2 * h );

    return ( y4 - 8 * y2 + 8 * y1 - y3 ) / ( 12 * h );
}

もはやなにがなんだかわからないので、こうすれば計算できることだけ覚えておいてありがたく使わせてもらうことにします。

偏微分

複数の変数でできた関数だと偏微分することになりますが、やることは同じです。
2点近似で$f(x,y,z)$の$y$を偏微分するなら

f_y(x_0,y_0,z_0) \approx \frac{f(x_0,y_0+h,z_0)-f(x_0,y_0,z_0)}{h}

これを実装することになります。