对统计学习方法附录A,B中的各种优化方法进行了实现
对着公式很容易编程实现
# f(x,y) = (x-2)*(x-2) + (y-1)*(y-1) * x*x*y*y
# f'x = 2*(x-2) + 2*x*y*y f'y = 2*(y-1) + 2*y*x*x
# f'xx = 2 + 2*y*y f'xy = 4*x*y
# f'yx = 4*x*y f'yy = 2 + 2*x*x
import numpy as np
#梯度下降
def gradient(e = 1e-6,lamb = 0.001):
x = np.array([0,0]).reshape(-1,1)
gradient = np.array([2 * (x[0][0] - 2) + 2 * x[0][0] * x[1][0] * x[1][0],2 * (x[1][0] - 1) + 2 * x[1][0] * x[0][0] * x[0][0]]).reshape(-1, 1)
while( np.sqrt( sum( [x[0]*x[0] for x in gradient] ))>e ):
x = x-lamb*gradient
gradient = np.array([2*(x[0][0]-2) + 2*x[0][0]*x[1][0]*x[1][0],2*(x[1][0]-1) + 2*x[1][0]*x[0][0]*x[0][0]]).reshape(-1,1)
print("梯度下降法:")
print(x)
print()
#牛顿法
def newton(e = 1e-1,lamb = 0.1):
x = np.array([0,0]).reshape(-1,1)
gradient = np.array([2 * (x[0][0] - 2) + 2 * x[0][0] * x[1][0] * x[1][0],2 * (x[1][0] - 1) + 2 * x[1][0] * x[0][0] * x[0][0]]).reshape(-1, 1)
H = np.array([2+2*x[1][0]*x[1][0],4*x[0][0]*x[1][0] ,4*x[1][0]*x[0][0] ,2+2*x[0][0]*x[0][0] ]).reshape(2,2)
while( np.sqrt( sum( [x[0]*x[0] for x in gradient] ) ) >e ):
print(np.sqrt( sum( [x[0]*x[0] for x in gradient] ) ))
x = x-lamb*np.dot(np.linalg.inv(H),gradient)
gradient = np.array([2*(x[0][0]-2) + 2*x[0][0]*x[1][0]*x[1][0],2*(x[1][0]-1) + 2*x[1][0]*x[0][0]*x[0][0]]).reshape(-1,1)
H = np.array([4 * x[0][0] * x[0][0], 4 * x[0][0] * x[1][0], 4 * x[1][0] * x[0][0], 4 * x[1][0] * x[1][0]]).reshape(2, 2)
print("牛顿法:")
print(x)
print()
# DFP 拟牛顿法
def dfp(e = 1e-6,lamb = 0.001):
x = np.array([0,0]).reshape(-1,1)
gradient = np.array([2 * (x[0][0] - 2) + 2 * x[0][0] * x[1][0] * x[1][0],2 * (x[1][0] - 1) + 2 * x[1][0] * x[0][0] * x[0][0]]).reshape(-1, 1)
G = np.array([1,0 ,0 ,1 ]).reshape(2,2)
while( np.sqrt( sum( [x[0]*x[0] for x in gradient] ) ) >e ):
x = x-lamb*np.dot(G,gradient)
delta_x = -lamb*np.dot(G,gradient)
yk = np.array([2*(x[0][0]-2) + 2*x[0][0]*x[1][0]*x[1][0],2*(x[1][0]-1) + 2*x[1][0]*x[0][0]*x[0][0]]).reshape(-1,1) - gradient
gradient = np.array([2*(x[0][0]-2) + 2*x[0][0]*x[1][0]*x[1][0],2*(x[1][0]-1) + 2*x[1][0]*x[0][0]*x[0][0]]).reshape(-1,1)
G = G + np.dot(delta_x,np.transpose(delta_x)) / np.dot(np.transpose(delta_x),yk) - np.dot( np.dot(G,yk),np.dot(np.transpose(yk),G) ) / np.dot(np.dot(np.transpose(yk),G),yk)
print("dfp:")
print(x)
print()
# BFGS 拟牛顿法
def bfgs(e = 1e-6,lamb = 0.001):
x = np.array([0,0]).reshape(-1,1)
gradient = np.array([2 * (x[0][0] - 2) + 2 * x[0][0] * x[1][0] * x[1][0],2 * (x[1][0] - 1) + 2 * x[1][0] * x[0][0] * x[0][0]]).reshape(-1, 1)
B = np.array([1,0 ,0 ,1 ]).reshape(2,2)
while( np.sqrt( sum( [x[0]*x[0] for x in gradient] ) ) >e ):
x = x-lamb*np.dot(np.transpose(B),gradient)
delta_x = -lamb*np.dot(np.transpose(B),gradient)
yk = np.array([2*(x[0][0]-2) + 2*x[0][0]*x[1][0]*x[1][0],2*(x[1][0]-1) + 2*x[1][0]*x[0][0]*x[0][0]]).reshape(-1,1) - gradient
gradient = np.array([2*(x[0][0]-2) + 2*x[0][0]*x[1][0]*x[1][0],2*(x[1][0]-1) + 2*x[1][0]*x[0][0]*x[0][0]]).reshape(-1,1)
B = B + np.dot(yk,np.transpose(yk)) / np.dot(np.transpose(yk),delta_x) - np.dot( np.dot(B,delta_x),np.dot(np.transpose(delta_x),B) ) / np.dot(np.dot(np.transpose(delta_x),B),delta_x)
print("bfgs:")
print(x)
print()
# byoyden类的 拟牛顿法
def broyden(e = 1e-6,lamb = 0.001,alpha=0.5):
x = np.array([0,0]).reshape(-1,1)
gradient = np.array([2 * (x[0][0] - 2) + 2 * x[0][0] * x[1][0] * x[1][0],2 * (x[1][0] - 1) + 2 * x[1][0] * x[0][0] * x[0][0]]).reshape(-1, 1)
G = np.array([1,0 ,0 ,1 ]).reshape(2,2)
while( np.sqrt( sum( [x[0]*x[0] for x in gradient] ) ) >e ):
x = x-lamb*np.dot(G,gradient)
delta_x = -lamb*np.dot(G,gradient)
yk = np.array([2*(x[0][0]-2) + 2*x[0][0]*x[1][0]*x[1][0],2*(x[1][0]-1) + 2*x[1][0]*x[0][0]*x[0][0]]).reshape(-1,1) - gradient
gradient = np.array([2*(x[0][0]-2) + 2*x[0][0]*x[1][0]*x[1][0],2*(x[1][0]-1) + 2*x[1][0]*x[0][0]*x[0][0]]).reshape(-1,1)
G_dfp = G + np.dot(delta_x,np.transpose(delta_x)) / np.dot(np.transpose(delta_x),yk) - np.dot( np.dot(G,yk),np.dot(np.transpose(yk),G) ) / np.dot(np.dot(np.transpose(yk),G),yk)
G_bfgs = np.dot ( np.dot( np.identity(2)-np.dot(delta_x,np.transpose(yk))/np.dot(np.transpose(delta_x),yk) , G ) , np.identity(2)-np.dot(delta_x,np.transpose(yk))/np.dot(np.transpose(delta_x),yk) ) + np.dot(delta_x,np.transpose(delta_x))/np.dot(np.transpose(delta_x) ,yk)
G = alpha*G_dfp +(1-alpha)*G_bfgs
print("broyden:")
print(x)
print()
if __name__ == "__main__":
gradient()
#newton() 牛顿法无法求逆,奇异
dfp()
bfgs()
broyden()