牛顿法与拟牛顿法
软件设计 2017-08-20
牛顿法与拟牛顿法
优化问题是机器学习中非常重要的部分,无论应用何种算法,构建何种模型,最终我们的目的都是找到最优解的. 那优化算法是无法回避的. 当今机器学习,特别是深度学习中, 梯度下降算法(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写的博文很清晰,故没有广泛查找文献, 虽有失严谨,但已足够.
