1 牛顿法
参考：http://blog.csdn.net/itplus/article/details/21896453
1.1 原始牛顿法
考虑数据是一维的优化问题：
x∗=minxf(x)(1)
若当前 x 已迭代到 xk ，得到的值是 f(xk) ，在 xk 处做二阶泰勒展开：
φ(x)=f(xk)+f′(xk)(x−xk)+12f′′(xk)(x−xk)2(2)
求的是 xk 附近的最值，所以令 φ′(x)=0 得：
f′(xk)+f′′(xk)(x−xk)=0(3)
得到 x 的更新公式：
x=xk−f′(xk)f′′(xk)(4)
若数据是多维：
φ(x)=f(xk)+∇f(xk)⋅(x−xk)+12(x−xk)T∇2f(xk)(x−xk)(5)
记： 梯度g=∇f=⎡⎣⎢⎢⎢⎢⎢⎢⎢∂f∂x1∂f∂x2⋮∂f∂xN⎤⎦⎥⎥⎥⎥⎥⎥⎥海森矩阵H=∇2f=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢∂2f∂x21∂2f∂x2∂x1⋮∂2f∂xN∂x1∂2f∂x1∂x2∂2f∂x22⋮∂2f∂xN∂x2⋯⋯⋱⋯∂2f∂x1∂xN∂2f∂x2∂xN⋯∂2f∂x2N⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
得到：
gk+Hk(x−xk)=0(6)
x=xk−H−1k⋅gk(7)
牛顿法的优点是，对于二次性较强的函数，求解速度较快。缺点是，没有步长因子，不能保证函数值稳定地下降。
1.2 阻尼牛顿法
由上面的分析可以看出，每次更新的方向为牛顿方向 −H−1k⋅gk ，如果给它加上步长，就是阻尼牛顿法：
x=xk−λkH−1k⋅gk(8)
其中， λk 是通过在牛顿方向上做搜索得出的：
λk=argminλf(xk−λkH−1k⋅gk)(9)
牛顿法是梯度法的发展，利用了二阶偏导的信息，即梯度下降的趋势。但是其要求函数具有连续一、二阶导数，海森矩阵正定。并且计算量很大。
2 拟牛顿法
上面牛顿法的计算量大就大在海森矩阵的计算上，不要提计算，就算海森矩阵是不是正定的都无法保证。所以人们发明了一些方法，搞出来一个矩阵来近似海森矩阵，这就是拟牛顿法。
2.1 拟牛顿条件
二阶泰勒展开式可以写成：
f(x)≈f(xk+1)+∇f(xk+1)⋅(x−xk+1)+12(x−xk+1)T∇2f(xk+1)(x−xk+1)(10)
两边作用梯度算子：
gx≈gk+1+Hk+1(x−xk+1)(11)
取 x=xk ，并记 sk=xk+1−xkyk=gk+1−gk：
yk≈Hk+1sk(12)
sk≈H−1k+1yk(13)
所以，我们找到的近似海森矩阵 Bk+1 或者近似逆矩阵 Dk+1，只要满足 yk=Bk+1sk 或sk=Dk+1yk即可。这就是拟牛顿条件。
2.2 近似算法框架
1、给定初始值 x0 和阈值 ϵ ，令 D0=I ；
2、确定搜索方向 dk=−Dk⋅gk 或 dk=−B−1k⋅gk；
3、确定步长 λk ， 令 sk=λkdk, xk+1=xk+sk ；
4、若 ||gk||<ϵ,算法结束；
5、计算 Dk+1 或 Bk+1，k=k+1 转至第2步。
由此可见，算法的关键，就在于怎么构造出Dk+1 或 Bk+1。
2.3 DFP算法
DFP算法通过迭代计算近似海森逆矩阵：
Dk+1=Dk+ΔDk(14)
问题就转化为怎么构造 ΔDk ，我们假设 ΔDk 有如下的形式：
ΔDk=αuuT+βvvT(15)
由(13)式可知：
sk=Dk+1yk=Dkyk+αuuTyk+βvvTyk=Dkyk+(αuTyk)u+(βvTyk)v(16)
小括号里面的是实数，假设：
αuTyk=1βvTyk=−1(17)
即：
α=1uTykβ=−1vTyk(18)
将(17)代入(16)：
u−v=sk−Dkyk(19)
直接令：
u=sk ,v=Dkyk(20)
则(18)变成：
α=1sTkykβ=−1(Dkyk)Tyk=−1yTkDkyk(21)
Dk 是对称矩阵。
将(20)(21)代入(15)得到：
ΔDk=sksTksTkyk−DkykyTkDkyTkDkyk(22)
再代入(14)得：
Dk+1=Dk+sksTksTkyk−DkykyTkDkyTkDkyk(23)
这样就得到了 Dk+1 的更新公式，代入2.2的算法框架就可以了。
2.4 BFGS算法
上面的推导中用的是 Dk ，如果我们换成是 Bk ,经过一顿推导，我们会得到：
Bk+1=Bk+ykyTkyTksk−BksksTkBksTkBksk(24)
但是每次确定搜索方向的时候都要计算 B−1k，更好的做法是用Sherman-Morrison公式：
设 A∈Rn 为非奇异方阵， u,v∈Rn， 若1+vTA−1u≠0， 则有：
(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u
得到：
B−1k+1=(I−skyTkyTksk)B−1k(I−yksTkyTksk)+sksTkyTksk(25)
换成 Dk+1 ：
Dk+1=(I−skyTkyTksk)Dk(I−yksTkyTksk)+sksTkyTksk(26)
代入2.2的算法框架就可以了。BFGS和DFP的不同就在于迭代公式的不同。
2.3 L-BFGS算法
当我们的数据维度很大的时候，存储 Dk 需要很大的空间。为解决这个问题，L-BFGS诞生了。
记 ρk=1yTksk，Vk=I−ρkyksTk，则(26)可以写成：
Dk+1=VTkDkVk+ρksksTk(27)
经过一顿递推推导我们可以得到通项公式：
Dk+1=(VTkVTk−1⋯VT1VT0)D0(V0V1⋯Vk−1Vk)+(VTkVTk−1⋯VT2VT1)ρ0s0sT0(V1V2⋯Vk−1Vk)+(VTkVTk−1⋯VT3VT2)ρ1s1sT1(V2V3⋯Vk−1Vk)+⋯+(VTkVTk−1)ρk−2sk−2sTk−2(Vk−1Vk)+(VTk)ρk−1sk−1sTk−1(Vk)+ρksksTk(28)
上式中，如果我们指定只存储 m 组 sk, yk， 后面计算的时候有新值的时候，依次删除开始的值，保证只有 m 组 s, y。取 m^=min(k,m−1) 则得到近似的计算公式：
Dk+1=(VTkVTk−1⋯VTk−m^+2VTk−m^+1)D0(Vk−m^+1Vk−m^+2⋯Vk−1Vk)+(VTkVTk−1⋯VTk−m^+3VTk−m^+2)ρ0s0sT0(Vk−m^+2Vk−m^+3⋯Vk−1Vk)+⋯+(VTkVTk−1)ρk−2sk−2sTk−2(Vk−1Vk)+(VTk)ρk−1sk−1sTk−1(Vk)+ρksksTk(29)
写到这我们知道，只要存储 m 组 s, y 进行近似计算就好了。之后又有人搞出了快速计算搜索方向 dk 的方法（毕竟计算Dk就是要计算搜索方向）：
1、初始化
δ={0k−mif k≤mif k>m;L={kmif k≤mif k>m;qL=gk
2、后向循环
FOR i=L−1,L−2,⋯,1,0 DO
{
j=i+δ
αi=ρjsTjqi+1
qi=qi+1−αiyj
}
3、前向循环
r0=D0⋅q0
FOR
i=0,1,⋯,L−2,L−1 DO
{
j=i+δ
βj=ρjyTjri
ri+1=ri−(αi−βi)sj
}
搜索方向： dk=−rk。
