Introduction to Optimization(四): 拟牛顿法

本节介绍：

• hessian matrix 近似
• DFP算法
• bfgs算法

hessian matrix 近似

牛顿法的基本思路是用二次函数来局部逼近目标函数 f 并解近似函数的极小点作为下一个迭代点,迭代公式

xxk+1=xxk−αFF−1(xxk)gk

但是牛顿法的缺陷是需要计算hessian矩阵，及其逆矩阵，但是实际上我们如果能保证f(xk+1)<f(xk), 就ok了，即我们想找到一个类似F−1 的矩阵Hk,使得迭代点

xxk+1=xxk−αHHkgk满足f(xk+1)<f(xk)

在xxk+1 出对f(xk+1) 一阶泰勒展开

f(xxk+1)=f(xk)+gTk(xk+1−xk)+o(||xk+1−xk||)=f(xk)−αgTkHkgk+o(||Hkgk||α)

为了保证f(xk+1)<f(xk),必须有
Hk正定

拟牛顿法

拟牛顿法的迭代思路是

dkαkxk+1=−Hkgk=argminf(xk+αdk)=xk+αdk

第二条件

对于正定矩阵H0,H1,…,Hn当目标函数为n 维二次型时，他们必须满足

Hk+1Δgi=Δxi,0≤i≤k

并且当它们满足以上条件的时候，可以证明(P135 定理 11.1) d0,d1,…,di,…,dn−1 是Q共轭的，即将拟牛顿法用到n 维二次型时只需要n 步迭代即可到达最优解

由以上可以看到，Hk 并不是唯一确定的因此，我们需要找到一种方式让Hk满足条件并且容易计算
下面讲两种算法，分别是DFP和BFGS

DFP算法

1. 选择初始点x0,对称正定矩阵H0
2. 计算 gk, 如果 ||gk||<ϵ 停止迭代
3. 计算 dk=−Hkgk
   计算αk=argminα≥0f(xk+αdk)
   xk+1=xk+αkdk
   Δxk=αkdk
   Δgk=gk+1−gk
   tmp=HkΔgk
   Hk+1=Hk+ΔxkΔxTkΔxTkΔgk−tmptmpTΔgTkHkΔgk

这里得证明这样构造的 Hk 满足两个条件

定理11.3
用DFP算法求解 二次型问题时，Hessian matrix Q=QT,有,Hk+1Δgi=Δxi,0≤i≤k
因证明过程较为难写，故这里直接给出书上证明:

定理11.4

假定 gk≠0,只要有 Hk正定 那么Hk+1 正定

书上说，DFP 算法处理某些较大的非二次型问题是容易被’卡’ 住，所以后面有四个大牛提出了一种BFGS算法解决这一问题

BFGS

BFGS取自4个人的名字
图片来源于码农场

前面提到的Hk满足的第二个条件即Hk+1Δgi=Δxi,0≤i≤k 是根据 Δgi=QΔxi, 构造的Q−1 的近似矩阵Hk, 同样，我们可以构造Q 的近似矩阵Bk+1, 那么Bk+1 满足

Δgi=Bk+1Δxi

利用互补的概念，应用DFP算法中的更新公式我们可以得到Bk的更新公式

Bk+1=Bk+ΔgkΔgTkΔgTkΔxk−(BkΔxk)(BkΔxk)TΔxTkBkΔxk

因此 Hk+1=B−1k+1

关于这个式子我们有如下引理

引理11.1
如果A 非奇异， 列向量u,v，s.t,1+vTA−1u≠0 , 那么其逆矩阵存在且

(A+uvT)−1=A−1−(A−1u)(vTA−1)1+uTA−1v

直接计算可证

因此应用两次引理，可得

HBFGSK+1=Hk+(1+ΔgTkHkΔgkΔgTkΔxk)ΔxkΔxTkΔxTkΔgk−HkΔgkΔxTk+(HkΔgkΔxTk)TΔxTkΔgk

BFGS在实践中较为常用scipy.optimize 里面也有实现.以下是我的实现

代码实现

def dfp(fun, grad, x0, args=(), g_args=(), tol=1e-8, max_iter=5000):    h0 = np.eye(len(x0))    g_0 = grad(*((x0,) + g_args))    alpha = lambda a, x, d: fun(*((x + a * d,) + args))    for i in range(max_iter):        if is_stop(g_0, np.zeros(g_0.shape), tol):            break        d = -h0.dot(g_0)        alp = minimize_scalar(alpha, bounds=(0, 10000), args=(x0, d), method='brent', tol=1e-4)        alp = alp.x        x_next = x0 + alp * d        delta_x = (alp * d).reshape((len(x0), 1))        g_next = grad(*((x_next,) + g_args))        delta_g = g_next - g_0        delta_g = delta_g.reshape((delta_x.shape))        tmp = h0.dot(delta_g)        h0 = h0 + delta_x.dot(delta_x.T) / (delta_x.T.dot(delta_g)) - tmp.dot(tmp.T) / (delta_g.T.dot(tmp))        x0 = x_next        g_0 = g_next    return OptimizeResult({'nit': i, 'x': x0, 'jac': g_0, 'fun': fun(*((x0,) + args))})def bfgs(fun, grad, x0, args=(), g_args=(), tol=1e-8, max_iter=5000):    h0 = np.eye(len(x0))    g_0 = grad(*((x0,) + g_args))    alpha = lambda a, x, d: fun(*((x + a * d,) + args))    for i in range(max_iter):        if is_stop(g_0, np.zeros(g_0.shape), tol):            break        d = -h0.dot(g_0)        alp = minimize_scalar(alpha, bounds=(0, 10000), args=(x0, d), method='brent', tol=1e-4)        alp = alp.x        x_next = x0 + alp * d        delta_x = (alp * d).reshape((len(x0), 1))        g_next = grad(*((x_next,) + g_args))        delta_g = g_next - g_0        delta_g = delta_g.reshape((delta_x.shape))        tmp = h0.dot(delta_g).dot(delta_x.T)        tmp1 = (1+delta_g.T.dot(h0).dot(delta_g)/(delta_g.T.dot(delta_x)))*delta_x.dot(delta_x.T)/(delta_x.T.dot(delta_g))        tmp2 = (tmp + tmp.T)/(delta_g.T.dot(delta_x))        h0 = h0 + tmp1 - tmp2        x0 = x_next        g_0 = g_next    return OptimizeResult({'nit': i, 'x': x0, 'jac': g_0, 'fun': fun(*((x0,) + args))})

DFP 与bfgs除了更新公式以外其他几乎完全一样

以下是在rosen 函数下的实验效果

bfgs result fun: 6.0875493993460361e-23 jac: array([ -2.87637469e-10,   1.46505030e-10]) nit: 17   x: array([ 1.,  1.])dfp result fun: 2.4549731009788525e-23 jac: array([ -1.82560633e-10,   9.29922805e-11]) nit: 17   x: array([ 1.,  1.])

两次均优于共轭梯度算法.

实际使用

在实际使用中我们用的最多的还是SGD这样的用梯度来优化的算法，因为可以看到二阶项的一个重要的弊端就是空间复杂度太高，需要存贮二阶项，这对于n维的数据来说，需要n2 的内存

本文链接http://blog.csdn.net/Dylan_Frank/article/details/78287680