ニュートン法による極値計算

テキスト「これなら分かる最適化数学」の86ページにあるニュートン法をPythonで実装した。通常のニュートン法は関数f(x)に対し、f(x)=0なるxを求める反復計算を行うものであるが、反復の式を少し変更するだけで、極値(f'(x)=0なるx)を求めることができる。詳しい理論はテキストに譲るが、x <- x - f'(x)/f''(x)とし、収束するまで繰り返せば、極値にたどり着くことができる。

ニュートン法で極値を求めるにあたり、注意しなければならないのは、コードの更新式からも見当がつくが、途中で2回微分が0にならないようにすることである。

以下、f(x)=x^3-2*x^2+x+3に対して、極値が求まるか試してみた。なお、f(x)が極値をとるxの値は解析的に解くこともでき、x=1、1/3の2つとなる。

# -*- encoding: utf-8 -*-

'''
ニュートン法による極値の計算
2回微分が存在する場合、勾配法よりも効率よく極値を計算できる
'''

from sympy import * # 微分操作のために必要
import itertools as it # mapで複数の引数を与えたいので利用する

def newton(init, vari, func, delta=10**-6):
    '''
    funcは極値を探す対象の関数、variは微分する変数、deltaは収束判定のための定数
    返す値は見つかった極値の1つである
    '''
    df = diff(func, vari)
    d2f = diff(df, vari)
    x = init
    while True:
        xvar = x
        # 2回微分の微分係数が0にならないように注意が必要
        x = xvar - float(limit(df, vari, xvar)) / float(limit(d2f, vari, xvar))
        if abs(x - xvar) < delta:
            break
    return round(x, 5) # deltaの値によって有効桁数が変わるが、簡単のため5で決めうち

x = Symbol('x')
f = x**3 - 2*x**2 + x + 3
start_num = range(-5, 5) # 初期値を複数設定しておく
print map(newton, start_num, it.repeat(x, len(start_num)),
          it.repeat(f, len(start_num)))

最後のmapを利用することで、一度に複数の結果を実行、表示することができるようになる。mapに与える引数については、繰り返し実行できる形式にしておく必要がある。

実行結果は、[0.33333, 0.33333, 0.33333, 0.33333, 0.33333, 0.33333, 1.0, 1.0, 1.0, 1.0]となる。それぞれ初期値を-5、-4、...、4に変更しながら実行した結果である。すべて1/3と1のうち、近いほうの極値に収束していることが確認できる。