Python NumPy サンプルコード: Gauss-Newton 法で非線形最小二乗問題を解く

今回は、非線形最小二乗問題を解く方法のうち、解の収束が早くて有用な Gauss-Newton 法 を紹介します。

[追記 15-06-25] こちら に SciPy の非線形最小二乗問題を解く関数についてまとめました。SciPy には Gauss-Newton 法のルーチンはありませんが、Levenberg-Marquardt 法が使えます。 本来、非線形性が強くない問題では Gauss-Newton 法の方が Levenberg-Marquardt 法より速いのですが、SciPy だと Fortran のライブラリを使っているので、Python で Gauss-Newton 法を実装するよりも高速に最適解を求められます。

Gauss-Newton 法

Gauss-Newton 法は非線形最小二乗問題を解く方法の一つで、計測データに対して理論式の未知パラメータを同定するのに使います。 Python なら以下のように書けます。Wikipedia に載っている例題をやってみました。式と見比べてみると Gauss-Newton 法が理解しやすいと思います。

# coding: utf-8
import numpy as np
import matplotlib.pyplot as plt


def gaussNewton(function, beta, tolerance, epsilon):
    # Gauss-Newton 法
    delta = 2*tolerance
    alpha = 1
    while np.linalg.norm(delta) > tolerance:
        F = function(beta)
        J = np.zeros((len(F), len(beta)))  # 有限差分ヤコビアン
        for jj in range(0, len(beta)):
            dBeta = np.zeros(beta.shape)
            dBeta[jj] = epsilon
            J[:, jj] = (function(beta+dBeta)-F)/epsilon
        delta = -np.linalg.pinv(J).dot(F)  # 探索方向
        beta = beta + alpha*delta
    return beta


def objectiveFunction(beta):
    # 目的関数
    r = y - theoreticalValue(beta)
    return r


def theoreticalValue(beta):
    # 理論値
    f = beta[0]*x / (beta[1]+x)
    return f

# サンプルデータの作成
x = np.array([0.038, 0.194, 0.425, 0.626, 1.253, 2.500, 3.740])
y = np.array([0.050, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317])

# 初期値を設定して Gauss-Newton 法で非線形最小二乗法を解く
initialValue = np.array([0.9, 0.2])
betaID = gaussNewton(objectiveFunction, initialValue, 1e-4, 1e-4)

# 結果を表示
print(betaID)
plt.figure()
plt.plot(x, y, 'r.')
plt.plot(x, theoreticalValue(betaID), 'b')
plt.show()

結果は以下のようになります。

>>> [ 0.36183738  0.55626934]

コメント

Comments powered by Disqus
書籍更新情報
2017-04-18
Pythonによる科学技術計算 基礎編
固定版:1.5版、リフロー版:1.2版への更新が可能になりました。
サポートページはこちら
電子書籍