L-BFGS算法/Broyden族/BFGS算法/阻尼牛顿法的Python实现代码

下面定义了三个Python语言编写的函数：函数表达式fun，梯度向量gfun，和海森矩阵hess。这三个表达式在后面各个算法的实现中会用到。

# 函数表达式fun
fun = lambda x:100*(x[0]**2-x[1])**2 + (x[0]-1)**2

# 梯度向量 gfun
gfun = lambda x:np.array([400*x[0]*(x[0]**2-x[1])+2*(x[0]-1), -200*(x[0]**2-x[1])])

# 海森矩阵 hess
hess = lambda x:np.array([[1200*x[0]**2-400*x[1]+2, -400*x[0]],[-400*x[0],200]])

最速下降法的Python实现。

def gradient(fun,gfun,x0):
    #最速下降法求解无约束问题
    # x0是初始点，fun和gfun分别是目标函数和梯度
    maxk = 5000
    rho = 0.5
    sigma = 0.4
    k = 0
    epsilon = 1e-5

    while k<maxk:
        gk = gfun(x0)
        dk = -gk
        if np.linalg.norm(dk) < epsilon:
            break
        m = 0
        mk = 0
        while m< 20:
            if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):
                mk = m
                break
            m += 1
        x0 += rho**mk*dk
        k += 1

    return x0,fun(x0),k

4.牛顿法
牛顿法的基本思想是用迭代点xkxk处的一阶导数（梯度gkgk）和二阶倒数（海森矩阵GkGk）对目标函数进行二次函数近似，然后把二次模型的极小点作为新的迭代点。

牛顿法推导

函数f(x)f(x)在xkxk处的泰勒展开式的前三项得
T(f,xk,3)=fk+gTk(x−xk)+12(x−xk)TGk(x−xk)
T(f,xk,3)=fk+gkT(x−xk)+12(x−xk)TGk(x−xk)

则其稳定点
∇T=gk+Gk(x−xk)=0
∇T=gk+Gk(x−xk)=0

若 GkGk非奇异，那么解上面的线性方程（标记其解为 xk+1xk+1）即得到牛顿法的迭代公式
xk+1=xk−G−1kgk
xk+1=xk−Gk−1gk

即优化方向为 dk=−G−1kgkdk=−Gk−1gk
求稳定点是用到了以下公式：

考虑y=xTAxy=xTAx，根据矩阵求导法则，有
dydx=Ax+ATxd2ydx2=A+AT
dydx=Ax+ATxd2ydx2=A+AT
注意实际中dkdk的是通过求解线性方程组Gkd=−gkGkd=−gk获得的。

阻尼牛顿法及其Python实现
阻尼牛顿法采用了牛顿法的思想，唯一的差别是阻尼牛顿法并不直接采用xk+1=xk+dkxk+1=xk+dk，而是用线搜索技术获得合适的步长因子αkαk,用公式xk+1=xk+δmdkxk+1=xk+δmdk计算xk+1xk+1。

阻尼牛顿法的算法描述：

step 1: 给定终止误差0≤ϵ≪1,δ∈(0,1),σ∈(0,0.5)0≤ϵ≪1,δ∈(0,1),σ∈(0,0.5). 初始点x0∈Rnx0∈Rn. 令k←0k←0
step 2: 计算gk=∇f(xk)gk=∇f(xk). 若||gk||≤ϵ||gk||≤ϵ，停止迭代，输出x∗≈xkx∗≈xk
step 3: 计算Gk=∇2f(xk)Gk=∇2f(xk)，并求解线性方程组得到解dkdk,
Gkd=−gk
Gkd=−gk

step 4: 记 mkmk是满足下列不等式的最小非负整数 mm.
f(xk+βmdk)≤f(xk)+δβmgTkdk
f(xk+βmdk)≤f(xk)+δβmgkTdk

step 5: 令 αk=δmk,xk+1=xk+αkdk,k←k+1αk=δmk,xk+1=xk+αkdk,k←k+1，转 step 2
阻尼牛顿法的Python实现：

def dampnm(fun,gfun,hess,x0):
    # 用阻尼牛顿法求解无约束问题
    #x0是初始点，fun，gfun和hess分别是目标函数值，梯度，海森矩阵的函数
    maxk = 500
    rho = 0.55
    sigma = 0.4
    k = 0
    epsilon = 1e-5
    while k < maxk:
        gk = gfun(x0)
        Gk = hess(x0)
        dk = -1.0*np.linalg.solve(Gk,gk)
        if np.linalg.norm(dk) < epsilon:
            break
        m = 0
        mk = 0
        while m < 20:
            if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):
                mk = m
                break
            m += 1
        x0 += rho**mk*dk
        k += 1

    return x0,fun(x0),k

性能展示

作图代码：

iters = []
for i in xrange(np.shape(data)[0]):
    rt = dampnm(fun,gfun,hess,data[i])
    if rt[2] <=200:
        iters.append(rt[2])
    if i%100 == 0:
        print i,rt[2]

plt.hist(iters,bins=50)
plt.title(u'阻尼牛顿法迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})

修正牛顿法法及其Python实现
牛顿法要求目标函数的海森矩阵G(x)=∇2f(x)G(x)=∇2f(x)在每个迭代点xkxk处是正定的，否则难以保证牛顿方向dk=−G−1kgkdk=−Gk−1gk是ff在xkxk处的下降方向。为克服这一缺陷，需要对牛顿法进行修正。
下面介绍两种修正方法，分别是“牛顿-最速下