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

昨日 に続いて、非線形最小二乗問題を解く方法として、Gauss-Newton 法と 黄金探索法 を組み合わせた手法を紹介します。

黄金探索法でステップ幅 \(\alpha\) の最適値を求める

Gauss-Newton 法は、初期値が真の解から大きく離れていたり、モデル関数が悪いと解が収束しないことがあります。 その場合、まずは初期値とモデルを考え直しますが、より安定して解が得られる方法を使ってみたりします。

その方法はいくつか提案されていますが、Newton 法の場合と同じでステップ幅 \(\alpha (0 \lt \alpha \leq {1})\) の最適値を黄金探索法で求める方法なんかもよく使われます。 しかし、計算量が増えて遅くなるので、通常の Gauss-Newton 法で十分なときは使う必要はありません。 Newton 法で黄金探索法を使用した場合 と見比べてみると共通性が見えて理解しやすいかと思います。

import numpy as np
import matplotlib.pyplot as plt


def goldenRatioSearch(function, range, tolerance):
    # 黄金探索法によるステップ幅の最適化
    gamma = (-1+np.sqrt(5))/2
    a = range[0]
    b = range[1]
    p = b-gamma*(b-a)
    q = a+gamma*(b-a)
    Fp = function(p)
    Fq = function(q)
    width = 1e8
    while width > tolerance:
        if Fp <= Fq:
            b = q
            q = p
            Fq = Fp
            p = b-gamma*(b-a)
            Fp = function(p)
        else:
            a = p
            p = q
            Fp = Fq
            q = a+gamma*(b-a)
            Fq = function(q)
        width = abs(b-a)/2
    alpha = (a+b)/2
    return alpha


def gaussNewton(function, beta, tolerance, epsilon, flag):
    # Gauss-Newton 法
    delta = 2*tolerance
    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)  # 探索方向
        if flag == 1:
            alpha = goldenRatioSearch(lambda alpha: np.linalg.norm(function(beta+alpha*delta)), [1e-2, 1], 1e-2)
        else:
            alpha = 1
        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, 1)

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

結果は以下のようになり変わりません。

>>> [ 0.36183687  0.55626645]

コメント

Comments powered by Disqus
書籍更新情報
2017-02-18
Pythonによる科学技術計算 基礎編
1.4版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
電子書籍
線形代数(1): Pythonによる科学技術計算 実践編