ニュートン法による最適化
ある関数$${f(x)}$$について、$${f(x) = 0}$$となるような$${x}$$を見つけたいとします。エラー関数がゼロになるような変数の値を知りたい、ということです。$${x}$$はベクトルであってもかまいません。
関数の$${x = a}$$でのテーラー展開が得られるものとすると
$$
f(x) = \displaystyle\sum_{0}^{\infty} \frac{f^{(n)}(a)}{n!}(x - a)^n …(1)
$$
(1)の１階微分までを用いると
$$
f(x) \fallingdotseq \frac{f(a)}{0!}(x - a)^0 + \frac{f^{(1)}(a)}{1!}(x - a)^1
$$
$$
f(x) \fallingdotseq f(a) + f'(a)(x - a) …(2)
$$
$${f(x)=0}$$となる$${x}$$を求めたいので、(2)は
$$
0 = f(a) + f'(a)(x - a)
$$
よって
$$
x = a - \frac{f(a)}{f'(a)} …(3)
$$
となります。
この(3)は、$${x=a}$$のとき、次に$${x}$$をいくらに修正すれば$${f(x) = 0}$$になるのか（近づくのか）を推定していると考えればよいので、$${a \to x_i , x \to x_{i+1}}$$と置き換えると
$$
x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)} …(4)
$$
$${x_i から x_{i+1}}$$を推定するときに、補正が大きくかかりすぎないようにするテクニックとして、係数$${\mu=0 ～ 1}$$を用いて
$$
x_{i+1} = x_i - \mu \frac{f(x_i)}{f'(x_i)} …(5)
$$
最終的に$${f(x)=0}$$となる$${x}$$を求めるということは、途中経過では
$$
|f(x_{i+1})|<|f(x_i)| …(6)
$$
であればよいことになります。ある$${x_i}$$のとき、$${\mu = 1}$$から始めて$${x_{i+1}}$$を求め、(6)が満たされるまで $${\mu = \frac{\mu}{2}}$$として小さくし繰り返します。
$${x_{i+1}}$$が決まれば、それを$${x_i}$$に入れて再度(6)を満たすように(5)を求めます。
以上を繰り返し、$${f(x)=0}$$となれば、その時の$${x}$$が求める値です。あるいは、あらかじめ決めておいた繰り返し回数に達したら、求まらなかったとして計算を止めます。
$${x}$$はスカラーである必要はなく、ベクトルであってもそれぞれ要素ごとに計算していくことで最適値の推定ができます。その場合の微分は偏微分です。関数の微分の計算方法については、別の記事に書いていますのでご覧ください（数値微分の検討、微小hの値について）。
例として$${x}$$が２次元ベクトルのエラー関数の場合を示します。
エラー関数は、$${f(\.{x}) = (x_0 - 1)^2 + (x_1 - 2)^2}$$ とします。$${f(\.{x}) = 0}$$となる最適値は、$${\.{x} = (1, 2)}$$です（Fig.1）。当然ですが、実際の問題の時は、エラー関数の形も最適値もわかっているわけではありません。エラー関数は$${\.{x}}$$が与えられれば一意に値が求まればよいです。
Excel VBAのソースコードを示します。これで計算した結果は、$${x_0 = 1.00000442461132, x_1= 1.9999950680878 }$$となりました。精度が十分かどうかは問題の性質によりますが、適切な結果が得られていると思います。
'ニュートン法のテスト
Sub testNewtonN()
Const TargetValue As Double = 0.0000000001 '狙いのエラー値
Const iMAX As Long = 1000 '試行回数
Const muMin As Double = 0.0000000001 'muの許容最小値
Dim mu As Double
Dim x1(0 To 1) As Double 'x初期値
Dim x2(0 To 1) As Double 'x推定値
Dim fx1 As Double 'x1での関数fxの値
Dim fx2 As Double 'x2での関数fxの値
'xの初期値
x1(0) = 10
x1(1) = -10
fx1 = funcN(x1) 'その時の関数fxの値
Dim df() As Double '微分ベクトル
Dim k As Long 'ベクトルの要素番号
Dim i As Long '試行回数
i = 0
Do
df = diffN(x1) '微分ベクトルを求める
fx1 = funcN(x1) '今のエラー関数値
For k = 0 To UBound(x1) 'ベクトルxの要素ごとに
mu = 1
Do
x2(k) = x1(k) - mu * fx1 / df(k) 'xの要素の最適推定値を求める
fx2 = funcN(x2) 'xの要素の最適推定値でのエラー関数値
If Abs(fx2) < Abs(fx1) Then 'エラーは小さくなった
x1(k) = x2(k) 'その要素の値を採用する
Exit Do
Else
mu = mu / 2 'エラーは小さくないのでmuを小さくする
End If
If mu < muMin Then 'muは許容範囲か
Exit Do
End If
Loop
Next
If Abs(fx2) < TargetValue Then
Exit Do
End If
If i >= iMAX Then
Exit Do
End If
i = i + 1
Loop
Debug.Print x2(0), x2(1)
End Sub'f(x)=(x(0)-1)^2 + (x(1)-2)^2
Function funcN(x() As Double) As Double
Dim res As Double
Dim i As Long
res = 0
For i = 0 To UBound(x)
res = res + (x(i) - (i + 1)) ^ 2
Next
funcN = res
End Function'偏微分を求める関数
Function diffN(x() As Double) As Double()
Const Accuracy = 0.000001 '微分の精度
Const iMAX = 1000 '最大試行回数
Dim df() As Double '微分の結果
ReDim df(UBound(x))
Dim dfpost As Double '直前の試行の微分の値
Dim h As Double '微小変化
Dim i As Long '繰り返しのカウント
Dim x2() As Double
ReDim x2(UBound(x))
Dim k As Long 'ベクトルxの要素番号
For k = 0 To UBound(x)
i = 0
h = 1
x2 = x
x2(k) = x2(k) + h
dfpost = (funcN(x2) - funcN(x)) / h
h = h / 2
Do
x2 = x
x2(k) = x2(k) + h
df(k) = (funcN(x2) - funcN(x)) / h
If Abs(df(k) - dfpost) < Accuracy Then
Exit Do
End If
h = h / 2
dfpost = df(k)
i = i + 1
Loop While i < iMAX
'注意。この状態ではi=iMaxに達して
'微分が正しくも止まっていなくても
'そのままになっている
Next
diffN = df
End Function応援してやろうということで、お気持ちをいただければ嬉しいです。もっと勉強したり、調べたりする糧にしたいと思います。