牛顿法与拟牛顿法
优化问题是机器学习中非常重要的部分,无论应用何种算法,构建何种模型,最终我们的目的都是找到最优解的. 那优化算法是无法回避的. 当今机器学习,特别是深度学习中, 梯度下降算法(gradient descent algorithm) 可谓炙手可热. 不过优化算法不只其一种,其他算法也很常用,也很优秀,比如今天介绍的牛顿法(Newton methods)与拟牛顿法(Quasi Newton methods).
关于此算法的文献很多, 这里算是一个分享,也算是对自己知识的复习. 本文主要参考 peghoty 的博文.
考虑如下的约束问题:
其中 (f(a): mathbb R^N ightarrow mathbb R) 为优化函数, 而 a 是待优化参数向量 $ a = (a_1,a_2,dots,a_N)^T in mathbb{R}^N$
(Ps: 这里回避了 用 x 表示待优参数, 因为x 在机器学习中一般指的是数据,而 实际问题中, 待优化的几乎都是模型参数而非数据). 这里要求 f 是凸的, 且是二阶连续可微的(否则后面用到的二阶偏导无法满足).
牛顿法
牛顿法的基本思想是 用泰勒展开式来替代原有函数: 在现有极小值点附近对 (f(x)) 二阶展开(忽略高阶项,因此是一种近似). 这样的做的好处是可以很容易找到优化条件,并且可以用到二阶偏导条件:
其中(psi(a)) 是原始待优化函数; (a^*) 是当前极小值点(估计值). 由极值条件:
即:
可见这里要求 二阶偏导矩阵(hessian 矩阵可逆(非奇异), 幸好在实际中这个经常满足).
为简化, 梯度向量(gradient vector) 与海森矩阵(hessian matrix) 分别用 (g_k = abla_{a^*} f^T , H_k = ( ( abla_{a^*}^2f )^{-1})^T) 表示.则上式可写成:
写成迭代形式:
这里的 (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)):
于是修上述代码:
# 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矩阵的逆,而是通过方程组:
求解.
拟牛顿法
牛顿法优点是收敛速度快.但计算量大( 二阶偏导),且限制较多(Hessian 矩阵需正定).
为解决此问题即出现了拟牛顿法,其基本思想: 不用二阶偏导,而是应用构造法,人为地构造出 近似Hessian 矩阵(或其逆阵)的正定对称阵. 依据不同的构造方法产生不同的拟牛顿法.
拟牛顿条件
构造近似Hessian矩阵(或其逆)的理论依据称为拟牛顿条件.
设B,D 分别表示对Hessian矩阵及其逆阵的近似:(B approx H, D approx H^{-1}).
考虑第 k+1 步:
两边同时对 a 求偏导:
令 (a = a_k) 并整理:
再令 (y_k = g_{k+1} - g_{k}, s_k = a_{k+1} - a_k):
上式即为拟牛顿条件.也即构造出的B,D 应满足如下条件:
BFGS算法
有了拟牛顿条件就可以构造近似矩阵了. 首先需明确,近似矩阵的构造采用迭代的式. 对Hessian 矩阵的逆阵做近似的算法叫做 DFP 算法,对Hessian矩阵直接构造的叫做BFGS算法, 都是分别以其开发者们的名字首字母命名的. 其中BFGS 的效果会更好,并且其也有较完善的局部收敛理论支持,因此本节介绍BFGS.
BFGS 要做的是 求得近似矩阵 (B_k) 使得: (B_k approx H_k). 设其迭代式:
其中 (D_0) 设为单位阵 (I). 这里的关键是(Delta B_k)的构造, 因为Hessian 矩阵是对称阵,很自然地,想到构造出的(Delta B_k) 也应该是对称阵, 由拟牛顿条件,其中有两个参数((s_k,y_k)), 推想(Delta B_k) 应与其二者有关, 也即需两种参数而得(Delta B_k), 因此不防构造成两个对称阵扔组合:
其中(alpha,eta) 为系数, 而u,v 为N维向量.代入(13)式:
上式中$(alpha u^T cdot s_k) ,quad (eta v^Tcdot s_k) $ 均为实数,不妨令其分别为 1, -1则:
(16)式也变成:
可令:
于是:
最终:
#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公式,直接给出:
上式中最的等式的部分中的括号是两个数.
Ps: Sherman-Morrison 公式, 设(Ain mathbb{R}^N)为非奇异方阵, u,v 也为N 阶方阵, 若(1 + v^TA^{-1}u e 0),则:
将 (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).
由上一个伪码中最后一步:
记( ho_k = frac{1}{y_k^Ts_k}quad V_k = I - ho_k y_k s_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写的博文很清晰,故没有广泛查找文献, 虽有失严谨,但已足够.