牛顿法与拟牛顿法
优化问题是机器学习中非常重要的部分,无论应用何种算法,构建何种模型,最终我们的目的都是找到最优解的. 那优化算法是无法回避的. 当今机器学习,特别是深度学习中, 梯度下降算法(gradient descent algorithm) 可谓炙手可热. 不过优化算法不只其一种,其他算法也很常用,也很优秀,比如今天介绍的牛顿法(Newton methods)与拟牛顿法(Quasi Newton methods).
关于此算法的文献很多, 这里算是一个分享,也算是对自己知识的复习. 本文主要参考 peghoty 的博文.
考虑如下的约束问题:
其中 \(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)\) 是原始待优化函数; \(a^*\) 是当前极小值点(估计值). 由极值条件:
即:
可见这里要求 二阶偏导矩阵(hessian 矩阵可逆(非奇异), 幸好在实际中这个经常满足).
为简化, 梯度向量(gradient vector) 与海森矩阵(hessian matrix) 分别用 \(g_k =\nabla_{a^*} f^T , H_k = ( (\nabla_{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,\beta\) 为系数, 而u,v 为N维向量.代入(13)式:
上式中$(\alpha u^T \cdot s_k) ,\quad (\beta v^T\cdot 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 公式, 设\(A\in \mathbb{R}^N\)为非奇异方阵, u,v 也为N 阶方阵, 若\(1 + v^TA^{-1}u \ne 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).
由上一个伪码中最后一步:
记\(\rho_k = \frac{1}{y_k^Ts_k}\quad V_k = I - \rho_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写的博文很清晰,故没有广泛查找文献, 虽有失严谨,但已足够.
