無約束演算法-最速下降，牛頓法，擬牛頓，共軛梯度求解二次函式極小值

阿新 • • 發佈：2018-12-21

import numpy as np
import math
a=np.random.randint(1,10,size=[100,1])
G=(a*a.T)+np.random.randint(1,5)*np.eye(100)
b=np.dot(G,np.ones((100,1)))
x0=np.random.random(size=[100])
def Df(x):#求出梯度函式
    return (np.dot(x,G)+b.T).reshape(100)
def hessian(df,z):#求hessian矩陣並賦值
    return np.array(G)
def fun(x):
    return 0.5*np.dot(np.dot(x,G),x.T)+np.dot(x,b)
def consice(x,dk):#精確線搜尋
    x=np.array(x)
    return (-np.dot(Df(x),dk)/(np.dot(np.dot(dk,G),dk)))
def stepest(x):#最速下降法
    return -Df(x)
def basicnewton(x):#牛頓法
    A=np.array(-np.dot(np.mat(G).I,Df(x)))
    A=A.reshape(100)
    return A
def Nnewton(x0,H0):#擬牛頓法
    return np.array(-np.dot(H0,Df(x0)))
def SR1(H0,sk,yk):#對稱秩1修正公式
    return np.array(H0+((sk.T-np.dot(H0,yk.T))*(sk-np.dot(H0,yk.T).T))/np.array((np.dot(sk-np.dot(H0,yk.T).T,yk.T))).reshape(1)[0])
def DFP(H0,sk,yk):
    return np.array(H0+(sk.T*sk)/(np.array(yk*sk.T).reshape(1)[0])-np.dot(np.dot(H0,yk.T)*yk,H0)/np.array(np.dot(np.dot(yk,H0),yk.T)).reshape(1)[0])
def BFGS(H0,sk,yk):
    return np.array(H0+(yk.T*yk)/(np.array(sk*yk.T).reshape(1)[0])-np.dot(np.dot(H0,sk.T)*sk,H0)/np.array(np.dot(np.dot(sk,H0),sk.T)).reshape(1)[0])
def solve(x0,eps):#x0為一n維向量
    k=0
    H0=np.eye(100)
    dk=-Df(x0)
    while np.linalg.norm(Df(x0))>eps:
        if 0:
            dk=basicnewton(x0)#最速下降法：stepest，牛頓法：basicnewton
            ak=consice(x0,dk)#精確線搜尋：consice，Armijo非精確線搜尋：Armijo，Wolfe非精確線搜尋：Wolfe
            x0=x0+np.dot(ak,dk)
        if 0:#擬牛頓法
            dk=Nnewton(x0,H0)#擬牛頓法
            ak=consice(x0,dk)#精確線搜尋：consice，Armijo非精確線搜尋：Armijo，Wolfe非精確線搜尋：Wolfe
            x1=x0
            y1=Df(x0)
            x0=x0+np.dot(ak,dk)
            sk=np.mat(x0-x1)
            yk=np.mat(Df(x0)-y1)
            H0=BFGS(H0,sk,yk)#修正公式選擇可選為SR1,DFP，BFGS
        if 1:#共軛梯度法
            ak=consice(x0,dk)#精確線搜尋：consice，Armijo非精確線搜尋：Armijo，Wolfe非精確線搜尋：Wolfe
            y1=Df(x0)
            x0=x0+np.dot(ak,dk)
            bk=np.dot(Df(x0),Df(x0))/np.dot(y1,y1)
            dk=-Df(x0)+bk*dk
        k+=1
        print("當前迭代點x"+str(k)+"="+str(x0)+",步長ak="+str(ak)+"，方向dk="+str(dk)+"梯度範數為"+str(np.linalg.norm(Df(x0))))
    return x0
solve(x0,0.01)#初始搜尋點，精度

引數說明：

G=[[ 5., 1., 7., ..., 5., 2., 2.],

[ 1., 5., 7., ..., 5., 2., 2.],

[ 7., 7., 53., ..., 35., 14., 14.],

...,

[ 5., 5., 35., ..., 29., 10., 10.],

[ 2., 2., 14., ..., 10., 8., 4.],

[ 2., 2., 14., ..., 10., 4., 8.]]

b=[ 518., 518., 3602., 3602., 3088., 4116., 2060., 2060., 2060., 2060., 3602., 1546., 2574., 2060., 1546., 3602., 4116., 518., 3602., 4116., 3602., 3088., 3088., 2060., 2574., 3088., 2574., 4630., 1546., 1032., 1032., 1546., 4630., 2574., 2574., 1546., 4116., 2574., 1032., 4630., 3602., 4116., 1546., 3088., 518., 1032., 4630., 4116., 2574., 3088., 4116., 3602., 518., 2574., 2574., 518., 1546., 1546., 1032., 4630., 3602., 2574., 4630., 518., 4630., 3088., 3088., 2574., 3088., 2574., 518., 3088., 3602., 3602., 3602., 3088., 2574., 2574., 4116., 1032., 2574., 2060., 3088., 4116., 1032., 3602., 518., 4116., 3602., 4116., 1032., 1032., 4116., 2574., 2060., 2060., 4630., 2574., 1032., 1032.]

初始點x0：[0.36287532, 0.56325161, 0.51364076, 0.56812451, 0.7888315 , 0.51851197, 0.71170918, 0.21314076, 0.29913342, 0.60791091, 0.66803629, 0.9909708 , 0.25767777, 0.55595337, 0.38623466, 0.62318896, 0.33560756, 0.46152935, 0.2218254 , 0.25434027, 0.98353889, 0.66276815, 0.7908568 , 0.61753831, 0.32039222, 0.71970002, 0.60152159, 0.57305435, 0.79107977, 0.94360134, 0.9204414 , 0.68385586, 0.15941899, 0.56939455, 0.70809533, 0.71181405, 0.44015248, 0.64435278, 0.97920569, 0.79014869, 0.43985713, 0.74625209, 0.40611885, 0.11572177, 0.13629511, 0.4138309 , 0.36313254, 0.83684369, 0.59573989, 0.31371232, 0.12146105, 0.95490561, 0.46864787, 0.06499357, 0.44101515, 0.20665268, 0.35984193, 0.92585616, 0.15035858, 0.02689437, 0.58740769, 0.9081886 , 0.55210143, 0.96850184, 0.86565383, 0.60416907, 0.26910545, 0.18905435, 0.04546626, 0.80063242, 0.78408039, 0.14940875, 0.11748599, 0.46092049, 0.63143222, 0.75883652, 0.98800527, 0.68956898, 0.52701569, 0.62221375, 0.40302016, 0.75470691, 0.53014519, 0.37138053, 0.43437477, 0.65774449, 0.72548463, 0.3557923 , 0.99207761, 0.9516513 , 0.32944806, 0.77629668, 0.11325534, 0.83355515, 0.66354271, 0.19271679, 0.29166004, 0.72486387, 0.05547915, 0.42431452]

精度=0.01

最速下降法：

執行結果：

第一次迭代資訊如下：

第二次迭代資訊如下：

第三次迭代資訊如下：

在x3處停止運算，經過三次迭代後找到極小值點

牛頓法：

第一次迭代資訊如下：

在x1處停止運算，經過一次迭代後找到極小值點

擬牛頓DFGS方法：

由於迭代次數過多，故只列出最後一次迭代資訊：

在x965處停止計算，經過965次迭代後找到極小值點

共軛梯度法

第一次迭代資訊如下：

第二次迭代資訊如下：

在x2處停止計算，經過兩次迭代找到極小值點

方法對比

最速下降法：3次迭代

牛頓法：1次迭代（精確解）

擬牛頓BFGS方法：965次迭代

共軛梯度法：2次迭代（精確解）