无约束最优化方法之牛顿法、拟牛顿法、BFGS、LBFGS及若干

好久没写博客了，今天打开一看****终于可以用latex，不用到处去粘贴标签，方便了许多。且先试试效果如何。先讲讲一些优化方法。

最速下降法
牛顿法
拟牛顿法
SR1
BFGS
DFP
LBFGS

【最速下降法】

无约束最优化方法不涉及约束条件，所以都是介绍如何寻找搜索方向以及搜索步长。
无约束最优化问题的目标函数：

min x \in R n f (x)

$\min_{x\in R^n}\quad f(x)$
感觉这latex还是有些别扭，稍不留意就直接当做字符处理了。
还是首先介绍一下梯度下降，梯度下降学过优化的都很清楚，一般叫最速下降法，这个方法有两点，首先是

x $x$ 更新的方向是负梯度方向，第二个是沿着该方向搜索，找到该方向的最小值所对应的

x $x$ 就是下次更新的值。梯度下降是最简单的一种方法，但是很多情况下却并不使用这种方法，原因是收敛速率比较慢，问题出在第二步上，由于搜索搜索时一直打到该方向的最小值，那么很显然，继续沿着该方向搜索会使函数值变小，函数梯度与搜索方向夹角大于九十度，所以该点的梯度和搜索方向在此时正交，这样相邻搜索点的梯度就会呈现锯齿状，函数沿着锯齿状下降，严重降低目标函数的收敛速率。
梯度下降的递推公式推导是根据函数的一阶泰勒展开近似得到的。将

f(x) $f(x)$ 在

x (k) $x^{(k)}$ 附近进行一阶泰勒展开：

f (x) \approx f (x (k)) + g T k (x - x (k))

$f(x)\approx f(x^{(k)})+g^T_k(x-x^{(k)})$
这里,

g k =g(x (k) )=∇f(x (k) ) $g_k=g(x^{(k)})=\nabla f(x^{(k)})$ 为

f(x) $f(x)$ 在

x (k) $x^{(k)}$ 的梯度。
那么第

k+1 $k+1$ 次的迭代值就可以通过：

x (k + 1) \leftarrow x (k) + λ k p k

$x^{(k+1)}\leftarrow x^{(k)}+\lambda_{k} p_k$ .
其中

p k $p_k$ 是搜索方向，取负梯度方向

p k =−∇f(x (k) ) $p_k=-\nabla f(x^{(k)})$ 可以使函数下降最快，

λ k $\lambda _k$ 是步长，并且取

λ k $\lambda_k$ 使得

f (x (k) + λ k p k) = min λ \geq 0 f (x (k) + λ p k)

$f(x^{(k)}+\lambda_kp_k)=\min_{\lambda\ge0} f(x^{(k)}+\lambda p_k)$
最速下降法就是这样，不断地寻找搜索方向以及确定搜索步长，直到达到终止条件，相邻函数值相遇某个阈值或是

x (k) $x^{(k)}$ 和

x (k+1) $x^{(k+1)}$ 小于某个阈值。但是产生的问题就是最速下降在接近终点的时候收敛速度较慢，容易之字形收敛。当然步长也不必是取该方向下降尽头的值，可以取固定值，但是太大容易发散，太小收敛速率比较慢。
关于随机梯度下降法与批量下降法，大多数用梯度下降是求无约束目标函数，例如求经验损失最小时函数的参数，含有大量的训练数据。批量下降法是同时使用所有数据对梯度进行更新，很显然需要好多次迭代。随机梯度下降是每次只使用一个数据对函数参数进行更新，这样往往只通过一部分数据更新参数就会收敛，但是由于每次根据一个数据跟新，容易造成噪音问题。

【牛顿法】

由于最速梯度下降收敛速度并不“最速”，局部搜索中最速下降的方向和全局的最小值方向并不一致，所以就有后来改进的方法，包括牛顿法以及拟牛顿法。
牛顿法要求 $f(x)$ 具有二阶连续可导性。
仍然考虑无约束最优化问题的目标函数：

min x \in R f (x)

$\min_{x\in R}\quad f(x)$
这里所不同的是进行二阶泰勒展开：

f (x) \approx f (x (k)) + g T k (x - x (k)) + 1 2 (x - x (k)) T H (x (k)) (x - x (k))

$f(x)\approx f(x^{(k)})+g^T_k(x-x^{(k)})+\dfrac{1}{2}(x-x^{(k)})^TH(x^{(k)})(x-x^{(k)})$
这里,

g k =g(x (k) )=∇f(x (k) ) $g_k=g(x^{(k)})=\nabla f(x^{(k)})$ 为

f(x) $f(x)$ 在

x (k) $x^{(k)}$ 的梯度。

H(x (k) ) $H(x^{(k)})$ 是

f(x) $f(x)$ 的海塞矩阵

H (x) = [\partial 2 f ( x ) \partial x i \partial x j] n \times n

$H(x)=\begin{matrix}[ {\dfrac{\partial ^2 f(x)}{\partial x_i \partial x_j}}]_{n\times n}\end{matrix}$
显然，

f(x) $f(x)$ 有极值的条件是在

x k $x^k$ 处的一阶导数为0，

∇f(x)=0 $\nabla f(x)=0$ ,所以，当我们从

x k $x^k$ 处开始搜索时，搜索终止处

x k+1 $x^{k+1}$ 应该满足

∇f(x (k+1) )=0 $\nabla f(x^{(k+1)})=0$ 。所以我们对二阶近似求导。

\nabla f (x) = g k + H k (x - x (k))

$\nabla f(x)=g_k+H_k(x-x^{(k)})$
所以

g k + H k (x - x (k)) = 0

$g_k+H_k(x-x^{(k)})=0$
then,

x (k + 1) = x (k) - H - 1 k g k

$x^{(k+1)}=x^{(k)}-H^{-1}_kg_k$
经典牛顿法虽然具有二次收敛性，但是要求初始点需要尽量靠近极小点，否则有可能不收敛。计算过程中需要计算目标函数的二阶偏导数，难度较大。更为复杂的是目标函数的Hesse矩阵无法保持正定，会导致算法产生的方向不能保证是f在Xk 处的下降方向，从而令牛顿法失效；特别的，如果Hesse矩阵奇异，牛顿方向可能根本是不存在的。

拟牛顿法

上面说了，虽然牛顿法能够具有二次收敛性，但是要求太高，个别情况下甚至无法求出牛顿法的迭代方向，所以就有了拟牛顿法，来对Hesse矩阵的逆进行近似。
通过泰勒二阶近似可以得到：

\nabla f (x k + 1) = \nabla f (x k) + H k (x (k + 1) - x k)

$\nabla{f(x^{k+1})}=\nabla{f(x^{k})}+H_k(x^{(k+1)}-x^k)$
令，

y k = \nabla f (x k + 1) - \nabla f (x k), s k = x (k + 1) - x k

$y_k=\nabla{f(x^{k+1})}-\nabla{f(x^{k})},s_k=x^{(k+1)}-x^k$
then,

y k = H k s k

$y_k=H_ks_k$
或者说，

H - 1 k y k = s k

$H^{-1}_ky_k=s_k$
注意到，

s k = x (k + 1) - x (k) = α d k

$s_k=x^{(k+1)}-x^{(k)}=\alpha d_k$ ,所以拟牛顿法模拟了牛顿的方向。
所以，拟牛顿法选取满足条件

B k s k =y k $B_ks_k=y_k$ ,

B k $B_k$ 作为Hesse矩阵

H k $H_k$ 的近似，或者

s k =G k y k $s_k=G_ky_k$

G k $G_k$ 作为hesse矩阵逆的近似，而且要使得计算简便。当有了

B k $B_k$ 之后，通过对

B k $B_k$ 进行低秩修改得到

B k+1 $B_{k+1}$ ,

B k + 1 = B k + Δ k

$B_{k+1}=B_{k}+\Delta_k$
使其仍满足近似条件。
一般，最初始

B k $B_k$ 都是使用单位矩阵或者随机初始化。

SR1

根据修改 $B_k$ 方法的不同，衍生出很多不同的方法，最简单的就是给 $B_{k-1}$ 加上一个秩为1的对称矩阵，由于秩为1的对称矩阵可以写成一个列向量和其转置相乘的形式，所以 $B_k$ 的约束条件可以写成：

(B k - 1 + β k u k u T k) s k = y k

$(B_{k-1}+\beta_ku_ku^T_k)s_k=y_k$
展开得到：

B k - 1 s k + β k u k u T k s k = y k

$B_{k-1}s_k+\beta_ku_ku^T_ks_k=y_k$
注意到

u T k s k $u^T_ks_k$ 是个常数，所以，

- B k - 1 s k + y k = (β k u T k s k) u k

$-B_{k-1}s_k+y_k=(\beta_ku^T_ks_k)u_k$
所以我们可以选

β k $\beta_k$ 使其满足

β k u T k s k =1 $\beta_ku^T_ks_k=1$

u k = y k - B k - 1 s k, β k = 1 u T k s k = 1 s T k u k = 1 s T k ( y k - B k - 1 s k )

$u_k=y_k-B_{k-1}s_k,\beta_k=\frac{1}{u^T_ks_k}=\frac{1}{s^T_ku_k}=\frac{1}{s^T_k(y_k-B_{k-1}s_k)}$
最后得到

B k+1 $B_{k+1}$ 的更新式子

B k + 1 = B k + ( y k - B k - 1 s k ) ( y k - B k - 1 s k ) T s T k ( y k - B k - 1 s k )

$B_{k+1}=B_k+\frac{(y_k-B_{k-1}s_k)(y_k-B_{k-1}s_k)^T}{s^T_k(y_k-B_{k-1}s_k)}$
当然，通过

G k $G_k$ 也能得到类似的式子，

BFGS

BFGS方法是一种秩2近似，至于为什么使用秩2近似这个暂时还不得而知。先讲一下是如何推导的。
BFGS是近似海瑟矩阵 $H$ ,首先，相应的牛顿条件是

B k + 1 s k = y k

$B_{k+1}s_k=y_k$
使用秩2近似，

B k + 1 = B k + P k + Q k = B k + α u k u T k + β v k v T k

$B_{k+1}=B_k+P_k+Q_k =B_k+\alpha u_ku^T_k+\beta v_kv^T_k$
所以，

B k + 1 s k = (B k + P k + Q k) s k = B k s k + α u k u T k s k + β v k v T k s k = y k

$B_{k+1}s_k=(B_k+P_k+Q_k)s_k =B_ks_k+\alpha u_ku^T_ks_k+\beta v_kv^T_ks_k =y_k$

B k + 1 s k = B k s k + (α u T k s k) u k + (β v T k s k) v k = y k

$B_{k+1}s_k=B_ks_k+(\alpha u^T_ks_k)u_k+(\beta v^T_ks_k)v_k=y_k$
由于满足条件的

α,β,u k ,v k $\alpha,\beta,u_k,v_k$ 相当多，所以可以这样设置，

α u T k s k = 1, β v T k s k = 1

$\alpha u^T_ks_k=1,\beta v^T_ks_k=1$

α = 1 u T k s k, β = 1 v T k s k

$\alpha=\frac{1}{u^T_ks_k} ,\beta=\frac{1}{ v^T_ks_k}$
这样式子就成了

B k + 1 s k = B k s k + u k + v k = y k

$B_{k+1}s_k=B_ks_k+u_k+v_k=y_k$
令

u k =y k ,B k s k +v k =0,v k =−B k s k $u_k=y_k,B_ks_k+v_k=0,v_k=-B_ks_k$
所以(

B k $B_k$ 是对称的)

B k = B k + α u k u T k + β v k v T k

$B_k=B_k+\alpha u_ku^T_k+\beta v_kv^T_k$

= B k + y k y T k y T k s k - B k s k s T k B k s T k B k s k

$=B_k+\frac{y_ky^T_k}{y^T_ks_k}-\frac{B_ks_ks^T_kB_k}{s^T_kB_ks_k}$
我们使用的

B k $B_k$ 的逆，所以这里还需要使用Sherman-Morrison公式，假设A是n阶可逆矩阵，

u,v $u,v$ 是n维向量，且

A+uv T $A+uv^T$ 也是可逆矩阵，则

(A + u v T) - 1 = A - 1 - A - 1 u v T A - 1 1 + v T A - 1 u

$(A+uv^T)^{-1}=A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}$
得到

B - 1 k + 1 = (I - s k y T k y T k s k) B - 1 k (I - y k s T k y T k s k) + s k s T k y T k s k

$B_{k+1}^{-1}=(I-\frac{s_ky_k^T}{y_k^Ts_k})B_k^{-1}(I-\frac{y_ks_k^T}{y_k^Ts_k})+\frac{s_ks_k^T}{y_k^Ts_k}$

DFP

DFP推导方法和BFGS类似，只不过是对hesse矩阵的逆进行近似,略。

LBFGS

关于LBFGS的推导，可以参考【3】和【4】，主要是通过BFGS的最后目标式子，不再保留完整的矩阵B_k^{-1}，因为当维度很大的时候（n>10^4），需要的空间非常大，所以保留了一些计算 $B_k^{-1}$ 需要的 $s_k,y_k$ 序列，而且只保存最近的m个序列。
这里不妨用 $H_k$ 表示 $B_k^{-1}$ ,非hesse矩阵.

H k + 1 = (I - s k y T k y T k s k) H k (I - y k s T k y T k s k) + s k s T k y T k s k

$H_{k+1}=(I-\frac{s_ky_k^T}{y_k^Ts_k})H_k(I-\frac{y_ks_k^T}{y_k^Ts_k})+\frac{s_ks_k^T}{y_k^Ts_k}$
define:

ρ k =1y T k s k $\rho_k=\frac{1}{y_k^Ts_k}$ ,

V k =I−ρ k y k s T k $V_k=I-\rho_ky_ks_k^T$ ,then the above formulation can be rewritten as:

H k + 1 = V T k H k V k + ρ k s k s T k

$H_{k+1}=V_k^TH_kV_k+\rho_ks_ks_k^T$
Then,recursively

H 1 = V T 0 H 0 V 0 + ρ 0 s 0 s T 0

$H_{1}=V_{0}^TH_{0}V_{0}+\rho_{0}s_{0}s_{0}^T$

H 2 = = = V T 1 H 1 V 1 + ρ 1 s 1 s T 1 V T 1 (V T 0 H 0 V 0 + ρ 0 s 0 s T 0) V 1 + ρ 1 s 1 s T 1 V T 1 V T 0 H 0 V 0 V 1 + V T 1 ρ 0 s 0 s T 0) V 1 + ρ 1 s 1 s T 1

$\begin{align} H_{2}=&V_{1}^TH_{1}V_{1}+\rho_{1}s_{1}s_{1}^T \\ =&V_{1}^T(V_{0}^TH_{0}V_{0}+\rho_{0}s_{0}s_{0}^T)V_{1}+\rho_{1}s_{1}s_{1}^T\\ =&V_{1}^TV_{0}^TH_{0}V_{0}V_{1}+V_{1}^T\rho_{0}s_{0}s_{0}^T)V_{1}+\rho_{1}s_{1}s_{1}^T \end{align}$
所以就有了这个公式：

H k + 1 = + + + + (V T k V T k - 1 . . . V T 1 V T 0) H 0 (V 0 V 1 . . . V k - 1 V k) (V T k V T k - 1 . . . V T 1) ρ 1 s 1 s T 1 (V 1 . . . V k - 1 V k) . . . (V T k) ρ k - 1 s k - 1 s T k - 1 (V k) ρ k s k s T k

$\begin{align} H_{k+1}=&(V_{k}^TV_{k-1}^T...V_{1}^TV_{0}^T)H_{0}(V_{0}V_{1}...V_{k-1}V_{k})\\ +&(V_{k}^TV_{k-1}^T...V_{1}^T)\rho_1s_1s_1^T(V_{1}...V_{k-1}V_{k})\\ +&...\\ +&(V_{k}^T)\rho_{k-1}s_{k-1}s_{k-1}^T(V_{k})\\ +&\rho_{k}s_{k}s_{k}^T \end{align}$
然后为了算这个式子，需要不断迭代LBFGS原著中给了一个两层的递推程序求这个式子，只保留最近m步：

H k + 1 = + + . . . + + (V T k V T k - 1 . . . V T k - m) H 0 (V k - m . . . V k - 1 V k) (V T k V T k - 1 . . . V T k - m + 1) ρ k - m s k - m s T k - m (V k - m + 1 . . . V k - 1 V k) (V T k) ρ k - 1 s k - 1 s T k - 1 (V k) ρ k s k s T k

$\begin{align} H_{k+1}=&(V_{k}^TV_{k-1}^T...V_{k-m}^T)H_{0}(V_{k-m}...V_{k-1}V_{k})\\ +&(V_{k}^TV_{k-1}^T...V_{k-m+1}^T)\rho_{k-m}s_{k-m}s_{k-m}^T(V_{k-m+1}...V_{k-1}V_{k})\\ +...\\ +&(V_{k}^T)\rho_{k-1}s_{k-1}s_{k-1}^T(V_{k})\\ +&\rho_{k}s_{k}s_{k}^T \end{align}$
更新的方向：

H k + 1 \nabla f (x) = + + + + (V T k V T k - 1 . . . V T k - m) H 0 (V k - m . . . V k - 1 V k) \nabla f (x) (V T k V T k - 1 . . . V T k - m + 1) ρ k - m s k - m s T k - m (V k - m + 1 . . . V k - 1 V k) \nabla f (x) . . . (V T k) ρ k - 1 s k - 1 s T k - 1 (V k) \nabla f (x) ρ k s k s T k \nabla f (x)

$\begin{align} H_{k+1}\nabla f(x)=&(V_{k}^TV_{k-1}^T...V_{k-m}^T)H_{0}(V_{k-m}...V_{k-1}V_{k})\nabla f(x)\\ +&(V_{k}^TV_{k-1}^T...V_{k-m+1}^T)\rho_{k-m}s_{k-m}s_{k-m}^T(V_{k-m+1}...V_{k-1}V_{k})\nabla f(x)\\ +&...\\ +&(V_{k}^T)\rho_{k-1}s_{k-1}s_{k-1}^T(V_{k})\nabla f(x)\\ +&\rho_{k}s_{k}s_{k}^T\nabla f(x)\\ \end{align}$
所谓的Two-loop算法：

$q_k\leftarrow \nabla f(x^k)$
对
$i=k-1$ to $k-m$
$\alpha_i=\rho_i s_i^Tq_{i+1}$
$q_i=q_{i+1}-\alpha_i y_i$
然后第二次循环，
根据 wiki LBFGS 【5】
$H_0=\frac{y^T_{k-1}s_{k-1}}{y^T_{k-1}y_{k-1}}$
初始化： $r_{k-m-1}=H_0 q_{k-m}$
对于 $i=k-m,k-m+1$ to $k-1$
$\beta_i=\rho_i y_i^Tr_{i-1}$
$r_i=r_{i-1}+s_i(\alpha_i-\beta_i)$
最后得到的 $r$ 即为所求。上面的q以及 r都只有最后一步结果，中间结果的可以用一个变量代替。

参考：
【1】http://blog.****.net/lilyth_lilyth/article/details/8973972
【2】统计学习方法
【3】http://blog.****.net/lansatiankongxxc/article/details/38801863
【4】http://blog.****.net/zhirom/article/details/38332111
【5】http://en.wikipedia.org/wiki/Limited-memory_BFGS

无约束最优化步骤之牛顿法、拟牛顿法、BFGS、LBFGS及若干

【最速下降法】

【牛顿法】

拟牛顿法

SR1

BFGS

DFP

LBFGS

相关推荐