统计学习方法之梯度下降，牛顿法，拟牛顿法

对统计学习方法附录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()