Pythonで勾配法を実装してみました。

今回はこの本に出てくる勾配法を実装します。
www.amazon.co.jp

1変数の場合

第3章のアルゴリズム3.1(1変数の勾配法)を、Pythonで実装してみました。gradient_descentは比較のために書いています。本に登場するのはgradient_descent_2の方です。

import numpy as np

def f(x):
    return - x ** 2

def f_prime(x):
    return  - 2 * x

def gradient_descent(x, step, prime, iteration, e):
    for i in range(iteration):
        x_new = x + prime(x) * step
        if abs(f(x_new) - f(x)) < e:
            break
        x = x_new
    return(x,i)

def gradient_descent_2(x, h, prime, iteration, e):
    for i in range(iteration):
        h = np.sign(prime(x)) * abs(h)
        X = x
        X_new = x + h
        
        if f(X) < f(X_new):
            while f(X) < f(X_new):
                if f(X) >= f(X + h):
                    break
                h = 2 * h
                X = X_new
                X_new = X + h
            x = X
            h = h * 0.5
        else:
            while f(X) > f(X_new):
                if f(X) <= f(X - h):
                    break
                h = h * 0.5
                X_new = X_new - h
            x = X_new
            h = 2 * h
        
        if abs(prime(x)) <= e:
            break
            
    return(x,i)

print( gradient_descent(-100,0.1,f_prime,1000,0.0001) )
print( gradient_descent_2(-100,0.1,f_prime,1000,0.0001) )

多変数の場合

次に、多変数の場合で実装しました。下記は2変数の場合に対応しています。この記事(Rで勾配法(勾配降下法・最急降下法)を実装 - Qiita)をかなり参考にさせて頂きました。表記なども基本的には合わせています。難易度はやはり1変数の場合と比べると難しいですね。教科書ではアルゴリズム3.2(多変数の勾配法)です。

import numpy as np

def f(x):
	return -x[0]**2 - 4 * x[1] **2 + 5 * x[0]

def dx(x):
	return -2*x[0] + 5

def dy(x):
	return -8 * x[1]

def search(x0,delta,f,df):
	x = x0
	h = 0.5
	iteration = 10000
	eps = 0.000001
	for i in range(iteration):
		oldx = x
		h = np.sign(df(x)) * abs(h)
		newx = x + h * delta
		if f(x) < f(newx):
			while f(x) < f(newx):
				if f(newx) >= f(newx+h*delta):
					break
				h = 2 * h
				x = newx
				newx = x + h * delta
			h = h * 0.5
		else:
			while f(x) > f(newx):
				if f(x) <= f(x + h * delta):
					break
				h = h * 0.5
				newx = x + h * delta
			x = newx
			h = 2 * h
		if abs(h) < eps:
			break
	return x

def hillclimb(x,f,dx,dy):
	eps = 0.00000001
	iteration = 10000
	for i in range(iteration):
		gradientf0 = np.array([dx(x),dy(x)])
		print('hillclimb i=', i, 'f(x)=', f(x), 'gradientf0=', gradientf0)
		def df(x):
			return np.dot( np.array([dx(x),dy(x)])  ,  gradientf0 )
		print df(x)
		delta = ( gradientf0 ) / np.linalg.norm(gradientf0)
		newx = search(x,delta,f,df)
		diff = newx -x
		x = newx
		if sum(diff * diff) < eps:
			break
		print(delta, sum(diff*diff))
	return (x)


print(hillclimb(np.array([1,0.5]),f,dx,dy))