問題設定

一変数関数$y=(x+1) \times x \times (x-2)$を与えて、ニュートン法によって根を導出する。

この解は$x = -1, 0, 2$だが、はたしてこの解がニュートン法で得られるだろうか？

実装

いずれかひとつの解を求める

import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import theano
import theano.tensor as T

input = np.linspace(-5, 5, 1000)

x = T.dscalar('x')
y = (x+1)*x*(x-2)
f = theano.function(inputs=[x], outputs=y)

update = y/T.grad(cost=y, wrt=x)
init = 3.0
a = theano.shared(init)
newton = theano.function(inputs=[], outputs=a-update, givens={x: a}, updates={a: a-update})

n_train = 5
output = [newton() for i in range(n_train)]

plt.plot(input, [f(input[i]) for i in range(len(input))])
plt.scatter(init, f(init), color="red")
plt.scatter(output, [f(output[i]) for i in range(len(output))], color="green")
plt.xlim(-5, 5)
plt.ylim(-15, 15)
plt.show()

output.insert(0, np.array(init))
print(output)

[array(3.0), array(2.3684210526315788), array(2.0771631246659057), array(2.00452016332014), array(2.0000169296346986), array(2.0000000002388387)]

吟味

無事、ひとつの解(x=2)で収束するようになった。

複数の解を同時に導く

ランダムで初期値を複数生成し、それらからニュートン法を実行する。

導かれた解を判定し、指定した数と異なる解の個数が同じだったら、成功とする。

失敗したらまた初期値を乱数生成して繰り返す。

x = T.dscalar('x')
y = (x+1)*x*(x-2)
f = theano.function(inputs=[x], outputs=y)

def newton(x, y, init=3.0, n_train=10):
    update = y/T.grad(cost=y, wrt=x)
    a = theano.shared(init)
    newton = theano.function(inputs=[], outputs=a-update, givens={x: a}, updates={a: a-update})
    output = [newton() for i in range(n_train)]
    return output[-1]

def make_init(n):
    
    def r_sign(x):
        return -x if np.random.randint(2)%2==0 else x
    
    return list(map(r_sign, np.random.rand(n)))

def executor(n=3, max_loop=1000, n_train=20):

    for i in range(max_loop):
        ans = [newton(x, y, init=i, n_train=n_train).tolist() for i in make_init(n)]
        if len(set(ans)) == n:
            break
    return ans

executor()

[0.0, 2.0, -1.0]

吟味

無事、３つの解すべてを得ることが出来た。