牛顿法与拟牛顿法

分手后的思念是犯贱 2021-09-29 10:08 575阅读 0赞

# 牛顿法与拟牛顿法 #

优化问题是机器学习中非常重要的部分,无论应用何种算法,构建何种模型,最终我们的目的都是找到最优解的. 那优化算法是无法回避的. 当今机器学习,特别是深度学习中, 梯度下降算法(gradient descent algorithm) 可谓炙手可热. 不过优化算法不只其一种,其他算法也很常用,也很优秀,比如今天介绍的牛顿法(Newton methods)与拟牛顿法(Quasi Newton methods).

关于此算法的文献很多, 这里算是一个分享,也算是对自己知识的复习. 本文主要参考 peghoty 的博文.

考虑如下的约束问题:  
\\\[ a^\* = \\arg \\max\_\{a\} f( a) \\\]  
其中 \$f(a): \\mathbb R^N \\rightarrow \\mathbb R\$ 为优化函数, 而 a 是待优化参数向量 $ a = (a\_1,a\_2,\\dots,a\_N)^T \\in \\mathbb\{R\}^N$

(Ps: 这里回避了 用 x 表示待优参数, 因为x 在机器学习中一般指的是数据,而 实际问题中, 待优化的几乎都是模型参数而非数据). 这里要求 f 是**凸的**, 且是**二阶连续可微的**(否则后面用到的二阶偏导无法满足).

## 牛顿法 ##

牛顿法的基本思想是 用**泰勒展开式来替代原有函数**: 在现有极小值点附近对 \$f(x)\$ 二阶展开(忽略高阶项,因此是一种近似). 这样的做的好处是可以很容易找到优化条件,并且可以用到二阶偏导条件:  
\\\[ \\psi(a) = f(a^\*) + \\nabla\_\{a^\*\} f \\cdot(a - a^\*) + \\frac\{1\}\{2\}(a - a^\*)^T\\cdot \\nabla^2\_\{a^\*\} f\\cdot (a - a^\*) \\\]  
其中\$\\psi(a)\$ 是原始待优化函数; \$a^\*\$ 是当前极小值点(估计值). 由极值条件:  
\\\[ \\nabla\_a\\psi = 0 \\\]  
即:  
\\\[ \\nabla\_\{a^\*\}f + (a - a^\*)^T\\cdot \\nabla\_\{a^\*\}^2 f = 0 \\\\ \\Longrightarrow a^T= \{a^\*\}^T - ( (\\nabla\_\{a^\*\}^2f )^\{-1\})^T \\cdot \\nabla\_\{a^\*\} f^T \\\]  
可见这里要求 二阶偏导矩阵(hessian 矩阵可逆(非奇异), 幸好在实际中这个经常满足).

为简化, 梯度向量(gradient vector) 与海森矩阵(hessian matrix) 分别用 \$g\_k =\\nabla\_\{a^\*\} f^T , H\_k = ( (\\nabla\_\{a^\*\}^2f )^\{-1\})^T\$ 表示.则上式可写成:  
\\\[ a = a^\* - H\_k \\cdot g\_k \\\]  
写成迭代形式:  
\\\[ a\_\{k + 1\} = a\_k - H\_k^\{-1\} \\cdot g\_k \\\]  
这里的 \$a\_k\$ 表示第 k 步极值. 一般 称上式的搜索方向\$d\_k = - g\_k \\cdot H\_k^\{-1\}\$ 为牛顿方向.

# pseudo code a=0, e= 1e-4 # just a example,人为给定. for k in loops: calculate g_k, H_k, if g_k <= e: break d_k = - H_k^(-1).g_k a_k += d_k

可见,牛顿法没有 步长因子(或 学习率,learning rate), 因此可能出现,越过极值点的情况,甚至可能使结果发散,从而得到非最优结果. 针对此问题, 可通过所谓的线性搜索(line search) 给出步长因子( 设为 \$\\lambda\_k\$):  
\\\[ \\lambda\_k = \\arg \\max\_\{\\lambda \\in \\mathbb\{R\}\} f(a\_k + \\lambda d\_k) \\\]  
于是修上述代码:

# pseudo code a=0, e= 1e-4 # just a example,人为给定. for k in loops: calculate g_k, H_k, if g_k <= e: break d_k = - H_k^(-1) . g_k lambda_ = arg_max(f(a_k,lambda d_k)) a_k += lambda_ . d_k

Ps: 牛顿法中的 \$d\_k\$ 的求解, 一般不会通过求解Hessian矩阵的逆,而是通过方程组:  
\\\[ H\_k \\cdot d\_k = g\_k \\\]  
求解.

## 拟牛顿法 ##

牛顿法优点是收敛速度快.但计算量大( 二阶偏导),且限制较多(Hessian 矩阵需正定).

为解决此问题即出现了拟牛顿法,其基本思想: **不用二阶偏导,而是应用构造法,人为地构造出 近似Hessian 矩阵(或其逆阵)的正定对称阵.** 依据不同的构造方法产生不同的拟牛顿法.

### 拟牛顿条件 ###

构造近似Hessian矩阵(或其逆)的理论依据称为拟牛顿条件.

设B,D 分别表示对Hessian矩阵及其逆阵的近似:\$B \\approx H, D \\approx H^\{-1\}\$.

考虑第 k+1 步:  
\\\[ f(a) \\approx f(a\_\{k+1\}) +\\nabla\_\{a\_\{k+1\}\} f\\cdot(a - a\_\{k+1\} )+ \\frac\{1\}\{2\} (a - a\_\{k+1\})^T\\cdot \\nabla^2 f\_\{a+1\}\\cdot (a- a\_\{k+1\}) \\\]  
两边同时对 a 求偏导:  
\\\[ \\nabla\_a f \\approx \\nabla\_\{a\_\{k+1\}\} f+ (a - a\_\{k+1\})^T\\cdot \\nabla^2 f\_\{a+1\} \\\]  
令 \$a = a\_k\$ 并整理:  
\\\[ g\_\{k+1\} - g\_k \\approx H\_\{k+1\} \\cdot(a\_\{k+1\} - a\_k) \\\]  
再令 \$y\_k = g\_\{k+1\} - g\_\{k\}, s\_k = a\_\{k+1\} - a\_k\$:  
\\\[ y\_k \\approx H\_\{k+1\} \\cdot s\_k\\\\ s\_k \\approx H^\{-1\}\_\{k+1\} \\cdot y\_k \\\]  
上式即为拟牛顿条件.也即构造出的B,D 应满足如下条件:  
\\\[ y\_k = B\_\{k+1\} \\cdot s\_k\\\\ s\_k = D\_\{k+1\} \\cdot y\_k \\\]

### BFGS算法 ###

有了拟牛顿条件就可以构造近似矩阵了. 首先需明确,近似矩阵的构造采用迭代的式. 对Hessian 矩阵的逆阵做近似的算法叫做 DFP 算法,对Hessian矩阵直接构造的叫做BFGS算法, 都是分别以其开发者们的名字首字母命名的. 其中BFGS 的效果会更好,并且其也有较完善的局部收敛理论支持,因此本节介绍BFGS.

BFGS 要做的是 求得近似矩阵 \$B\_k\$ 使得: \$B\_k \\approx H\_k\$. 设其迭代式:  
\\\[ B\_\{k+1\} = B\_k +\\Delta B\_k,\\qquad k= 1,2,\\cdots \\\]  
其中 \$D\_0\$ 设为单位阵 \$I\$. 这里的关键是\$\\Delta B\_k\$的构造, 因为Hessian 矩阵是对称阵,很自然地,想到构造出的\$\\Delta B\_k\$ 也应该是对称阵, 由拟牛顿条件,其中有两个参数(\$s\_k,y\_k\$), 推想\$\\Delta B\_k\$ 应与其二者有关, 也即需两种参数而得\$\\Delta B\_k\$, 因此不防构造成两个对称阵扔组合:  
\\\[ \\Delta B\_k = \\alpha uu^T + \\beta vv^T \\\]  
其中\$\\alpha,\\beta\$ 为系数, 而u,v 为N维向量.代入(13)式:  
\\\[ \\begin\{array\}\\\\ y\_k &=& (B\_k+\\Delta B\_k)\\cdot s\_k\\\\ &=& (B\_k + \\alpha uu^T + \\beta vv^T)\\cdot s\_k\\\\ &=& B\_k \\cdot s\_k + u(\\alpha u^T \\cdot s\_k)+ v(\\beta v^T\\cdot s\_k)\\\\ &=& B\_k \\cdot s\_k + (\\alpha u^T \\cdot s\_k) u+ (\\beta v^T\\cdot s\_k) v\\\\ \\end\{array\} \\\]  
上式中$(\\alpha u^T \\cdot s\_k) ,\\quad (\\beta v^T\\cdot s\_k) $ 均为实数,不妨令其分别为 1, -1则:  
\\\[ \\alpha = \\frac\{1\}\{u^T s\_k\},\\qquad \\beta = -\\frac\{1\}\{v^T s\_k\} \\\]  
(16)式也变成:  
\\\[ y\_k - B\_k s\_k = u - v \\\]  
可令:  
\\\[ u = y\_k,\\qquad v = B\_k s\_k \\\]  
于是:  
\\\[ \\alpha = \\frac\{1\}\{y\_k^T s\_k\},\\qquad \\beta = -\\frac\{1\}\{s^T\_k B\_k s\_k\} \\\]  
最终:  
\\\[ \\Delta B\_k = \\frac\{y\_k y\_k^T\}\{y\_k^T s\_k\} - \\frac\{B\_k s\_k s\_k^T B\_k\}\{s\_k^T B\_k s\_k\} \\\]

#pseudo code a = 0 # 初始化 e = 1e-4 #阈值, 人为给定 B = I # 初始化 k = 0 for k in loops: calculate g_k if g_k < e: break d_k = - B_k^(-1) g_k lambda_ = arg_max(f(a_k,lambda d_k)) s_k = lambda_ d_k a_k += s_k y_k = g_k -g_(k-1) g_(k-1) = g_k B_k += (y_k * y_k.T)/(y_k.T * s_k) - (B_k s_k s_k.T B_k)/(s_k.T * B_k *s_k)

上述伪码中,\$B\_k^\{-1\}\$ 的求解,通过对最后的迭代式应用 Sherman-Morrison公式,直接给出:  
\\\[ \\begin\{array\}\\\\ B\_\{k+1\}^\{-1\} & =& (I - \\frac\{s\_k y\_k^T\}\{y\_k^T s\_k\} )B\_k^\{-1\}(I - \\frac\{y\_k s\_k^T\}\{y\_k^T s\_k\} )+\\frac\{s\_k s\_k^T\}\{y\_k^T s\_k\}\\\\ & =& B\_k^\{-1\} (\\frac\{1\}\{s\_k^T y\_k\}+\\frac\{y\_k^TB\_k^\{-1\}y\_k\}\{(s\_k^T y\_k)^2\} )s\_ks\_k^T - \\frac\{1\}\{s\_k^T y\_k\}(B\_k^\{-1\}y\_ks\_k^T + s\_ky\_k^TB\_k^\{-1\})\\\\ \\end\{array\} \\\]  
上式中最的等式的部分中的括号是两个数.

Ps: Sherman-Morrison 公式, 设\$A\\in \\mathbb\{R\}^N\$为非奇异方阵, u,v 也为N 阶方阵, 若\$1 + v^TA^\{-1\}u \\ne 0\$,则:  
\\\[ (A + uv^T)^\{-1\} = A^\{-1\} - \\frac\{A^\{-1\}uv^T A^\{-1\}\}\{1 + v^TA^\{-1\}u\} \\\]  
将 \$B\_k^\{-1\}\$ 换用\$D\_k\$表示,则:

#pseudo code a = 0 # 初始化 e = 1e-4 #阈值, 人为给定 B = I # 初始化 k = 0 for k in loops: calculate g_k if g_k < e: break d_k = - D_k * g_k lambda_ = arg_max(f(a_k,lambda d_k)) s_k = lambda_ d_k a_k += s_k y_k = g_k -g_(k-1) g_(k-1) = g_k D_k = (I - (s_k*y_k.T)/(y_k.T* s_k))*D_k * (I - (y_k - s_k.T)/y_k.T * s_k) + (s_k*s_k.T)/(y_k.T* s_k)

## L-BFGS算法 ##

当待估参数很多,或待估向量维数很高,需很大存储空间. 为减少内丰开销, L-BFGS(Limited memory BFGS) 应运而生. 其是BFGS的近似,其基本思想: 不再存储完整的矩阵D\_k, 而是存储计算过程中的向量序列\$\\\{s\_k\\\},\\quad \\\{y\_k\\\}\$; 而且只存储最新的 m 个(人为给定). 每次计算D\_k时,只利用最新的 m 个\$\\\{s\_k\\\},\\quad \\\{y\_k\\\}\$, 经此处理,原来的O(N2)变成了O(mN).

由上一个伪码中最后一步:  
\\\[ D\_\{k+1\} = (I - \\frac\{s\_k y\_k^T\}\{y\_k^Ts\_k\})D\_k(I - \\frac\{y\_ks\_k^T\}\{y\_k^Ts\_k\}) + \\frac\{s\_ks\_k^T\}\{y\_k^Ts\_k\} \\\]  
记\$\\rho\_k = \\frac\{1\}\{y\_k^Ts\_k\}\\quad V\_k = I - \\rho\_k y\_k s\_k^T\$ 则:  
\\\[ D\_\{k+1\} = V\_k^TD\_kV\_k + \\rho\_ks\_ks\_k^T \\\]  
上面给出逻辑,下面给出计算方式:

# pseudo code to get D_k * g_k if k <= m: delta = 0 L = k else: delta = k - m L = m q_L = g_k a_list = [] for i in range(L-1,-1,-1): # loop: L-1 to 0 j = i + delta alpha = rho_j* s_j.T* q a_list.append(alpha) q -= alpha* y_i a_list = a_list[::-1] r = D_0 * q for i in range(L): j = i+ delta beta = rho_j*y_j.T * r r +=r +(alpha_i - beta)* s_j # alpha_i = a_list[i]

## 参考文献 ##

牛顿法与拟牛顿法学习笔记一 至 五, peghoty, 2014.

Ps: 本文主要参考上述文献,因主要目的是对此法做一些记录,且peghoty写的博文很清晰,故没有广泛查找文献, 虽有失严谨,但已足够.

转载于:https://www.cnblogs.com/vpegasus/p/newton.html