ヘッシアン(Hessian)

n変数関数のヘッシアンを求める関数を書きました.

# f(x1, ..., xn)のxiによる偏導関数を返す
def partial(f, n, i, dx=0.001):
    I = np.identity(n)
    def part(x):
        return (f(x+I[i]*dx) - f(x)) / dx
        #return (f(x+I[i]*2*dx) - f(x)) / 2*dx
    return part

# f(x1, ..., xn)のHessianを返す
def Hessian(f, n, dx=0.01):
    def h(x):
        return np.array([np.array([H[i][j](x) for j in np.arange(n)], dtype=float) for i in np.arange(n)])
    H = [[partial(partial(f, n, i), n, j) for j in np.arange(n)] for i in np.arange(n)]
    return h

コメントに書いてある通りでございます. n変数関数fのある点\bf{x}でのヘッシアンの値を求めたかったら,以下のようにすると求まります.

H = Hessian(f, n)
print H(x)

明示的にヘッシアンを与えなくても良い代わりに当然誤差というものがついてきます. 関数の極限を計算機で扱う宿命でしょうか. 極限をうまく計算できるアルゴリズムはどんなのがあるのかな.

numpy.gradientってのがある.
僕のソースでは,fからf'を求めるイメージだけど, numpy.gradientはf(x)からf'(x)を求めることができます.
僕:関数→導関数
numpy.gradient:fの値の列→f’の値の列
という感じかな.

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中