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]