梯度-牛顿-拟牛顿优化算法和实现
2015-05-30 10:19
316 查看
作者：金良（golden1314521@gmail.com） csdn博客： http://blog.csdn.net/u012176591
要求解的问题
线搜索技术和Armijo准则
最速下降法及其Python实现
牛顿法
阻尼牛顿法及其Python实现
修正牛顿法法及其Python实现
拟牛顿法
DFP算法及其Python实现
BFGS算法及其Python实现
Broyden族算法及其Python实现
L-BFGS算法及其Python实现
参考文献
minx∈R2f(x)=100(x21−x2)2+(x1−1)2\min_{x\in\mathbb{R}^2}f(x)=100(x_1^2-x_2)^2+(x_1-1)^2
该问题有精确解x∗=(1,1)T,f(x∗)=0x^*=(1,1)^T,f(x^*)=0
其梯度向量g(x)=(400x1(x21−x2)+2(x1−1),−200(x21−x2))g(x)=(400x_1(x_1^2-x_2)+2(x_1-1),-200(x_1^2-x_2))
其海森矩阵G(x)=[1200x21−400x2+2−400x1−400x1200]G(x) = \begin{bmatrix}
1200x_1^2-400x_2+2&-400x_1 \\[0.3em]
-400x_1&200
\end{bmatrix}
下面定义了三个Python语言编写的函数：函数表达式fun，梯度向量gfun，和海森矩阵hess。这三个表达式在后面各个算法的实现中会用到。
符号约定:
gkg_k : ∇f(xk)\nabla f(x_k)，即第目标函数关于kk次迭代值xkx_k的导数。
GkG_k : G(xk)=∇2f(xk)G(x_k)=\nabla^2f(x_k)，即海森矩阵。
dkd_k : 第kk次迭代的优化方向。在最速下降算法中,有dk=−gkd_k = -g_k
αk\alpha_k : 第kk次迭代的步长因子，有xk+1=xk+αkdkx_{k+1}=x_k+\alpha_k d_k
在精确线性搜索中，步长因子αk\alpha_k由下面的式子确定：
αk=argminαf(xk+αdk)\alpha_k=\arg\min_{\alpha} f(x_k+\alpha d_k)
而对于非精确线性搜索，选取的αk\alpha_k只要使得目标函数ff得到可接受的下降量，即
Δfk=f(xk)−f(xk+αkdk)>0\Delta f_k = f(x_k)-f(x_k+\alpha_kd_k)>0
Armijo准则用于非精确线性搜索中步长因子α\alpha的确定，内容如下：
Armijo准则：
已知当前位置xkx_k和优化的方向dkd_k ，参数β∈(0,1)\beta\in(0,1),δ∈(0,0.5)\delta\in(0,0.5).令步长因子αk=βmk\alpha_k=\beta^{m_k},其中mkm_k为满足下列不等式的最小非负整数mm：
f(xk+βmdk)≤f(xk)+δβmgTkdkf(x_k+\beta^m d_k)\leq f(x_k)+\delta\beta^mg_k^Td_k
由此确定的下一个位置xk+1=xk+αkdkx_{k+1}=x_k+\alpha_k d_k
该准则在接下来的介绍的几个算法中多次使用。
step 1 :选取初始点x0∈Rnx_0\in\mathbb{R}^n,容许误差0≤ϵ≪10\leq\epsilon\ll 1 .令 k←1k\leftarrow1.
step 2 :计算gk=∇f(xk)g_k=\nabla f(x_k). 若||gk||≤ϵ||g_k||\leq\epsilon，停止迭代，输出xkx_k作为近似最优解。
step 3 :取方向dk=−gkd_k=-g_k.
step 4 :由线搜索技术确定步长因子αk\alpha_k.
step 5 :令xk+1←xk+αkdkx_{k+1}\leftarrow x_k+\alpha_kd_k，转step 2.
上述step 3中步长因子αk\alpha_k的确定既可以使用精确线搜索方法，也可以使用非精确线搜索方法，在理论上都能保证其全局收敛性。若采用精确线搜索方法，则αk\alpha_k满足
ddαf(xk+αdk)|α=αk=∇f(xk+αkdk)Tdk=0\frac{\mathrm{d}}{\mathrm{d}\alpha}f(x_k+\alpha d_k)|_{\alpha=\alpha_k}=\nabla f(x_k+\alpha_kd_k)^Td_k=0
从而有∇f(xk+1)T∇f(xk)=0\nabla f(x_{k+1})^T\nabla f(x_k)=0
即xk+1x_{k+1}处的梯度与其前驱xkx_k处的梯度是正交的，也就是说，迭代点序列所走的路线是锯齿形的，这造成其收敛速度很缓慢。如下图所示，其中绿色折线所显示的路线就是由最速下降法得到的，红色曲线是由共轭梯度下降法确定的，通过对比可以看出梯度下降法所走的路线是锯齿形的，经测试，其收敛速度非常慢。
最速下降法的Python实现。
性能测试
可以看到大约有一半的样本的迭代次数要超过2000次。
牛顿法推导
函数f(x)f(x)在xkx_k处的泰勒展开式的前三项得
T(f,xk,3)=fk+gTk(x−xk)+12(x−xk)TGk(x−xk)T(f,x_k,3)=f_k+g_k^T(x-x_k)+\frac{1}{2}(x-x_k)^TG_k(x-x_k)
则其稳定点∇T=gk+Gk(x−xk)=0\nabla T = g_k+G_k(x-x_k)=0
若GkG_k非奇异，那么解上面的线性方程（标记其解为xk+1x_{k+1}）即得到牛顿法的迭代公式xk+1=xk−G−1kgkx_{k+1}=x_k-G_k^{-1}g_k
即优化方向为dk=−G−1kgkd_k=-G_k^{-1}g_k
求稳定点是用到了以下公式：
考虑y=xTAxy=x^TAx，根据矩阵求导法则，有
dydx=Ax+ATxd2ydx2=A+AT\frac{\mathrm{d}y}{\mathrm{d}x}=Ax+A^Tx\\\frac{\mathrm{d}^2y}{\mathrm{d}x^2}=A+A^T
注意实际中dkd_k的是通过求解线性方程组Gkd=−gkG_kd=-g_k获得的。
阻尼牛顿法的算法描述：
step 1: 给定终止误差0≤ϵ≪1,δ∈(0,1),σ∈(0,0.5)0\leq\epsilon\ll 1,\delta\in(0,1),\sigma\in(0,0.5). 初始点x0∈Rnx_0\in\mathbb{R}^n. 令k←0k\leftarrow0
step 2: 计算gk=∇f(xk)g_k=\nabla f(x_k). 若||gk||≤ϵ||g_k||\leq\epsilon，停止迭代，输出x∗≈xkx^*\approx x_k
step 3: 计算Gk=∇2f(xk)G_k=\nabla^2f(x_k)，并求解线性方程组得到解dkd_k,Gkd=−gkG_kd=-g_k
step 4: 记mkm_k是满足下列不等式的最小非负整数mm.f(xk+βmdk)≤f(xk)+δβmgTkdkf(x_k+\beta^m d_k)\leq f(x_k)+\delta\beta^mg_k^Td_k
step 5: 令αk=δmk,xk+1=xk+αkdk,k←k+1\alpha_k=\delta^{m_k},x_{k+1}=x_k+\alpha_kd_k,k\leftarrow k+1，转step 2
阻尼牛顿法的Python实现：
性能展示
作图代码：
下面介绍两种修正方法，分别是“牛顿-最速下降混合算法”和“修正牛顿法”。
牛顿-最速下降混合算法
该方法将牛顿法和最速下降法结合起来，基本思想是：当海森矩阵正定时，采用牛顿法确定的优化方向作为搜索方向；否则，即海森矩阵为非正定矩阵时，则采用负梯度方向作为搜索方向。
修正牛顿法
引入阻尼因子μk≥0\mu_k\geq0，即在每一迭代步适当选取参数μk\mu_k，使得矩阵Ak=Gk+μkIA_k=G_k+\mu_k I正定，用AkA_k代替GkG_k来求解dkd_k。
算法描述：
step 1: 给定终止误差0≤ϵ≪1,δ∈(0,1),σ∈(0,0.5)0\leq\epsilon\ll 1,\delta\in(0,1),\sigma\in(0,0.5). 初始点x0∈Rnx_0\in\mathbb{R}^n. 令k←0k\leftarrow0
step 2: 计算gk=∇f(xk)g_k=\nabla f(x_k)，μk=||gk||1+τ\mu_k = ||g_k||^{1+\tau}. 若||gk||≤ϵ||g_k||\leq \epsilon，停止迭代，输出x∗≈xkx^*\approx x_k
step 3: 计算Gk=∇2f(xk)G_k=\nabla^2f(x_k)，并求解线性方程组得到解dkd_k,(Gk+μkI)d=−gk(G_k+\mu_k I)d=-g_k
step 4: 记mkm_k是满足下列不等式的最小非负整数mm.f(xk+βmdk)≤f(xk)+δβmgTkdkf(x_k+\beta^m d_k)\leq f(x_k)+\delta\beta^mg_k^Td_k
step 5: 令αk=δmk,xk+1=xk+αkdk,k←k+1\alpha_k=\delta^{m_k},x_{k+1}=x_k+\alpha_kd_k,k\leftarrow k+1，转step 2
修正牛顿法的Python实现代码：
性能展示
但当海森矩阵G(xk)=∇2f(x)G(x_k)=\nabla^2f(x) 不正定时，不能保证所产生的方向是目标函数在xkx_k处的下降方向。
特别地，当G(xk)G(x_k)奇异时，算法就无法继续进行下去。尽管修正牛顿法可以克服这一缺陷，但修正参数的取值很难把握，过大或过小都会影响到收敛速度。
牛顿法的每一步迭代都需要目标函数的海森矩阵G(xk)G(x_k)，对于大规模问题其计算量是惊人的。
拟牛顿法的基本思想是用海森矩阵GkG_k的某个近似矩阵BkB_k取代GkG_k. BkB_k通常具有下面三个特点：
在某种意义下有Bk≈GkB_k\approx G_k ,使得相应的算法产生的方向近似于牛顿方向，确保算法具有较快的收敛速度。
对所有的kk，BkB_k是正定的，从而使得算法所产生的方向是函数ff在xkx_k处下降方向。
矩阵BkB_k更新规则比较简单
设 f:Rn→Rf:\mathbb{R}^n\rightarrow\mathbb{R}在开集D⊂RnD\subset \mathbb{R}^n上二次连续可微，那么ff在xk+1x_{k+1}处的二次近似模型为：
f(x)≈f(xk+1)+gTk+1(x−xk+1)+12(x−xk+1)TGk+1(x−xk+1)f(x)\approx f(x_{k+1})+g_{k+1}^T(x-x_{k+1})+\frac{1}{2}(x-x_{k+1})^TG_{k+1}(x-x_{k+1})
对上式求导得g(x)≈gk+1+Gk+1(x−xk+1)g(x)\approx g_{k+1}+G_{k+1}(x-x_{k+1})
令x=xkx=x_k，位移sk=xk+1−xks_k=x_{k+1}-x_k，梯度差yk=gk+1−gky_k=g_{k+1}-g_k，则有Gk+1sk≈ykG_{k+1}s_k\approx y_k.
因此，拟牛顿法中近似矩阵BkB_k要满足关系式Bk+1sk=ykB_{k+1}s_k=y_k
令Hk+1=B−1k+1H_{k+1}=B_{k+1}^{-1}，得到拟牛顿法的另一形式Hk+1yk=skH_{k+1}y_k=s_k
其中Hk+1H_{k+1}为海森矩阵逆矩阵的近似。上面两个公式称为拟牛顿方程。
搜索方向由dk=−Hkgkd_k=-H_kg_k或Bkdk=−gkB_kd_k=-g_k确定。根据近似矩阵的第三个特点，有
Bk+1=Bk+Ek,Hk+1=Hk+DkB_{k+1}=B_k+E_k,\qquad H_{k+1}=H_k+D_k
其中Ek,DkE_k,D_k为秩1或秩2矩阵。该公式称为校正规则。
通常将上面的三个式子（拟牛顿方程和校正规则）所确立的方法称为拟牛顿法。从下面的拟牛顿法的几个变种可以看出，不同的拟牛顿法的主要区别在于更新公式的不同。
Hk+1=Hk−HkykyTkHkyTkHkyk+sksTksTkykH_{k+1}=H_k-\frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\frac{s_ks_k^T}{s_k^Ty_k}
注意公式中的两个分数结构，分母yTkHkyky_k^TH_ky_k和sTkyks_k^Ty_k是标量，分子HkykyTkHkH_ky_ky_k^TH_k和sksTks_ks_k^T是与HkH_k同型的矩阵，而且都是对称矩阵。若HkH_k正定，且sTkyk>0s_k^Ty_k>0，则Hk+1H_{k+1}也正定。
当采用精确线搜索时，矩阵序列HkH_k的正定新条件sTkyk>0s_k^Ty_k>0可以被满足。但对于ArmijoArmijo搜索准则来说，不能满足这一条件，需要做如下修正：
Hk+1=⎧⎩⎨⎪⎪⎪⎪HkHk−HkykyTkHkyTkHkyk+sksTksTkyksTkyk≤0sTkyk>0⎫⎭⎬⎪⎪⎪⎪\begin{equation}
H_{k+1} = \left\{
\begin{array}{l l}
H_k & \quad \text{$s_k^Ty_k\leq 0$}\\
H_k-\frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\frac{s_ks_k^T}{s_k^Ty_k}& \quad \text{$s_k^Ty_k>0$}
\end{array} \right\}
\end{equation}
DFP算法描述：
step 1: 给定参数δ∈(0,1),σ∈(0,0.5)\delta\in(0,1),\sigma\in(0,0.5)，初始点x0∈Rx_0\in\mathbb{R},终止误差0≤ϵ≪10\leq \epsilon\ll1.初始对称正定矩阵H0H_0(通常取为G(x0)−1G(x_0)^{-1}或单位矩阵InI_n).令k←0k\leftarrow0
step 2: 计算gk=∇f(xk)g_k=\nabla f(x_k). 若||gk||≤ϵ||g_k||\leq\epsilon，停止迭代，输出x∗≈xkx^*\approx x_k
step 3: 计算搜索方向dk=−Hkgkd_k=-H_kg_k
step 4: 记mkm_k是满足下列不等式的最小非负整数mm.f(xk+βmdk)≤f(xk)+δβmgTkdkf(x_k+\beta^m d_k)\leq f(x_k)+\delta\beta^mg_k^Td_k令αk=δmk,xk+1=xk+αkdk\alpha_k=\delta^{m_k},x_{k+1}=x_k+\alpha_kd_k
step 5: 由校正公式确定Hk+1H_{k+1}
step 6： k←k+1k\leftarrow k+1,转 step 2
DFP算法的代码实现
由以上代码可知，拟牛顿法只需要在迭代开始前计算一次海森矩阵，更一般的，根本就不计算海森矩阵，而是初始化为单位矩阵，在每次迭代过程中是不需按传统方法（二阶导数）计算海森矩阵的。
实际性能
相关代码：
Bk+1=⎧⎩⎨⎪⎪⎪⎪BkBk−BkykyTkBkyTkBkyk+sksTksTkyksTkyk≤0sTkyk>0⎫⎭⎬⎪⎪⎪⎪\begin{equation}
B_{k+1} = \left\{
\begin{array}{l l}
B_k & \quad \text{$s_k^Ty_k\leq 0$}\\
B_k-\frac{B_ky_ky_k^TB_k}{y_k^TB_ky_k}+\frac{s_ks_k^T}{s_k^Ty_k}& \quad \text{$s_k^Ty_k>0$}
\end{array} \right\}
\end{equation}
BFGS算法的Python实现
实际性能
相关代码：
Hϕk+1=ϕkHBFGSk+1+(1−ϕk)HDFPk+1=Hk−HkykyTkHkyTkHkyk+sksTksTkyk+ϕkvkvTk\begin{equation}
\begin{array} {}
H_{k+1}^{\phi}=\phi_kH_{k+1}^{BFGS}+(1-\phi_k)H_{k+1}^{DFP}\\
\qquad= H_k-\frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\frac{s_ks_k^T}{s_k^Ty_k}+\phi_kv_kv_k^T
\end{array}
\end{equation}
其中ϕk∈[0,1]\phi_k\in[0,1]，vkv_k由下式定义：
vk=yTkHkyk−−−−−−−√(skyTksk−HkykyTkHkyk)v_k=\sqrt{y_k^TH_ky_k}\left(\frac{s_k}{y_k^Ts_k}-\frac{H_ky_k}{y_k^TH_ky_k}\right)
Broyden族算法Python实现
实际性能
相关代码：
基本BFGS算法的校正公式可以改写成
Hk+1=(I−skyTksTkyk)Hk(I−yksTksTkyk)+sksTksTkykH_{k+1}=\left(I-\frac{s_ky_k^T}{s_k^Ty_k}\right)H_k\left(I-\frac{y_ks_k^T}{s_k^Ty_k}\right)+\frac{s_ks_k^T}{s_k^Ty_k}
记ρk=1/sTkyk\rho_k = 1/s_k^Ty_k，以及Vk=(I−ρkyksTk)V_k=(I-\rho_ky_ks_k^T)，则上式可以写成
Hk+1=VTkHkVk+ρksksTkH_{k+1}=V_k^TH_kV_k+\rho_ks_ks_k^T
给定一个初始矩阵H0H_0(其取值后面有提到)，则由上式的递推关系可以得到
H1=VT0H0V0+ρ0s0sT0\begin{equation}
\begin{array} {}
H_1 = V_0^TH_0V_0+\rho_0s_0s_0^T\qquad\qquad\qquad\qquad\qquad
\end{array}
\end{equation}
H2=VT1H1V1+ρ1s1sT1=VT1(VT0H0V0+ρ0s0sT0)V1+ρ1s1sT1=VT1VT0H0V0V1+VT1ρ0s0sT0V1+ρ1s1sT1\begin{equation}
\begin{array} {}
H_2 = V_1^TH_1V_1+\rho_1s_1s_1^T\\
\qquad =V_1^T\left(V_0^TH_0V_0+\rho_0s_0s_0^T\right)V_1+\rho_1s_1s_1^T\\
\qquad =V_1^TV_0^TH_0V_0V_1+V_1^T\rho_0s_0s_0^TV_1+\rho_1s_1s_1^T
\end{array}
\end{equation}
⋯⋯⋯\cdots \cdots \cdots
更一般的，我们有
Hm+1=(VTKVTm−1⋯VT1VT0)H0(V0V1⋯Vm−1Vk)+(VTmVTm−1⋯VT2VT1)ρ0s0sT0(V1V2⋯Vm−1Vk)+(VTmVTm−1⋯VT3VT2)ρ1s1sT1(V2V3⋯Vm−1Vk)+⋯+(VTmVTm−1)ρm−2sm−2sTm−2(Vm−1Vk)VTkρm−1sm−1sTm−1Vk+ρmsmsTm
H_{m+1}=\left(V_K^TV_{m-1}^T\cdots V_1^TV_0^T\right) H_0 \left(V_0V_1\cdots V_{m-1}V_k\right)\\
\qquad+\left(V_{m}^TV_{m-1}^T\cdots V_2^TV_1^T\right)\rho_0s_0s_0^T\left(V_1V_2\cdots V_{m-1}V_k\right)\\
\qquad+\left(V_{m}^TV_{m-1}^T\cdots V_3^TV_2^T\right)\rho_1s_1s_1^T\left(V_2V_3\cdots V_{m-1}V_k\right)\\
\qquad + \cdots\qquad \qquad \qquad\qquad \qquad\qquad\qquad\qquad \qquad\\
\qquad+\left(V_{m}^TV_{m-1}^T\right)\rho_{m-2}s_{m-2}s_{m-2}^T\left(V_{m-1}V_k\right)\qquad\qquad\qquad\\
\qquad V_k^T\rho_{m-1}s_{m-1}s_{m-1}^TV_k\\
\qquad+\rho_ms_ms_m^T\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
在上式的右边中，H0H_0是由我们设置的，其余变量可由保存的最近mm次迭代所形成的向量序列，如位移向量序列{sk}\{s_k\}和一阶导数差所形成的向量序列{yk}\{y_k\}获得。也就是说，可由最近mm次迭代的信息计算当前的海森矩阵的近似矩阵。
补充几点：
H0=I⋅sTmymyTmymH_0=I\cdot \frac{s_m^Ty_m}{y_m^Ty_m}，第一次迭代时，序列{sks_k}和{yky_k}为空，则H0=IH_0=I
最初的若干次迭代， 序列{sks_k}和{yky_k}的元素个数较小，会有n<mn，则将Hm+1H_{m+1}计算公式右边的mm换成nn即可。
随着迭代次数增加， 序列{sks_k}和{yky_k}的元素个数增大，会有n>mn>m。由于我们只需要最近mm次迭代产生的序列元素，所以序列{sks_k}和{yky_k}只需要保存最新的mm个元素即可，如果再有新的元素进入，则同时扔掉最旧的元素，以保证序列元素个数恒为mm。
这就是L-BFGS算法的思想。聪明的同学会问：你上面的公式不还是要计算海森矩阵的近似矩阵吗？那岂不是还是需要n2n^2规模的内存？
其实在实际中是不计算该矩阵的, 而且不是采用上面的方法，而是直接得到HkgkH_kg_k，而搜索方向就是dk=−Hkgkd_k=-H_kg_k。下面给出了计算的函数twoloop(s,y,ρ\rho,gk)的伪代码，可知该函数内部的计算仅限于标量和向量，未出现矩阵。
函数参数s,ys,y即向量序列{sks_k}和{yky_k}，ρ\rho为元素序列{ρk\rho_k}，其元素ρk=1/sTkyk\rho_k = 1/s_k^Ty_k，gkg_k是向量，为当前的一阶导数，输出为向量z=Hkgkz=H_kg_k，即搜索方向的反方向
Functiontwoloop(s,y,ρ,gk)q=gkFori=k−1,k−2,…,k−mαi=ρisTiqq=q−αiyiHk=yTk−1sk−1/yTk−1yk−1z=HkqFori=k−m,k−m+1,…,k−1βi=ρiyTizz=z+si(αi−βi)Returnz
Function\,\,twoloop(s, y, \rho,gk)\qquad\qquad\qquad\qquad\\
q = g_k\,\, \qquad\qquad\qquad\qquad\qquad\qquad\quad\\
For \,i=k-1, k-2, \ldots, k-m\qquad\quad\\
\alpha_i = \rho_i s^{\rm T}_i q\,\!\qquad\qquad\quad\\
q = q - \alpha_i y_i\,\!\qquad\qquad\\
H_k=y^{\rm T}_{k-1} s_{k-1}/y^{\rm T}_{k-1} y_{k-1}\qquad \qquad\quad\\
z = H_k q\qquad\qquad\qquad\qquad\qquad\qquad\\
For\,\, i=k-m, k-m+1, \ldots, k-1\\
\beta_i = \rho_i y^{\rm T}_i z\,\!\qquad\qquad\\
z = z + s_i (\alpha_i - \beta_i)\,\!\\
Return \,\, z\qquad\qquad\qquad\qquad\qquad\qquad\,\!
给出L-BFGS的算法
可以看到其搜索方向dkd_k是根据函数twolooptwoloop计算得到的。
L-BFGS算法Python实现
实际性能
相关代码：
http://papers.nips.cc/paper/5333-large-scale-l-bfgs-using-mapreduce.pdf
L-BFGS的原理及在回归分析中的应用
http://blog.csdn.net/hlx371240/article/details/39970727
L-BFGS的原理和MATLAB实现
deep learning Softmax分类器（L-BFGS,CG,SD）
http://blog.csdn.net/hlx371240/article/details/40015395
机器学习中导数最优化方法
http://www.cnblogs.com/daniel-D/p/3377840.html
Newton’s method in optimization
http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization
Quasi-Newton method
http://en.wikipedia.org/wiki/Quasi-Newton_method
Davidon–Fletcher–Powell formula
http://en.wikipedia.org/wiki/Davidon%E2%80%93Fletcher%E2%80%93Powell_formula
Broyden–Fletcher–Goldfarb–Shanno algorithm
http://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm
Limited-memory BFGS
http://en.wikipedia.org/wiki/Limited-memory_BFGS
some books
《Basic Concepts and Stationary Iterative Methods》
《Iterative Methods for Optimization》
《Non-Linear Least-Squares Minimization》
some papers
《A comparison of algorithms for maximum entropy parameter estimation》
《Updating Quasi-Newton Matrices with Limited Storage》
《Scalable Training of L1-Regularized Log-Linear Models》
《最优化方法及其MATLAB程序设计》，马昌凤，科学出版社
《优化方法》，李春明，东南大学出版社
《应用最优化方法及MATLAB实现》，刘兴高，胡云卿，科学出版社
《优化技术与MATLAB优化工具箱》，赵继俊，机械工业出版社
http://hub.darcs.net/amgarching/pts/browse/bfgs.py
markdown 公式，可直接复制得到特性公式
http://www.aleacubase.com/cudalab/cudalab_usage-math_formatting_on_markdown.html
CSDN-markdown语法之如何插入图片
/article/1369460.html
要求解的问题
线搜索技术和Armijo准则
最速下降法及其Python实现
牛顿法
阻尼牛顿法及其Python实现
修正牛顿法法及其Python实现
拟牛顿法
DFP算法及其Python实现
BFGS算法及其Python实现
Broyden族算法及其Python实现
L-BFGS算法及其Python实现
参考文献
1.要求解的问题
求解无约束非线性优化问题minx∈R2f(x)=100(x21−x2)2+(x1−1)2\min_{x\in\mathbb{R}^2}f(x)=100(x_1^2-x_2)^2+(x_1-1)^2
该问题有精确解x∗=(1,1)T,f(x∗)=0x^*=(1,1)^T,f(x^*)=0
其梯度向量g(x)=(400x1(x21−x2)+2(x1−1),−200(x21−x2))g(x)=(400x_1(x_1^2-x_2)+2(x_1-1),-200(x_1^2-x_2))
其海森矩阵G(x)=[1200x21−400x2+2−400x1−400x1200]G(x) = \begin{bmatrix}
1200x_1^2-400x_2+2&-400x_1 \\[0.3em]
-400x_1&200
\end{bmatrix}
下面定义了三个Python语言编写的函数：函数表达式fun，梯度向量gfun，和海森矩阵hess。这三个表达式在后面各个算法的实现中会用到。
# 函数表达式fun fun = lambda x:100*(x[0]**2-x[1])**2 + (x[0]-1)**2 # 梯度向量 gfun gfun = lambda x:np.array([400*x[0]*(x[0]**2-x[1])+2*(x[0]-1), -200*(x[0]**2-x[1])]) # 海森矩阵 hess hess = lambda x:np.array([[1200*x[0]**2-400*x[1]+2, -400*x[0]],[-400*x[0],200]])
2.线搜索技术和Armijo准则
线搜索技术是求解许多优化问题下降算法的基本组成部分，但精确线性搜索需要计算很多的函数值和梯度值，从而耗费较多资源。特别是当迭代点远离最优点时，精确线搜索技术通常不是有效和合理的。对于许多优化算法，其收敛速度并不依赖于精确搜索过程。因此既能保证目标函数具有可接受的下降量，又能使最终形成的迭代序列收敛的非精确先搜索越来越流行。符号约定:
gkg_k : ∇f(xk)\nabla f(x_k)，即第目标函数关于kk次迭代值xkx_k的导数。
GkG_k : G(xk)=∇2f(xk)G(x_k)=\nabla^2f(x_k)，即海森矩阵。
dkd_k : 第kk次迭代的优化方向。在最速下降算法中,有dk=−gkd_k = -g_k
αk\alpha_k : 第kk次迭代的步长因子，有xk+1=xk+αkdkx_{k+1}=x_k+\alpha_k d_k
在精确线性搜索中，步长因子αk\alpha_k由下面的式子确定：
αk=argminαf(xk+αdk)\alpha_k=\arg\min_{\alpha} f(x_k+\alpha d_k)
而对于非精确线性搜索，选取的αk\alpha_k只要使得目标函数ff得到可接受的下降量，即
Δfk=f(xk)−f(xk+αkdk)>0\Delta f_k = f(x_k)-f(x_k+\alpha_kd_k)>0
Armijo准则用于非精确线性搜索中步长因子α\alpha的确定，内容如下：
Armijo准则：
已知当前位置xkx_k和优化的方向dkd_k ，参数β∈(0,1)\beta\in(0,1),δ∈(0,0.5)\delta\in(0,0.5).令步长因子αk=βmk\alpha_k=\beta^{m_k},其中mkm_k为满足下列不等式的最小非负整数mm：
f(xk+βmdk)≤f(xk)+δβmgTkdkf(x_k+\beta^m d_k)\leq f(x_k)+\delta\beta^mg_k^Td_k
由此确定的下一个位置xk+1=xk+αkdkx_{k+1}=x_k+\alpha_k d_k
该准则在接下来的介绍的几个算法中多次使用。
3.最速下降法及其Python实现
算法描述：step 1 :选取初始点x0∈Rnx_0\in\mathbb{R}^n,容许误差0≤ϵ≪10\leq\epsilon\ll 1 .令 k←1k\leftarrow1.
step 2 :计算gk=∇f(xk)g_k=\nabla f(x_k). 若||gk||≤ϵ||g_k||\leq\epsilon，停止迭代，输出xkx_k作为近似最优解。
step 3 :取方向dk=−gkd_k=-g_k.
step 4 :由线搜索技术确定步长因子αk\alpha_k.
step 5 :令xk+1←xk+αkdkx_{k+1}\leftarrow x_k+\alpha_kd_k，转step 2.
上述step 3中步长因子αk\alpha_k的确定既可以使用精确线搜索方法，也可以使用非精确线搜索方法，在理论上都能保证其全局收敛性。若采用精确线搜索方法，则αk\alpha_k满足
ddαf(xk+αdk)|α=αk=∇f(xk+αkdk)Tdk=0\frac{\mathrm{d}}{\mathrm{d}\alpha}f(x_k+\alpha d_k)|_{\alpha=\alpha_k}=\nabla f(x_k+\alpha_kd_k)^Td_k=0
从而有∇f(xk+1)T∇f(xk)=0\nabla f(x_{k+1})^T\nabla f(x_k)=0
即xk+1x_{k+1}处的梯度与其前驱xkx_k处的梯度是正交的，也就是说，迭代点序列所走的路线是锯齿形的，这造成其收敛速度很缓慢。如下图所示，其中绿色折线所显示的路线就是由最速下降法得到的，红色曲线是由共轭梯度下降法确定的，通过对比可以看出梯度下降法所走的路线是锯齿形的，经测试，其收敛速度非常慢。
最速下降法的Python实现。
def gradient(fun,gfun,x0):
#最速下降法求解无约束问题
# x0是初始点，fun和gfun分别是目标函数和梯度
maxk = 5000
rho = 0.5
sigma = 0.4
k = 0
epsilon = 1e-5
while k<maxk:
gk = gfun(x0)
dk = -gk
if np.linalg.norm(dk) < epsilon:
break
m = 0
mk = 0
while m< 20:
if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):
mk = m
break
m += 1
x0 += rho**mk*dk
k += 1
return x0,fun(x0),k性能测试
可以看到大约有一半的样本的迭代次数要超过2000次。
4.牛顿法
牛顿法的基本思想是用迭代点xkx_k处的一阶导数（梯度gkg_k）和二阶倒数（海森矩阵GkG_k）对目标函数进行二次函数近似，然后把二次模型的极小点作为新的迭代点。牛顿法推导
函数f(x)f(x)在xkx_k处的泰勒展开式的前三项得
T(f,xk,3)=fk+gTk(x−xk)+12(x−xk)TGk(x−xk)T(f,x_k,3)=f_k+g_k^T(x-x_k)+\frac{1}{2}(x-x_k)^TG_k(x-x_k)
则其稳定点∇T=gk+Gk(x−xk)=0\nabla T = g_k+G_k(x-x_k)=0
若GkG_k非奇异，那么解上面的线性方程（标记其解为xk+1x_{k+1}）即得到牛顿法的迭代公式xk+1=xk−G−1kgkx_{k+1}=x_k-G_k^{-1}g_k
即优化方向为dk=−G−1kgkd_k=-G_k^{-1}g_k
求稳定点是用到了以下公式：
考虑y=xTAxy=x^TAx，根据矩阵求导法则，有
dydx=Ax+ATxd2ydx2=A+AT\frac{\mathrm{d}y}{\mathrm{d}x}=Ax+A^Tx\\\frac{\mathrm{d}^2y}{\mathrm{d}x^2}=A+A^T
注意实际中dkd_k的是通过求解线性方程组Gkd=−gkG_kd=-g_k获得的。
阻尼牛顿法及其Python实现
阻尼牛顿法采用了牛顿法的思想，唯一的差别是阻尼牛顿法并不直接采用xk+1=xk+dkx_{k+1}=x_k+d_k，而是用线搜索技术获得合适的步长因子αk\alpha_k,用公式xk+1=xk+δmdkx_{k+1}=x_k+\delta^md_k计算xk+1x_{k+1}。阻尼牛顿法的算法描述：
step 1: 给定终止误差0≤ϵ≪1,δ∈(0,1),σ∈(0,0.5)0\leq\epsilon\ll 1,\delta\in(0,1),\sigma\in(0,0.5). 初始点x0∈Rnx_0\in\mathbb{R}^n. 令k←0k\leftarrow0
step 2: 计算gk=∇f(xk)g_k=\nabla f(x_k). 若||gk||≤ϵ||g_k||\leq\epsilon，停止迭代，输出x∗≈xkx^*\approx x_k
step 3: 计算Gk=∇2f(xk)G_k=\nabla^2f(x_k)，并求解线性方程组得到解dkd_k,Gkd=−gkG_kd=-g_k
step 4: 记mkm_k是满足下列不等式的最小非负整数mm.f(xk+βmdk)≤f(xk)+δβmgTkdkf(x_k+\beta^m d_k)\leq f(x_k)+\delta\beta^mg_k^Td_k
step 5: 令αk=δmk,xk+1=xk+αkdk,k←k+1\alpha_k=\delta^{m_k},x_{k+1}=x_k+\alpha_kd_k,k\leftarrow k+1，转step 2
阻尼牛顿法的Python实现：
def dampnm(fun,gfun,hess,x0):
# 用阻尼牛顿法求解无约束问题
#x0是初始点，fun，gfun和hess分别是目标函数值，梯度，海森矩阵的函数
maxk = 500
rho = 0.55
sigma = 0.4
k = 0
epsilon = 1e-5
while k < maxk:
gk = gfun(x0)
Gk = hess(x0)
dk = -1.0*np.linalg.solve(Gk,gk)
if np.linalg.norm(dk) < epsilon:
break
m = 0
mk = 0
while m < 20:
if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):
mk = m
break
m += 1
x0 += rho**mk*dk
k += 1
return x0,fun(x0),k性能展示
作图代码：
iters = []
for i in xrange(np.shape(data)[0]):
rt = dampnm(fun,gfun,hess,data[i])
if rt[2] <=200:
iters.append(rt[2])
if i%100 == 0:
print i,rt[2]
plt.hist(iters,bins=50)
plt.title(u'阻尼牛顿法迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})修正牛顿法法及其Python实现
牛顿法要求目标函数的海森矩阵G(x)=∇2f(x)G(x)=\nabla ^2f(x)在每个迭代点xkx_k处是正定的，否则难以保证牛顿方向dk=−G−1kgkd_k=-G_k^{-1}g_k是ff在xkx_k处的下降方向。为克服这一缺陷，需要对牛顿法进行修正。下面介绍两种修正方法，分别是“牛顿-最速下降混合算法”和“修正牛顿法”。
牛顿-最速下降混合算法
该方法将牛顿法和最速下降法结合起来，基本思想是：当海森矩阵正定时，采用牛顿法确定的优化方向作为搜索方向；否则，即海森矩阵为非正定矩阵时，则采用负梯度方向作为搜索方向。
修正牛顿法
引入阻尼因子μk≥0\mu_k\geq0，即在每一迭代步适当选取参数μk\mu_k，使得矩阵Ak=Gk+μkIA_k=G_k+\mu_k I正定，用AkA_k代替GkG_k来求解dkd_k。
算法描述：
step 1: 给定终止误差0≤ϵ≪1,δ∈(0,1),σ∈(0,0.5)0\leq\epsilon\ll 1,\delta\in(0,1),\sigma\in(0,0.5). 初始点x0∈Rnx_0\in\mathbb{R}^n. 令k←0k\leftarrow0
step 2: 计算gk=∇f(xk)g_k=\nabla f(x_k)，μk=||gk||1+τ\mu_k = ||g_k||^{1+\tau}. 若||gk||≤ϵ||g_k||\leq \epsilon，停止迭代，输出x∗≈xkx^*\approx x_k
step 3: 计算Gk=∇2f(xk)G_k=\nabla^2f(x_k)，并求解线性方程组得到解dkd_k,(Gk+μkI)d=−gk(G_k+\mu_k I)d=-g_k
step 4: 记mkm_k是满足下列不等式的最小非负整数mm.f(xk+βmdk)≤f(xk)+δβmgTkdkf(x_k+\beta^m d_k)\leq f(x_k)+\delta\beta^mg_k^Td_k
step 5: 令αk=δmk,xk+1=xk+αkdk,k←k+1\alpha_k=\delta^{m_k},x_{k+1}=x_k+\alpha_kd_k,k\leftarrow k+1，转step 2
修正牛顿法的Python实现代码：
def revisenm(fun,gfun,hess,x0):
# 用修正牛顿法求解无约束问题
#x0是初始点，fun，gfun和hess分别是目标函数值，梯度，海森矩阵的函数
maxk = 1e5
n = np.shape(x0)[0]
rho = 0.55
sigma = 0.4
tau = 0.0
k = 0
epsilon = 1e-5
while k < maxk:
gk = gfun(x0)
if np.linalg.norm(gk) < epsilon:
break
muk = np.power(np.linalg.norm(gk),1+tau)
Gk = hess(x0)
Ak = Gk + muk*np.eye(n)
dk = -1.0*np.linalg.solve(Ak,gk)
m = 0
mk = 0
while m < 20:
if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):
mk = m
break
m += 1
x0 += rho**mk*dk
k += 1
return x0,fun(x0),k性能展示
5.拟牛顿法
牛顿法的优点是具有二阶收敛速度，缺点是：但当海森矩阵G(xk)=∇2f(x)G(x_k)=\nabla^2f(x) 不正定时，不能保证所产生的方向是目标函数在xkx_k处的下降方向。
特别地，当G(xk)G(x_k)奇异时，算法就无法继续进行下去。尽管修正牛顿法可以克服这一缺陷，但修正参数的取值很难把握，过大或过小都会影响到收敛速度。
牛顿法的每一步迭代都需要目标函数的海森矩阵G(xk)G(x_k)，对于大规模问题其计算量是惊人的。
拟牛顿法的基本思想是用海森矩阵GkG_k的某个近似矩阵BkB_k取代GkG_k. BkB_k通常具有下面三个特点：
在某种意义下有Bk≈GkB_k\approx G_k ,使得相应的算法产生的方向近似于牛顿方向，确保算法具有较快的收敛速度。
对所有的kk，BkB_k是正定的，从而使得算法所产生的方向是函数ff在xkx_k处下降方向。
矩阵BkB_k更新规则比较简单
设 f:Rn→Rf:\mathbb{R}^n\rightarrow\mathbb{R}在开集D⊂RnD\subset \mathbb{R}^n上二次连续可微，那么ff在xk+1x_{k+1}处的二次近似模型为：
f(x)≈f(xk+1)+gTk+1(x−xk+1)+12(x−xk+1)TGk+1(x−xk+1)f(x)\approx f(x_{k+1})+g_{k+1}^T(x-x_{k+1})+\frac{1}{2}(x-x_{k+1})^TG_{k+1}(x-x_{k+1})
对上式求导得g(x)≈gk+1+Gk+1(x−xk+1)g(x)\approx g_{k+1}+G_{k+1}(x-x_{k+1})
令x=xkx=x_k，位移sk=xk+1−xks_k=x_{k+1}-x_k，梯度差yk=gk+1−gky_k=g_{k+1}-g_k，则有Gk+1sk≈ykG_{k+1}s_k\approx y_k.
因此，拟牛顿法中近似矩阵BkB_k要满足关系式Bk+1sk=ykB_{k+1}s_k=y_k
令Hk+1=B−1k+1H_{k+1}=B_{k+1}^{-1}，得到拟牛顿法的另一形式Hk+1yk=skH_{k+1}y_k=s_k
其中Hk+1H_{k+1}为海森矩阵逆矩阵的近似。上面两个公式称为拟牛顿方程。
搜索方向由dk=−Hkgkd_k=-H_kg_k或Bkdk=−gkB_kd_k=-g_k确定。根据近似矩阵的第三个特点，有
Bk+1=Bk+Ek,Hk+1=Hk+DkB_{k+1}=B_k+E_k,\qquad H_{k+1}=H_k+D_k
其中Ek,DkE_k,D_k为秩1或秩2矩阵。该公式称为校正规则。
通常将上面的三个式子（拟牛顿方程和校正规则）所确立的方法称为拟牛顿法。从下面的拟牛顿法的几个变种可以看出，不同的拟牛顿法的主要区别在于更新公式的不同。
DFP算法及其Python实现
DFP算法的校正公式如下：Hk+1=Hk−HkykyTkHkyTkHkyk+sksTksTkykH_{k+1}=H_k-\frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\frac{s_ks_k^T}{s_k^Ty_k}
注意公式中的两个分数结构，分母yTkHkyky_k^TH_ky_k和sTkyks_k^Ty_k是标量，分子HkykyTkHkH_ky_ky_k^TH_k和sksTks_ks_k^T是与HkH_k同型的矩阵，而且都是对称矩阵。若HkH_k正定，且sTkyk>0s_k^Ty_k>0，则Hk+1H_{k+1}也正定。
当采用精确线搜索时，矩阵序列HkH_k的正定新条件sTkyk>0s_k^Ty_k>0可以被满足。但对于ArmijoArmijo搜索准则来说，不能满足这一条件，需要做如下修正：
Hk+1=⎧⎩⎨⎪⎪⎪⎪HkHk−HkykyTkHkyTkHkyk+sksTksTkyksTkyk≤0sTkyk>0⎫⎭⎬⎪⎪⎪⎪\begin{equation}
H_{k+1} = \left\{
\begin{array}{l l}
H_k & \quad \text{$s_k^Ty_k\leq 0$}\\
H_k-\frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\frac{s_ks_k^T}{s_k^Ty_k}& \quad \text{$s_k^Ty_k>0$}
\end{array} \right\}
\end{equation}
DFP算法描述：
step 1: 给定参数δ∈(0,1),σ∈(0,0.5)\delta\in(0,1),\sigma\in(0,0.5)，初始点x0∈Rx_0\in\mathbb{R},终止误差0≤ϵ≪10\leq \epsilon\ll1.初始对称正定矩阵H0H_0(通常取为G(x0)−1G(x_0)^{-1}或单位矩阵InI_n).令k←0k\leftarrow0
step 2: 计算gk=∇f(xk)g_k=\nabla f(x_k). 若||gk||≤ϵ||g_k||\leq\epsilon，停止迭代，输出x∗≈xkx^*\approx x_k
step 3: 计算搜索方向dk=−Hkgkd_k=-H_kg_k
step 4: 记mkm_k是满足下列不等式的最小非负整数mm.f(xk+βmdk)≤f(xk)+δβmgTkdkf(x_k+\beta^m d_k)\leq f(x_k)+\delta\beta^mg_k^Td_k令αk=δmk,xk+1=xk+αkdk\alpha_k=\delta^{m_k},x_{k+1}=x_k+\alpha_kd_k
step 5: 由校正公式确定Hk+1H_{k+1}
step 6： k←k+1k\leftarrow k+1,转 step 2
DFP算法的代码实现
def dfp(fun,gfun,hess,x0):
#功能：用DFP族算法求解无约束问题：min fun(x)
#输入：x0是初始点，fun,gfun分别是目标函数和梯度
#输出：x,val分别是近似最优点和最优解,k是迭代次数
maxk = 1e5
rho = 0.55
sigma = 0.4
epsilon = 1e-5
k = 0
n = np.shape(x0)[0]
#海森矩阵可以初始化为单位矩阵
Hk = np.linalg.inv(hess(x0)) #或者单位矩阵np.eye(n)
while k < maxk :
gk = gfun(x0)
if np.linalg.norm(gk) < epsilon:
break
dk = -1.0*np.dot(Hk,gk)
m=0;
mk=0
while m < 20: # 用Armijo搜索求步长
if fun(x0+rho**m*dk) < fun(x0)+sigma*rho**m*np.dot(gk,dk):
mk = m
break
m += 1
#print mk
#DFP校正
x = x0 + rho**mk*dk
sk = x - x0
yk = gfun(x) - gk
if np.dot(sk,yk) > 0:
Hy = np.dot(Hk,yk)
print Hy
sy = np.dot(sk,yk) #向量的点积
yHy = np.dot(np.dot(yk,Hk),yk) # yHy是标量
#表达式Hy.reshape((n,1))*Hy 中Hy是向量，生成矩阵
Hk = Hk - 1.0*Hy.reshape((n,1))*Hy/yHy + 1.0*sk.reshape((n,1))*sk/sy
k += 1
x0 = x
return x0,fun(x0),k#分别是最优点坐标，最优值，迭代次数由以上代码可知，拟牛顿法只需要在迭代开始前计算一次海森矩阵，更一般的，根本就不计算海森矩阵，而是初始化为单位矩阵，在每次迭代过程中是不需按传统方法（二阶导数）计算海森矩阵的。
实际性能
相关代码：
n = 50
x = y = np.linspace(-10,10,n)
xx,yy = np.meshgrid(x,y)
data = [[xx[i][j],yy[i][j]] for i in range(n) for j in range(n)]
iters = []
for i in xrange(np.shape(data)[0]):
rt = dfp(fun,gfun,hess,np.array(data[i]))
if rt[2] <=200: # 提出迭代次数过大的
iters.append(rt[2])
if i%100 == 0:
print i,rt[2]
plt.hist(iters,bins=50)
plt.title(u'DFP迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})BFGS算法及其Python实现
BFGS算法与DFP步骤基本相同，区别在于更新公式的差异Bk+1=⎧⎩⎨⎪⎪⎪⎪BkBk−BkykyTkBkyTkBkyk+sksTksTkyksTkyk≤0sTkyk>0⎫⎭⎬⎪⎪⎪⎪\begin{equation}
B_{k+1} = \left\{
\begin{array}{l l}
B_k & \quad \text{$s_k^Ty_k\leq 0$}\\
B_k-\frac{B_ky_ky_k^TB_k}{y_k^TB_ky_k}+\frac{s_ks_k^T}{s_k^Ty_k}& \quad \text{$s_k^Ty_k>0$}
\end{array} \right\}
\end{equation}
BFGS算法的Python实现
def bfgs(fun,gfun,hess,x0):
#功能：用BFGS族算法求解无约束问题：min fun(x)
#输入：x0是初始点，fun,gfun分别是目标函数和梯度
#输出：x,val分别是近似最优点和最优解,k是迭代次数
maxk = 1e5
rho = 0.55
sigma = 0.4
epsilon = 1e-5
k = 0
n = np.shape(x0)[0]
#海森矩阵可以初始化为单位矩阵
Bk = eye(n)#np.linalg.inv(hess(x0)) #或者单位矩阵np.eye(n)
while k < maxk:
gk = gfun(x0)
if np.linalg.norm(gk) < epsilon:
break
dk = -1.0*np.linalg.solve(Bk,gk)
m = 0
mk = 0
while m < 20: # 用Armijo搜索求步长
if fun(x0+rho**m*dk) < fun(x0)+sigma*rho**m*np.dot(gk,dk):
mk = m
break
m += 1
#BFGS校正
x = x0 + rho**mk*dk
sk = x - x0
yk = gfun(x) - gk
if np.dot(sk,yk) > 0:
Bs = np.dot(Bk,sk)
ys = np.dot(yk,sk)
sBs = np.dot(np.dot(sk,Bk),sk)
Bk = Bk - 1.0*Bs.reshape((n,1))*Bs/sBs + 1.0*yk.reshape((n,1))*yk/ys
k += 1
x0 = x
return x0,fun(x0),k#分别是最优点坐标，最优值，迭代次数实际性能
相关代码：
iters = []
for i in xrange(np.shape(data)[0]):
rt = bfgs(fun,gfun,hess,np.array(data[i]))
if rt[2] <=200:
iters.append(rt[2])
if i%100 == 0:
print i,rt[2]
plt.hist(iters,bins=50)
plt.title(u'BFGS迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})Broyden族算法及其Python实现
之前的BFGS和DFP校正都是由yky_k和BkskB_ks_k(或 sks_k和HkykH_ky_k)组成的秩2矩阵。而Broyden族算法采用了BFGS和DFP校正公式的凸组合，如下：Hϕk+1=ϕkHBFGSk+1+(1−ϕk)HDFPk+1=Hk−HkykyTkHkyTkHkyk+sksTksTkyk+ϕkvkvTk\begin{equation}
\begin{array} {}
H_{k+1}^{\phi}=\phi_kH_{k+1}^{BFGS}+(1-\phi_k)H_{k+1}^{DFP}\\
\qquad= H_k-\frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\frac{s_ks_k^T}{s_k^Ty_k}+\phi_kv_kv_k^T
\end{array}
\end{equation}
其中ϕk∈[0,1]\phi_k\in[0,1]，vkv_k由下式定义：
vk=yTkHkyk−−−−−−−√(skyTksk−HkykyTkHkyk)v_k=\sqrt{y_k^TH_ky_k}\left(\frac{s_k}{y_k^Ts_k}-\frac{H_ky_k}{y_k^TH_ky_k}\right)
Broyden族算法Python实现
def broyden(fun,gfun,hess,x0):
#功能：用Broyden族算法求解无约束问题：min fun(x)
#输入：x0是初始点，fun,gfun分别是目标函数和梯度
#输出：x,val分别是近似最优点和最优解,k是迭代次数
x0 = np.array(x0)
maxk = 1e5
rho = 0.55;
sigma = 0.4;
epsilon = 1e-5
phi = 0.5;
k=0;
n = np.shape(x0)[0]
Hk = np.linalg.inv(hess(x0))
while k<maxk :
gk = gfun(x0)
if np.linalg.norm(gk) < epsilon:
break
dk = -1*np.dot(Hk,gk)
m=0;mk=0
while m < 20: # 用Armijo搜索求步长
if fun(x0+rho**m*dk) < fun(x0)+sigma*rho**m*np.dot(gk,dk):
mk = m
break
m += 1
#Broyden族校正
x = x0 + rho**mk*dk
sk = x - x0
yk = gfun(x) - gk
Hy = np.dot(Hk,yk)
sy = np.dot(sk,yk)
yHy = np.dot(np.dot(yk,Hk),yk)
if(sy < 0.2 *yHy):
theta = 0.8*yHy/(yHy-sy)
sk = theta*sk + (1-theta)*Hy
sy = 0.2*yHy
vk = np.sqrt(yHy)*(sk/sy-Hy/yHy)
Hk = Hk - Hy.reshape((n,1))*Hy/yHy +sk.reshape((n,1))*sk/sy + phi*vk.reshape((n,1))*vk
k += 1
x0 = x
return x0,fun(x0),k #分别是最优点坐标，最优值，迭代次数实际性能
相关代码：
iters = []
for i in xrange(np.shape(data)[0]):
rt = broyden(fun,gfun,hess,np.array(data[i]))
if rt[2] <=200:
iters.append(rt[2])
if i%100 == 0:
print i,rt[2]
plt.hist(iters,bins=50)
plt.title(u'Broyden迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})L-BFGS算法及其Python实现
拟牛顿法(如BFGS算法)需要计算和存储海森矩阵，其空间复杂度是n2n^2，当nn很大时，其需要的内存量是非常大的。为了解决该内存问题，有限内存BFGS（即传说中的L-BFGS算法）横空出世。基本BFGS算法的校正公式可以改写成
Hk+1=(I−skyTksTkyk)Hk(I−yksTksTkyk)+sksTksTkykH_{k+1}=\left(I-\frac{s_ky_k^T}{s_k^Ty_k}\right)H_k\left(I-\frac{y_ks_k^T}{s_k^Ty_k}\right)+\frac{s_ks_k^T}{s_k^Ty_k}
记ρk=1/sTkyk\rho_k = 1/s_k^Ty_k，以及Vk=(I−ρkyksTk)V_k=(I-\rho_ky_ks_k^T)，则上式可以写成
Hk+1=VTkHkVk+ρksksTkH_{k+1}=V_k^TH_kV_k+\rho_ks_ks_k^T
给定一个初始矩阵H0H_0(其取值后面有提到)，则由上式的递推关系可以得到
H1=VT0H0V0+ρ0s0sT0\begin{equation}
\begin{array} {}
H_1 = V_0^TH_0V_0+\rho_0s_0s_0^T\qquad\qquad\qquad\qquad\qquad
\end{array}
\end{equation}
H2=VT1H1V1+ρ1s1sT1=VT1(VT0H0V0+ρ0s0sT0)V1+ρ1s1sT1=VT1VT0H0V0V1+VT1ρ0s0sT0V1+ρ1s1sT1\begin{equation}
\begin{array} {}
H_2 = V_1^TH_1V_1+\rho_1s_1s_1^T\\
\qquad =V_1^T\left(V_0^TH_0V_0+\rho_0s_0s_0^T\right)V_1+\rho_1s_1s_1^T\\
\qquad =V_1^TV_0^TH_0V_0V_1+V_1^T\rho_0s_0s_0^TV_1+\rho_1s_1s_1^T
\end{array}
\end{equation}
⋯⋯⋯\cdots \cdots \cdots
更一般的，我们有
Hm+1=(VTKVTm−1⋯VT1VT0)H0(V0V1⋯Vm−1Vk)+(VTmVTm−1⋯VT2VT1)ρ0s0sT0(V1V2⋯Vm−1Vk)+(VTmVTm−1⋯VT3VT2)ρ1s1sT1(V2V3⋯Vm−1Vk)+⋯+(VTmVTm−1)ρm−2sm−2sTm−2(Vm−1Vk)VTkρm−1sm−1sTm−1Vk+ρmsmsTm
H_{m+1}=\left(V_K^TV_{m-1}^T\cdots V_1^TV_0^T\right) H_0 \left(V_0V_1\cdots V_{m-1}V_k\right)\\
\qquad+\left(V_{m}^TV_{m-1}^T\cdots V_2^TV_1^T\right)\rho_0s_0s_0^T\left(V_1V_2\cdots V_{m-1}V_k\right)\\
\qquad+\left(V_{m}^TV_{m-1}^T\cdots V_3^TV_2^T\right)\rho_1s_1s_1^T\left(V_2V_3\cdots V_{m-1}V_k\right)\\
\qquad + \cdots\qquad \qquad \qquad\qquad \qquad\qquad\qquad\qquad \qquad\\
\qquad+\left(V_{m}^TV_{m-1}^T\right)\rho_{m-2}s_{m-2}s_{m-2}^T\left(V_{m-1}V_k\right)\qquad\qquad\qquad\\
\qquad V_k^T\rho_{m-1}s_{m-1}s_{m-1}^TV_k\\
\qquad+\rho_ms_ms_m^T\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
在上式的右边中，H0H_0是由我们设置的，其余变量可由保存的最近mm次迭代所形成的向量序列，如位移向量序列{sk}\{s_k\}和一阶导数差所形成的向量序列{yk}\{y_k\}获得。也就是说，可由最近mm次迭代的信息计算当前的海森矩阵的近似矩阵。
补充几点：
H0=I⋅sTmymyTmymH_0=I\cdot \frac{s_m^Ty_m}{y_m^Ty_m}，第一次迭代时，序列{sks_k}和{yky_k}为空，则H0=IH_0=I
最初的若干次迭代， 序列{sks_k}和{yky_k}的元素个数较小，会有n<mn，则将Hm+1H_{m+1}计算公式右边的mm换成nn即可。
随着迭代次数增加， 序列{sks_k}和{yky_k}的元素个数增大，会有n>mn>m。由于我们只需要最近mm次迭代产生的序列元素，所以序列{sks_k}和{yky_k}只需要保存最新的mm个元素即可，如果再有新的元素进入，则同时扔掉最旧的元素，以保证序列元素个数恒为mm。
这就是L-BFGS算法的思想。聪明的同学会问：你上面的公式不还是要计算海森矩阵的近似矩阵吗？那岂不是还是需要n2n^2规模的内存？
其实在实际中是不计算该矩阵的, 而且不是采用上面的方法，而是直接得到HkgkH_kg_k，而搜索方向就是dk=−Hkgkd_k=-H_kg_k。下面给出了计算的函数twoloop(s,y,ρ\rho,gk)的伪代码，可知该函数内部的计算仅限于标量和向量，未出现矩阵。
函数参数s,ys,y即向量序列{sks_k}和{yky_k}，ρ\rho为元素序列{ρk\rho_k}，其元素ρk=1/sTkyk\rho_k = 1/s_k^Ty_k，gkg_k是向量，为当前的一阶导数，输出为向量z=Hkgkz=H_kg_k，即搜索方向的反方向
Functiontwoloop(s,y,ρ,gk)q=gkFori=k−1,k−2,…,k−mαi=ρisTiqq=q−αiyiHk=yTk−1sk−1/yTk−1yk−1z=HkqFori=k−m,k−m+1,…,k−1βi=ρiyTizz=z+si(αi−βi)Returnz
Function\,\,twoloop(s, y, \rho,gk)\qquad\qquad\qquad\qquad\\
q = g_k\,\, \qquad\qquad\qquad\qquad\qquad\qquad\quad\\
For \,i=k-1, k-2, \ldots, k-m\qquad\quad\\
\alpha_i = \rho_i s^{\rm T}_i q\,\!\qquad\qquad\quad\\
q = q - \alpha_i y_i\,\!\qquad\qquad\\
H_k=y^{\rm T}_{k-1} s_{k-1}/y^{\rm T}_{k-1} y_{k-1}\qquad \qquad\quad\\
z = H_k q\qquad\qquad\qquad\qquad\qquad\qquad\\
For\,\, i=k-m, k-m+1, \ldots, k-1\\
\beta_i = \rho_i y^{\rm T}_i z\,\!\qquad\qquad\\
z = z + s_i (\alpha_i - \beta_i)\,\!\\
Return \,\, z\qquad\qquad\qquad\qquad\qquad\qquad\,\!
给出L-BFGS的算法
可以看到其搜索方向dkd_k是根据函数twolooptwoloop计算得到的。
L-BFGS算法Python实现
def twoloop(s, y, rho,gk):
n = len(s) #向量序列的长度
if np.shape(s)[0] >= 1:
#h0是标量，而非矩阵
h0 = 1.0*np.dot(s[-1],y[-1])/np.dot(y[-1],y[-1])
else:
h0 = 1
a = empty((n,))
q = gk.copy()
for i in range(n - 1, -1, -1):
a[i] = rho[i] * dot(s[i], q)
q -= a[i] * y[i]
z = h0*q
for i in range(n):
b = rho[i] * dot(y[i], z)
z += s[i] * (a[i] - b)
return z
def lbfgs(fun,gfun,x0,m=5):
# fun和gfun分别是目标函数及其一阶导数,x0是初值,m为储存的序列的大小
maxk = 2000
rou = 0.55
sigma = 0.4
epsilon = 1e-5
k = 0
n = np.shape(x0)[0] #自变量的维度
s, y, rho = [], [], []
while k < maxk :
gk = gfun(x0)
if np.linalg.norm(gk) < epsilon:
break
dk = -1.0*twoloop(s, y, rho,gk)
m0=0;
mk=0
while m0 < 20: # 用Armijo搜索求步长
if fun(x0+rou**m0*dk) < fun(x0)+sigma*rou**m0*np.dot(gk,dk):
mk = m0
break
m0 += 1
x = x0 + rou**mk*dk
sk = x - x0
yk = gfun(x) - gk
if np.dot(sk,yk) > 0: #增加新的向量
rho.append(1.0/np.dot(sk,yk))
s.append(sk)
y.append(yk)
if np.shape(rho)[0] > m: #弃掉最旧向量
rho.pop(0)
s.pop(0)
y.pop(0)
k += 1
x0 = x
return x0,fun(x0),k#分别是最优点坐标，最优值，迭代次数实际性能
相关代码：
n = 50
x = y = np.linspace(-10,10,n)
xx,yy = np.meshgrid(x,y)
data = [[xx[i][j],yy[i][j]] for i in range(n) for j in range(n)]
iters = []
for i in xrange(np.shape(data)[0]):
rt = lbfgs(fun,gfun,data[i])
if rt[2] <=1000:
iters.append(rt[2])
if i%100 == 0:
print i,rt[2]
plt.hist(iters,bins=100)
plt.title(u'L-BFGS法迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})参考文献
Large-scale L-BFGS using MapReducehttp://papers.nips.cc/paper/5333-large-scale-l-bfgs-using-mapreduce.pdf
L-BFGS的原理及在回归分析中的应用
http://blog.csdn.net/hlx371240/article/details/39970727
L-BFGS的原理和MATLAB实现
deep learning Softmax分类器（L-BFGS,CG,SD）
http://blog.csdn.net/hlx371240/article/details/40015395
机器学习中导数最优化方法
http://www.cnblogs.com/daniel-D/p/3377840.html
Newton’s method in optimization
http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization
Quasi-Newton method
http://en.wikipedia.org/wiki/Quasi-Newton_method
Davidon–Fletcher–Powell formula
http://en.wikipedia.org/wiki/Davidon%E2%80%93Fletcher%E2%80%93Powell_formula
Broyden–Fletcher–Goldfarb–Shanno algorithm
http://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm
Limited-memory BFGS
http://en.wikipedia.org/wiki/Limited-memory_BFGS
some books
《Basic Concepts and Stationary Iterative Methods》
《Iterative Methods for Optimization》
《Non-Linear Least-Squares Minimization》
some papers
《A comparison of algorithms for maximum entropy parameter estimation》
《Updating Quasi-Newton Matrices with Limited Storage》
《Scalable Training of L1-Regularized Log-Linear Models》
《最优化方法及其MATLAB程序设计》，马昌凤，科学出版社
《优化方法》，李春明，东南大学出版社
《应用最优化方法及MATLAB实现》，刘兴高，胡云卿，科学出版社
《优化技术与MATLAB优化工具箱》，赵继俊，机械工业出版社
http://hub.darcs.net/amgarching/pts/browse/bfgs.py
markdown 公式，可直接复制得到特性公式
http://www.aleacubase.com/cudalab/cudalab_usage-math_formatting_on_markdown.html
CSDN-markdown语法之如何插入图片
/article/1369460.html
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- HDU 5239 上海大都会 D题（线段树+数论）
- 阅读8.9.10章
- 2013斯坦福大学iOS应用开发学习笔记 11 Table View and iPad
- 与人说话
- poj1061青蛙的约会 扩展欧几里得
- 常用正则表达式
- 1.CentOS 6.5 静态网站制作学习
- ZOJ 3713 In 7-bit
- ASP.NET MVC中如何使用PartialView
- 一步一步教你实现CTreeCtrl 自绘
- mysql之group by，order by
- IE8~实现未知节点高度实现垂直居中
- 使用Dojo（一） widget的生命周期管理
- SSDT-BI之七：循环任务（容器）
- 闲谈MongoDb+GridFS+Nginx
- Spring与Mybatis整合的MapperScannerConfigurer处理过程源码分析
- 通过代码自定义cell
- 前端常用工具整理
- 接入微信SDK64位包报错
- 某著名公司2015暑期实习招聘试题及相关内容复习