本章进入凸优化问题的求解、算法阶段。





无约束优化问题

本文讨论一下无约束问题：



 $$\min f(x)$$ 


其中，$f$是二次可微凸函数。假定该问题可解，即存在最优解$x^\star$，用$p^\star$表示最优值为：$\inf _xf(x) = f(x^\star)$。

因为$f$可微，则最优点$x^\star$满足以下条件：



 $$\triangledown f(x^\star)= 0$$ 


在特殊情况下，我们可以通过解析法求解最优性方程，但大多数情况下没有办法求得解析解。因此，最常见的方法是使用迭代法。

我们需要计算点列：$x^{(0)},x^{(1)},..$，使得$k\rightarrow \infty ,f(x^{(k)}) \rightarrow p^\star$。使用$\epsilon$表示容许误差值，当$f(x^{(k)}) - p^\star \le \epsilon$时，算法终止。



例子$_{p_{438}}$



二次优化

 $$\min (1/2)x^TPx + q^Tx +r$$ 


很容易我们可以对其求导得到$Px^\star +q= 0$。

因此，若$P$为正定矩阵，则存在唯一解$x^\star = -P^{-1}q$。

若此方程无解，则优化问题无下界。



最小二乘

作为二次优化的特例：



 $$\min ||Ax - b||^2_2=  x^T(A^TA)x - 2(A^Tb)^Tx + b^Tb$$ 


其最优性解为其正规方程：



 $$A^TAx^\star = A^Tb$$ 




强凸性

强凸是指：



 $$\triangledown^2 f(x) \succeq mI$$ 


对任意$x\in S$都成立。这里的$\triangledown^2 f(x)$表示Hessian矩阵。

我们可以通过强凸性推导出有意义的方程。



次优性条件

 $$\begin{align}f(y) &= f(x) + \triangledown f(x)^T(y-x)+\frac{1}{2}(y-x)^T \triangledown^2 f(z)(y-x) \\& \ge f(x)+ \triangledown f(x)^T(y-x)+\frac{m}{2}||y-x||_2^2  \end{align}$$ 


当$m =0$时，上式变为凸性的基本不等式（一阶可微），当$m\ge 0$时，对$f(y)$的下界得到了更好的估计结果。

对上式进行求导可得：
 $$\triangledown f(x) + m(y-x) = 0$$ 


因此，带入原式可得：



 $$\begin{align}f(y) &\ge f(x) + \triangledown f(x)^T(y-x)+\frac{m}{2}||y-x||_2^2  \\& \ge f(x) - \frac{1}{2m}||\triangledown f(x)||^2_2 \end{align}$$ 


既然该式子对任意$y\in S$都成立，则：



 $$p^\star \ge f(x) - \frac{1}{2m}||\triangledown f(x)||^2_2$$ 


由于$||\triangledown f(x)||_2 \le (2m\epsilon )^{1/2}$，因此带入可得到次优性条件：



 $$f(x) - p^\star \le \epsilon$$ 


我们也可以得到$x$与任意最优解$x^\star$之间的距离与$||\triangledown f(x)||_2$的关系：



 $$||x- x^\star||_2 \le \frac {2}{m}||\triangledown f(x)||_2$$ 




关于$\triangledown^2 f(x)$的上界

由于$||\triangledown^2 f(x)||$的最大特征值是$x$在$S$上的连续函数，因此他在$S$上有界，即存在常数$M$，使得(没懂)：



 $$||\triangledown^2 f(x)|| \preceq MI$$ 


与上面类似，我们可以得到$p^\star$的上界：



 $$p^\star \le f(x) - \frac{1}{2M}||\triangledown^2 f(x)||_2^2$$ 


对比$p^\star$的上界：



 $$p^\star \ge f(x) - \frac{1}{2m}||\triangledown f(x)||^2_2$$ 




下水平集的条件数

从之前的分析我们可以得到：



 $$mI \preceq ||\triangledown^2 f(x)|| \preceq MI$$ 


因此，比值$M/m$是矩阵$\triangledown^2 f(x)$的条件数的上界，这是影响其计算效率的重要因素。



下降方法

我们讨论的所有方法都是下降方法，只要不是最优点则应该满足：



 $$f(x^{(k+1)}) < f(x^{(k)})$$ 


而优化点列为：



 $$x^{(k+1)} = x^{(k)} + t^{(k)}\triangle x^{(k)},\quad t^{(k)} > 0$$ 


由凸性可知$\triangledown f(x^{(k)})^T(y-k^{(k)}) \ge 0$意味着$f(y)\ge f(x^{(k)})$，因此，一个下降方法的搜索方向必须满足：



 $$\triangledown f(x^{(k)})^T\Delta x^{(k)} <0$$ 


这将作为我们判断后面下降算法的条件之一。



通用下降算法


确定下降方向$\Delta x$
直线搜索。选择步长$t > 0$
修改。$x := x+t\Delta x$
直到满足停止条件




精确直线搜索$_{p_{444}}$

$t$值是通过沿着射线$\{x+t\triangle x| t\ge 0\}$优化$f$而确定的：



 $$t = \arg\min_{s\ge0} f(x+s\triangle x)$$ 


当求解式中的单变量优化问题的成本比计算搜索方向的成本低时，采用精确直线搜索。

特殊情况可以用解析的方法确定最优解。



回溯直线搜索

很多时候我们并不需要找到一个最小的$t$，只需要$f$有“足够的”减少即可。

常用的方法为：


给定下降方向$\Delta x$,参数$\alpha\in (0,0.5),\beta\in (0,1)$
设定初始$t = 1$
$if \quad f(x+t\triangle x )> f(x) + \alpha t \triangledown f(x)^T \triangle x,then \quad t :=\beta t$


回溯算法从单位步长开始，按比例逐渐减小，直到满足停止条件。

由于$\Delta x$是下降方向，$\triangledown f(x)^T  \Delta x< 0$，所以只要$t$足够小，一定有：



 $$f(x+t\Delta x) \approx f(x) + t\triangledown f(x)^T\Delta x< f(x) + t\alpha\triangledown f(x)^T\Delta x$$ 


因此回溯算法一定会停止。



梯度下降方法

梯度下降利用$\triangle x = -\triangledown f(x)$，是一种自然的选择。



算法过程


给定初始点$x\in dom f$，$\Delta x := \triangledown f(x)$
直线搜索。通过精确或回溯直线搜索确定步长$t$
修改$x:=x+t\Delta x$
直到满足停止准则




收敛性分析

我们可以推导得出：



 $$f(x^{(k)}) - p^\star \le c^k(f(x^{(0)}) - p^\star)$$ 


其中$c = 1-m/M <1$，因此最多经过$\frac{\log (f(x^{(0)}) -p^\star)/ \epsilon }{\log (1/c)}$次迭代，可以收敛到次优性条件。



例子



$R^2$空间中的二次问题

考虑二次目标函数



 $$f(x)  = \frac{1}{2}(x_1^2 + \gamma x_2^2)$$ 


其中$\gamma$大于0。很容易得到其Hessian矩阵为常数，特征值为1和$\gamma$，因此其下水平集的条件数都等于：



 $$\frac{\max(1,\gamma)}{\min(1,\gamma)} = \max (\gamma ,1/\gamma )$$ 


我们选取初始点$x^{(0)} = (\gamma ,1)$，可以计算：





最速下降方法$_{p_{454}}$

对$f(x+v)$在$x$处进行一阶Taylor展开：



 $$f(x+v) \approx \hat f(x+v) = f(x) +\triangledown f(x)^Tv$$ 


其中$\triangledown f(x)^Tv$是$f$在$x$在沿着方向$v$的方向导数，若其为负数，则$v$为下降方向。

如何选择$v$使得其方向导数尽可能小（下降最快），我们定义一个规范化的最速下降方向：



 $$\Delta x_{nsd} = \arg\min \{\triangledown f(x)^T v\space | \space ||v|| =1 \}$$ 


我们也可以考虑将最速下降方向乘以一个特殊的比例因子，从而考虑非规范化的最速下降方向：



 $$\Delta x_{sd} = ||\triangledown f(x)|| _\star \Delta x_{nsd}$$ 


其中$|| \odot ||_\star$定义为对偶范数（$||z||_\star = \sup\{z^Tx | ||x||\le1\}$）。

对于这种最速下降步径，我们有：



 $$\triangledown f(x)^T\Delta x_{sd} =||\triangledown f(x)|| _\star \triangledown f(x)^T\Delta x_{nsd} = - ||\triangledown f(x)|| _\star^2<0$$ 




Newton方法$_{p_{462}}$



Newton步径

$\Delta x_{nt} = -\triangledown^2 f(x)^{-1}\triangledown f(x)$

由$\triangledown^2f(x)$的正定性可知（凸函数的二阶条件），除非$\triangledown f(x) = 0$，否则：



 $$\triangledown f(x)^T\Delta x = -\triangledown f(x)^T \triangledown^2 f(x)^{-1}\triangledown f(x) <0$$ 


因此，$\Delta x_{nt}$与负梯度方向为锐角，$\Delta$x_{nt}是下降方向。

可以从以下几个方面来了解Newton步径。



二阶近似的最优解

考虑函数$f$在$x$处的二阶Taylor近似为：



 $$\hat{f(x+v)} = f(x)+\triangledown f(x)^Tv+\frac{1}{2}v^Tf(x)v$$ 


这是$v$的二次凸函数，在$v =\Delta x_{nt}$处达到最小值。

因此，将$x$加上Newton步径$\Delta x_{nt}$能够极小化$f$在$x$处的二阶近似。

如果函数是二次的，那么使用Newton步径是$f$的精确最优解。若函数近似二次，则$x+ \Delta x_{nt}$是$f$的最优解，即$x^\star$的很好的估计值。



Hessian范数下的最速下降方向

若定义Hessian矩阵$\triangledown^2 f(x)$定义的二次范数，即：



 $$||u||_{\triangledown^2 f(x)} = (u^T\triangledown^2 f(x)u)^{1/2}$$ 


那么可以通过最速下降方向推出$\Delta x_{nt}$。



线性化最优性条件的解

若我们在$x$附近对最优性条件$\triangledown f(x^\star)=0$进行线性化，可得到：



 $$\triangledown f(x+v) \approx \triangledown f(x) + \triangledown^2 f(x)v = 0$$ 


其解就是我们的Newton步径。



Newton减量$_{p_{464}}$

 $$\lambda(x) = (\triangledown f(x) ^T\triangledown^2 f(x)^{-1}\triangledown f(x))^{1/2}$$ 


$\lambda(x)$称为Newton减量，可以用来设计停止准则。

将$\Delta x_{nt}$带入$f(x+t\Delta x)$得：$f(x+t\Delta x) = f(x)- \frac{1}{2}\triangledown ^Tf(x)\triangledown^2 f(x)^{-1}\triangledown f(x)$

因此：$f(x) - f(x+t\Delta x)  = \frac{1}{2}\triangledown ^Tf(x)\triangledown^2 f(x)^{-1}\triangledown f(x) = \frac{1}{2}\lambda(x)^2$

因此，$\lambda^/2$是$f$在$x$处的二阶近似对$f(x) - p^\star$作出的估计。

我们也可以将Newton减量表示为$\lambda = (\Delta x_{nt}^T\triangledown ^2f(x)\Delta x_{nt})^{1/2}$，表明$\lambda$是Newton步径的二次范数，该范数由Hessian矩阵定义。

Newton减量也出现在回溯直线搜索中，因为我们可以得到：



 $$\triangledown ^Tf(x)\Delta x_{nt} = -\lambda(x)^2$$ 




Newton方法


给定初始点$x \in dom f$，误差阈值$\epsilon >0$
计算Newton步径和减量
停止准则：如果$\lambda^2/2\le \epsilon$，退出
直线搜索，根据回溯直线搜索确定步长$t$
改进：$x:=x + t\Delta x_{nt}$

9.1 无约束优化问题

例子
	
强凸性及其含义
无约束优化问题

其中是二次可微凸函数（dom(f)是开集），假设该问题可解，存在最优点，这里用表示最优值。

由于f是二次可微凸函数，最优点应满足：

所以无约束优化问题的求解变成了求解上述方程的解。一般情况下，必须采用迭代算法求解此方程，即计算点列使得时，，这样的点列被称为优化问题的极小化点列。当时，算法将终止，其中是设定的容许误差值。

初始点和下水平集

迭代的方法需要一个适当的初始点，该初始点必须满足两个条件：

初始点必须在dom(f)内
	
下水平集必须是闭集。
条件１很容易满足，条件２却不容易满足。除非在全部的下水平集都是闭集的时候，条件２才容易满足。

全部下水平集都是闭集的情况：

epi f （ｆ的上境图）是闭集
	

	
当，即x趋于dom(f)的边界时，
下水平集的条件数

对任意满足的方向向量q，我们定义凸集的宽度如下：



再定义C的最小宽度和最大宽度：



于是凸集C的条件数可以表示成：



例子

无约束几何规划



其最优性条件为：



但一般情况下，此方程组没有解析解，于是采用迭代方法，这里，所以人和店都可以是初始点。

线性不等式的解析中心



其中

采用迭代的方法，也可以找到初始点，因为当时，，满足下水平集为闭集的条件。

强凸性及其含义

假设目标函数在Ｓ上是强凸的，这是指存在使得，对任意的都成立。

对于任意的都有

由强凸性可推出





可知，对固定的x，不等式右边是y的二次凸函数，记为Ｌ，令其对y的一阶导数为０





可知Ｌ的最小值为，所以可推出



由于y是任意的，又可推出即

 

来源：https://blog.csdn.net/wangchy29/article/details/87944620

9.1无约束优化问题

例子
	
强凸性及其含义
无约束优化问题

其中是二次可微凸函数（dom(f)是开集），假设该问题可解，存在最优点，这里用表示最优值。

由于f是二次可微凸函数，最优点应满足：

所以无约束优化问题的求解变成了求解上述方程的解。一般情况下，必须采用迭代算法求解此方程，即计算点列使得时，，这样的点列被称为优化问题的极小化点列。当时，算法将终止，其中是设定的容许误差值。

初始点和下水平集

迭代的方法需要一个适当的初始点，该初始点必须满足两个条件：

初始点必须在dom(f)内
	
下水平集必须是闭集。
条件１很容易满足，条件２却不容易满足。除非在全部的下水平集都是闭集的时候，条件２才容易满足。

全部下水平集都是闭集的情况：

epi f （ｆ的上境图）是闭集
	

	
当，即x趋于dom(f)的边界时，
下水平集的条件数

对任意满足的方向向量q，我们定义凸集的宽度如下：



再定义C的最小宽度和最大宽度：



于是凸集C的条件数可以表示成：

​​​​​​​

例子

无约束几何规划



其最优性条件为：



但一般情况下，此方程组没有解析解，于是采用迭代方法，这里，所以人和店都可以是初始点。

线性不等式的解析中心



其中

采用迭代的方法，也可以找到初始点，因为当时，，满足下水平集为闭集的条件。

强凸性及其含义

假设目标函数在Ｓ上是强凸的，这是指存在使得，对任意的都成立。

对于任意的都有

由强凸性可推出





可知，对固定的x，不等式右边是y的二次凸函数，记为Ｌ，令其对y的一阶导数为０





可知Ｌ的最小值为，所以可推出



由于y是任意的，又可推出即​​​​​​​

无约束优化 
本篇文档主要介绍无约束优化问题，同时初步介绍解该类问题目前常用的一种算法即
Quasi-Newton Method (拟牛顿法)。在介绍无约束优化问题之前，我们首先会从直观上引入无约束优化的概念，并在此基础上引入解这类问题的两个重要概念：步长和方向。由步长的选择引入重要概念 line search，由方向的选择引入重要概念 Quasi-Newton Method。因此本篇介绍文档主要分为以下几个部分：无约束优化问题引入，Line Search，Quasi-Newton Method 和算法总结。 
1．无约束最优化 
 对无约束优化不熟悉的读者也许要问，什么是无约束优化。这里以一个例子来说明该问题。 
上图所示为一元函数 f(x)的图像，无约束优化问题，即不对定义域或值域做任何限制的情况下，求解函数 f(x)的小值，上面显示两个小值点：一个为全局小值点，另一个为局部小值点。受限于算法复杂度等问题，目前大部分无约束优化算法只能保证求取局部小值点。这时读者不免要问，既然只能求取局部小值点，那为什么这类算法还能应用呢。这是因为实际应用中，许多情形被抽象为函数形式后均为凸函数，对于凸函数来说局部小值点即为全局小值点，因此只要能求得这类函数的一个小值点，该点一定为全局小值点。 
理解了上面的无约束优化问题之后，我们就可以开始介绍无约束优化的求解过程
了，对于无约束优化的求解首先我们需要选择一个初始点 x_0，如下所示： 
初始点选择好之后，就可以按照各种不同的无约束优化求解算法，求解小值点了。求解过程中主要涉及两个概念，即从初始点开始沿“哪个方向”以及“走多远”到达下一个点处。所谓“走多远”即之前提的“步长”的概念，“哪个方向”即方向概念。 
2．Line Search 
Line search 主要用于解决之前提到的步长的概念，即方向确定好之后，需要确定从当前点 x_k 沿着该方向走多远，以便走到下一个合适的点 x_k+1。若用 p_k 代表从第 k 个点走向第 k+1 点的方向，x_k 代表当前点，x_k+1 代表下一个点，a_k 代表步长，则存在如下的等式： x_k+1 = x_k + a_k * p_k        (1) 
这里简要介绍一下 p_k，大部分 line search 方法要求 p_k 为下降方向，即从当前点沿着 p_k 方向移动后导致函数值减少。由于目标是求取一个函数的小值，因此优的情况是求取沿 p_k 方向满足 f(x_k+1)为全局小的 a_k 值，可用下式表示为： 
Ø (a_k) = f(x_k+1) = f(x_k + a_k * p_k)                          	(2) 
由于直接求取满足Ø (a_k)为全局小值的 a_k 涉及到大量 f(x_k + a_k * p_k)的计算，若从求导角度计算小值，还会涉及到▽f_k+1的计算，计算量较大。因此从计算量角度考虑，可以采用如下较为折中的策略。 方向确定好之后，每一步的 line Search 主要涉及两个问题：1）什么样的 a_k 是合理的 2）
如何选择 a_k 的长度。下面将沿这两个面展开讨论，首先讨论“什么样的 a_k 是合理的”，确定了该问题之后，我们就可以在此基础上选择 a_k 了。 
2.1 a_k 合理性讨论 如下将要讨论关于a_k需要满足的两个条件，当a_k满足这两个条件后，就可以认为从x_k 点移动到x_k+1 点的步长已经确定下来了。第一个条件为sufficient decrease condition，从直
观角度来看，该条件主要要用保证x_k+1 点的函数值要小于x_k点的函数值，满足该条件后，才有全局收敛[1]的可能性。第二个条件为curvature condition，从直观角度来看，该条件主要用于保证x_k点经过步长a_k的移动到达x_k+1 后，▽f_k+1小于▽f_k。 
2.1.1 	sufficient decrease condition 
a_k 的选择一定要使得函数值满足 sufficient decrease condition，该条件可以用如下不等式描述： 
f(x_k+1) ≤ f(x_k) + c_1* a_ k * ▽f_k T *  p_k 	 	 	 	 	 	(3) 
将公式（1）代入上式，可得： 
f(x_k + a_k * p_k) ≤ f(x_k) + c_1* a_ k * ▽f_k T *  p_k  	 	 	(4) 这里有必要对上面的不等式做一些解释： 
f(x_k) 代表函数在第 x_k 点的值 
▽f_k代表函数在第 x_k 点的梯度 
p_k 代表从第 x_k 点的走到 x_k+1 点的方向 
a_ k 代表从第 x_k 点沿着 p_k 方向走到 x_k+1 点步长 
c_1[2]为常量，需满足 0< c_1 < 1，一般取c_1 为 1E-4 当 p_k 为函数下降方向时，有： 
             ▽f_k * p_k  <  0                                               (5) 因此不等式 4，即要求： 
f(x_k+1) ＜  f(x_k)                                      	 	(6) 从图形角度看，函数位于第 k 点时，以上各参数中只有 a_k 为变量，其他均为常量，因
此(4)可以用以下不等式重新描述为： 
Ø(a_k) ≤ l(a_k)              (7) 
其中： 
Ø (a_k) = f(x_k + a_k * p_k)  l(a_k) = f(x_k) + c_1* a_ k * ▽f_k T *  p_k 
以下为不等式（7）的图形化表示： 
因此只要步长 a_k 的选择使得函数 Ø(a_k)位于 acceptable 区间，就满足 sufficient decrease condition。 
2.1.2 	curvature condition 
a_k 的选择一定要使得函数梯度值满足 curvature condition，该条件可以用如下不等式描
述： 
     ▽f(x_k + a_k * p_k) T *p_k ≥ c_2 * ▽f(x_k) T * p_k                         (8) 即为 
 	  ▽f _k+1 T *p_k ≥ c_2 * ▽f_k T * p_k  	 	 	 	 	 	            (9) 这里有必要对上面的不等式做一些解释： 
▽f_k+1代表函数在第 x_k+1 点的梯度 
▽f_k代表函数在第 x_k 点的梯度 
p_k 代表从第 k 点的走到 k+1 点的方向 
c_2 为常量，需满足 0< c_1 < c_2 < 1，一般取 c_2 为 0.9 当 p_k 为函数下降方向时，有： 
             ▽f_k * p_k  <  0  
因此不等式 9，即要求： 
              ▽f _k+1 ≥  c_2 * ▽f_k                                       (10) 从图形角度看，不等式10即要求函数在第x_k+1点的变化速度要低于x_k点的变化速度，这一点可以从这两点处的梯度处看出，如下图所示。 
2.1.3 Wolfe conditions 所谓 Wolfe conditions 即 sufficient decrease condition 和 curvature condition 的综合，即 a_k
需要同时满足如下两个条件： 
f(x_k + a_k * p_k) ≤ f(x_k) + c_1* a_ k * ▽f_k T *  p_k  	 	 	       (11) 
▽f(x_k + a_k * p_k) T *p_k ≥ c_2 * ▽f(x_k) T * p_k 	                       (12) 
图形化表示后，如下图所示： 
实际应用中，常常会用到由 Wolfe conditions 引申出的 strong Wolfe conditions，即 a_k 需要
同时满足如下两个条件： 
f(x_k + a_k * p_k) ≤ f(x_k) + c_1* a_ k * ▽f_k T *  p_k           (13) |▽f(x_k + a_k * p_k) T *p_k| ≤ c_2 * |▽f(x_k) T * p_k|                    (14) 与 Wolfe condtions 唯一的区别是 strong Wolfe condtions 避免了▽f(x_k + a_k * p_k)取较大的正值情况。 
2.2 a_k 步长的选择 
了解了 a_k 的合理性之后，就相当于获得了标尺，在此基础上我们可以选择合适的策略来求取 a_k。所有的 line search 过程在计算每一步的 a_k 时，均需要提供一个初始点 a_0，然后再此基础上生成一系列的{a_i}，直到 a_i 满足 2.1 节所规定的条件为止，此时该 a_k 即被确定为 a_i，或者未找到一个合适的 a_k。这里我们仅介绍目前常用的策略平方插值和立方插值法。因此本节内容分为两部分，2.2.1 节介绍选择 a_k 常用的平方插值和立方插值法，
2.2.2 节介绍由 x_k 点到 x_k+1 点，方向确定为 p_k 后，步长 a_k 具体计算过程。 
2.2.1 平方立方插值法 当给定一个初始步长 a_0，若该初始步长满足 Wolfe conditions（或 strong Wolfe 
conditions），则 a_k 被确定为 a_0，当前点 x_k 步长计算过程结束，否则，我们可以在此基础上，利用已知的三个信息 Ø(a_0)、Ø(0)、Ø’(0)，构建一个二次（平方）插值多项式来拟合Ø(a_k)。该二次插值多项式如下所示： 
Ø_q(a) =( ( Ø(a_0) - Ø(0) - a_0*Ø’(0) ) / a_0 2 )*a2 + Ø’(0)*a + Ø(0)                   (15) 观察上面的二次插值多项式可知，其满足如下插值条件： 
Ø_q(0) = Ø(0)    Ø_q’(0) = Ø’(0)     Ø_q(a_0) = Ø(a_0) 
对二次插值多项式(15)求导并令其为零，即可获得使得该多项式取得小值的 a 值，如下所示： 
a_1 = -1 * Ø’(0) * a_0 2 / ( 2 *( Ø(a_0) - Ø(0) – Ø’(0)*a_0 ) ) 	 	 	 	   (16) 若 a_1 满足 Wolfe conditions（或 strong Wolfe conditions），则 a_k 被确定为 a_1。否则在此基础上构建一个三次（立方）插值多项式，并求得使得该多项式取小值的 a 值，该 a 值的计算公式如下所示： 
a_i+ 1= a_i – (a_i – a_i-1)*(Ø’(a_i) + d_2 – d_1) / (Ø’(a_i) – Ø’(a_i-1) + 2*d_2)         (17) d_1 = Ø’(a_i-1) + Ø’(a_i) – 3*( Ø(a_i-1) - Ø(a_i) ) / ( a_i-1 – a_i) d_2 = sign(a_i - a_i-1) * ( d_12 - Ø’(a_i-1) * Ø’(a_i) ) 1/2 
若 a_i+1 满足 Wolfe conditions（或 strong Wolfe conditions），则 a_k 被确定为 a_i+1，否则一直使用三次插值多项式进行插值拟合。并且选择使用 a_i+1 对应的Ø (a_i+1)和Ø’( a_i +1) 替换 a_i-1 或 a_i 相应的值，一旦确定好 a_i-1 或 a_i 中的一种后，每次都替换对应的值，如替换 a_i-1，则每次都使用新的 a_i+1 对应的值替换 a_i-1，直到找到满足 Wolfe conditions （或 strong Wolfe conditions）的 a_i+1 为止，此时 a_k 被确定为 a_i+1，或者没有找到合适的 a_i+1。 
这里有必要简略解释一下为何只使用到三次插值多项式，而没有使用更高阶的插值多项式。原因是三次插值多项式对函数某个点处的具体值有较好的拟合效果，同时又有较好的抗过拟合作用。 
后有必要解释一下初始步长 a_0 的选择，对于 Newton 或 quasi-Newton methods 来说，初始步长 a_0 总是确定为 1，该选择确保当满足 Wolfe conditions（或 strong Wolfe conditions），我们总是选择单位 1 步长，因为该步长使得 Newton 或 quasi-Newton methods 达到较快的收敛速度。计算第零步长时，将初始步长 a_0 使用如下公式确定： 
a_0 = 1.0 / (p_0 T* p_0) 1/2   (18) 
第 1 步及其以后的初始步长 a_0，直接设定为 1。 
2.2.2 步长 a_k 计算 
本节主要做一个总结，即综合前面步长需要满足的条件及步长迭代计算公式给出步长计算的具体过程。下面假设我们处于 x_k 点，因此要从该点选择一个满足 Wolfe (或 strong 
Wolfe) conditions 的步长 a_k，以便走到下一点 x_k+1。因此我们选择初始点 a_0 = 1，Newton
或 quasi-Newton method 中初始步长总是选择为 1。因此有： 
1) 初始化 a_x  a_y  0   Ø(a_x) 、 Ø(a_y) 、Ø’(a_x) 、 Ø’(a_y)、a_min、a_max 2) 初始化 a_i  a_0、Ø’(a_i)、Ø(a_i) 
若Ø(a_i) > Ø(a_x) 
  说明步长 a_i 选择的过大，因此使得Ø(a_i)满足 Wolfe(或 strong Wolfe) conditions 的步长 a_k 应位于 a_x 和 a_i 之间，使用平方和立方插值法对 a_x 和 a_i 分别进行插值，求取两个新的步长 a_quadratic 和 a_cubic，并取两者之中和 a_x 较接近的步长，这里假设为
a_quadratic，将新的步长设为 a_qudratic，同时进行以下操作：    a_y      a_i 
Ø(a_y)   Ø(a_i) Ø’(a_y)   Ø’(a_i) a_i+1      a_quadrati 
  同时在此情况下，说明 a_k 位于 a_x 和 a_y 之间 
若Ø’(a_i)* Ø’(a_x) < 0 
   说明步长 a_i 选择的过大，因此使得Ø(a_i)满足 Wolfe(或 strong Wolfe) conditions 的步长 a_k 应位于 a_x 和 a_i 之间，使用平方和立方插值法对 a_x 和 a_i 分别进行插值，求取两个
新的步长a_quadratic和a_cubic，并取两者之中和a_i较接近的步长，这里假设为a_quadratic，
将新的步长设为 a_qudratic，同时进行以下操作： 
a_y      a_x 
Ø(a_y)   Ø(a_x) 
Ø’(a_y)   Ø’(a_x) 
a_x       a_i 
Ø(a_x)   Ø(a_i) 
Ø’(a_x)   Ø’(a_i) 
a_i+1      a_quadrati 
同时在此情况下，说明 a_k 位于 a_x 和 a_y 之间 
5) 若|Ø’(a_i)| <|Ø’(a_x) | 
   说明步长 a_i 选择的过小，因此使得Ø(a_i)满足 Wolfe(或 strong Wolfe) conditions 的步长 a_k 应位于 a_i 和 a_y 之间，使用平方和立方插值法对 a_x 和 a_i 分别进行插值，求取两个
新的步长a_quadratic和a_cubic，并取两者之中和a_i较接近的步长，这里假设为a_quadratic，
将新的步长设为 a_qudratic，同时进行以下操作：   a_x     a_i 
Ø(a_x)   Ø(a_i) 
Ø’(a_x)   Ø’(a_i) 
a_i+1      a_quadrati 
同时在此情况下，说明 a_k 位于 a_x 和 a_y 之间 
5) 若|Ø’(a_i)| ≥|Ø’(a_x) | 
说明步长a_i选择的过小，若之前已经确定出了a_k所属的范围a_x和a_y，则使得Ø( a_i ) 满足 Wolfe(或 strong Wolfe) conditions 的步长 a_k 应位于 a_i 和 a_y 之间，使用立方插值法对 a_i 和 a_y 进行插值，求取新的步长 a_cubic，新的步长设为 a_cubic；否则若 a_x 小于 a_i，说明的新的步长不够大，因此将新的步长设置为 a_max，若 a_x 大于等于 a_i，则说明新的步长不够小，将其设为 a_min，同时进行以下操作： a_x     a_i 
Ø(a_x)   Ø(a_i) 
Ø’(a_x)  Ø’(a_i) 
if 之前已经确定出了 a_k 所属的范围 a_x 和 a_y 
a_i+1    a_cubic 
else if  a_x  ＜  a_i a_i+1    a_max 
 else if  a_x  ≥  a_i         a_i+1    a_min 
6) 计算判断若 a_i+1 使得以下两式(strong Wolfe condition)均成立：    f(x_k + a_i+1 * p_k) ≤ f(x_k) + c_1* a_i+1 * ▽f_k T *  p_k 	 	 	 	    
|▽f(x_k + a_i+1 * p_k) T *p_k| ≤ c_2 * |▽f(x_k) T * p_k| 
   则找到合理的步长 a_k，将其设定为 a_i+1       a_k  a_i+1 
   则 x_k 点步长 a_k 计算结束。    否则转到 2）继续计算合理步长。 
3．Quasi-Newton Method 
    在第 2 节中我们了解了步长的概念，以及从 x_k 走到 x_k+1 点使用 line search 方法计算步长的方法。不过我们在那里忽略了一个重要的概念，即“方向”。从第 2 节，我们了解到从每一点 x_k 走到下一点 x_k+1 时，需要给出要走的“方向”，只有“方向”确定好之后，才能在此基础上应用 line search 方法找到对应的“步长”，因此在解决了“步长”计算问题之后，这里我们将和大家一起了解一下每一步的“方向”如何确定。本节分为 2 大部分，首先我们通过 newton method 引入方向的概念，在此基础上引入 quasi-newton method。然后引入 quasi-newton method 中的一种重要方法 BFGS method，并在 BFGS method 的基础上介绍用于大规模计算的 LBFGS method 算法，同时以此结束本节的所有内容。 
3.1 Newton Method 
我想我们还依稀记得，微积分中用于求一元函数根的 Newton 法。该方法可用如下所示
的图形来描述： 
首先我们选择一个初始点 x_0 并计算获得其对应的 f(x_0)，然后该方法用曲线 y = f(x)在 (x_0,f(x_0))的切线近似该曲线，把切线和 x 轴的交点记作 x_1，点 x_1 通常 x_0 更接近所要求得根，同样方法用点 x_1 处的切线近似该曲线，并求取切线和 x 轴的交点 x_2，一直迭代下去，直到找到满足我们所需的充分接近真实跟的解为止。因此 Newton 法从第 k 次近似值 x_k 求得第 k+1 次近似值 x_k+1 即为求 x_k 点处的切线和 x 轴的交点。 切线方程为： 
y – f(x_k) = ▽f (x_k)(x – x_k) 
=>  0 – f(x_k) =▽ f(x_k)(x – x_k) 
=>  ▽f(x_k)*x = ▽f(x_k)*x_k – f(x_k) =>  x = x_k – f(x_k)/▽f (x_k) 因此 x_k+1 = x_k – f(x_k)/▽f(x_k)                                           (19) 从上面使用 Newton Method 求函数根的过程可以发现，首先需要选择一个初始点，并
在点处构建一个模型来近似该函数。 
切线模型： M_k(x_k+p_k) = f(x_k) + ▽f (x_k) T *p_k 
上面使用了相应点处的切线模型来近似函数，然后求取该近似模型的根以便求得更接近函数根的下一个点，该过程一直迭代下去，直到找到根为止。从上面构建每个点处的近似模型可以发现，该模型相对原函数来说简化了很多，因此求解要容易一些。 
现在来考虑求取函数的小值问题，方法类似。首先开始我们选择一个初始点 x_0，并构建函数在该点处的一个近似模型，上面求函数根时，我们构建的近似模型为切线模型。这里我们构建一个抛物线模型： 
抛物线模型：M_k(x_k + p_k) = f(x_k) + ▽f(x_k) T *p_k + 1/2 * p_k T *▽2f(x_k) *p_k 
并求解该模型的梯度，同时令其为零，即：▽M_k(x+) = 0，在此基础上求得 x+，该值
即使得近似模型取得小值的点。 
对抛物线模型 M_k 求导并令其为零后，可得以下公式： 
p_k = –▽2f(x_k) -1*▽f(x_k)                    (20) 因此 x_k+1 = x_k –▽2f(x_k) -1*▽f(x_k)               (21) 
从上式可见，使用Newton法时，每一步的方向为p_k，而步长为 1。从Newton法方向公式我们不难看出，若使用Newton法，则从每个小值近似点x_k走到下一个近似点x_k+1 的过程中将涉及函数Hessian矩阵▽▽f(x_k) (二阶导)计算，而该Hessian矩阵无法保证在每个点处都是正定[3]的，对正定矩阵来说存在如下不等式： 
p_k T *▽2f(x_k) *p_k  >  0 	 	 	 	 	 	 	 	 	 	   (22) 由于矩阵无法保证为正定矩阵，因此下式中 
–▽f(x_k) T *p_k = p_k T *▽2f(x_k) *p_k 	 	 	 	 	 	 	 	   (23) p_k 无法保证总是为下降方向，即负梯度方向。 从上面的分析可见，Newton Method 面临两个问题： 
Hessian 矩阵计算量较大 
Hessian 矩阵无法保证在迭代的过程中始终处于正定状态 
为了解决这两个问题，我们引入 Quasi-Newton Method。 
3.2 Quasi-Newton Method 
   Quasi-Newton Method 每一步计算过程中仅涉及到函数值和函数梯度值计算，这样有效避免了 Newton Method 中涉及到的 Hessian 矩阵计算问题。于 Newton Method 不同的是 Quasi-Newton Method 在每点处构建一个如下的近似模型： 
M_k(x_k + p_k) = f(x_k) + ▽f(x_k) T *p_k + 1/2 * p_k T *B_k*p_k 
从上面的近似模型我们可以看出，该模型用 B_k 代替了 Newton Method 中近似模型中涉及到的 Hessian 矩阵。因此 Quasi-Newton Method 中方向计算公式如下所示： 
p_k = –B_k -1*▽f(x_k)  (24) 
这里有必要解释一下用于近似 Hessian 矩阵的 B_k 可行性，及一个指导性方案。根据
Taylor(泰勒)级数可知如下公式： 
▽f(x_k + p_k) = ▽f(x_k) +▽2f(x_k) T *p_k + ∫[▽2f(x_k+t*p_k)-▽2f(x_k)]p_kdt 
由于函数▽f(.)连续，因此上式可以表示为： 
▽f(x_k+1) = ▽f(x_k) + ▽2f(x_k)*(x_k+1 – x_k) + o(||x_k+1 – x_k||) 
▽2f(x_k)*(x_k+1 – x_k) ≈ ▽f(x_k+1) - ▽f(x_k) 	 	 	 	 	 	   (25) 
因此每一选择Hessian矩阵的近似B_ k+1时，可以像式(24)那样模仿真实的Hessian
矩阵的性质。得到下式： 
B_k+1*s_k = y_k  	 	 	 	 	 	 	 	 	   (26) 
其中： 
     s_k = x_k+1 – x_k     y_k = ▽f(x_k+1) - ▽f(x_k)        (27） 
同时要求 B_k+1 为对称正定矩阵。 
3.2.1 	BFGS Method 
   从 Quasi-Newton Method 方向公式 (24) 中，可以看到每一步计算方向的过程中均涉及到 B_k+1 矩阵求逆的问题，为了避免该计算，通过分析公式（26）可知，我们可以构建一个近似 H_k+1，该近视满足如下方程： 
H_k+1*y_k = s_k 	 	         	 	 	   	 	 	   (28) 
同时要求 H_k+1 为对称正定矩阵。因此 BFGS Method 中，每个点处的方向由如下公式计算：                   p_k = –H_k*▽f(x_k)  	 	 	 	 	 	 	 	       (29) 
在此基础上，BFGS 方向迭代公式如下所示： 
H_k+1 = (1 – ρ_k *s_k*y_k T) H_k (1-ρ_k * y_k *s_k T) +ρ_k * s_k* s_k T             (30) 
其中ρ_k 为一个标量： 
          ρ_k  = 1 / (y_k T * s_k) 
有了上面(30)的 H_k 迭代公式后，还有一个问题就是初始的 H_0 如何计算，目前常用的方法是初始的 H_0 直接设为单位矩阵 I。因此 BFGS Method 用于解无约束优化的过程可以表示为如下过程： 
选择一个初始点 x_0，并选择收敛判断条件 ε> 0，及初始 H_0 (单位矩阵 I) 
k  0 
while ||▽f(x_k)||  > ε 
计算从当前点 x_k 走到下一个点 x_k+1 的方向 
 p_k = –H_k*▽f(x_k) 
采用 line search 策略计算步长 a_k 
x_k+1 = x_k + a_k * p_k 
s_k = x_k+1 – x_k     y_k = ▽f(x_k+1) - ▽f(x_k) 
H_k+1 = (1 – ρ_k *s_k*y_k T) H_k (1-ρ_k * y_k *s_k T) +ρ_k * s_k* s_k T 
k  k+1 
3.2.2 	LBFGS Method 
   上一节所介绍的 BFGS Method 比较适合解决中小规模无约束优化问题，但是 BFGS 算法产生的 Hessian 近似矩阵 H_k 为 n * n 的，同时该矩阵非稀疏，因此当 n 的规模较大时将面临两个问题： 
1） 存储问题：n 规模较大时，n*n 矩阵对内存的消耗将较大； 2） 计算问题：n 规模较大，同时 n*n 矩阵非稀疏时，计算复杂度将较高； 
为了解决以上问题，引申出了 Limited-Memory Quasi-Newton Method，目前使用较多的
LBFGS 算法即属于该类算法。为了减少 H_k 矩阵的存储，LBFGS 算法利用近几代的
curvature 信息来构建 Hessian 矩阵的近似。由 BFGS Method 我们知道： 
x_k+1 = x_k + a_k * H_k*▽f(x_k) 
其中 a_k 为步长，H_k 为 Hessian 矩阵的近似，可以通过如下迭代公式计算： 
H_k+1 = V_k* H_k*V_k+ρ_k * s_k* s_k          (31) 
其中： 
           ρ_k  = 1 / (y_k T * s_k)   V_k = 1-ρ_k * y_k *s_k T            s_k = x_k+1 – x_k     y_k = ▽f(x_k+1) - ▽f(x_k) 
从上面的 H_k 的迭代计算公式可知，H_k 会慢慢由稀疏矩阵转变为稠密矩阵，因此存储该矩阵以及进行该矩阵和向量的相乘运算的消耗将较大。为了避免该问题，LBFGS 算法在 BFGS 算法的基础上从两点进行了改进： 
1）估算每一步对应的 Hessian 近似矩阵时，给出一个当前步的初始 Hessian 矩阵估计 H_k0 
2) 利用过去当前代及过去 m-1 代的 curvature 信息修正初始 Hessian 矩阵估计 H_k0，得到终的 Hessian 矩阵近似估计 H_k。 计算式如下所示： 
H_k = (V_k-1 T … V_k-m T) *H_k0*(V_k-m … V_k-1) 
     +ρ_k-m * (V_k-1 T … V_k-m+1 T)* s_k-m*s_k-m T *(V_k-m+1 … V_k-1) 
     +ρ_k-m+1*(V_k-1 T … V_k-m+2 T)*s_k-m+1*s_k-m+1 T *(V_k-m+2 … V_k-1) 
     + … 
     +ρ_k-1*s_k-1*s_k-1 T  (32) 
上述计算式(32），可以通过公式(31)递归计算获取。公式(32)可以用以下算法表示： 
q   ▽f(x_k) 
for i = k-1, k-2, …, k-m 
α_i  ρ_i * s_i T * q; 
q     q -α_i *y_i 
end 
r   H_k0 * q; 
for i = k-m, k-m+1, …, k-1 
β ρ_i * y_i T * r 
r  r + s_i*(α_i –β) 
end 
stop with result H_k*▽f(x_k) = r 
从上面计算 H_k 的公式(32)可知，要估算每个点 x_k 处的 Hessian 矩阵近似，需要给出
初始估计 H_k0，H_k0 一般通过以下公式计算： 
H_k0 = r_k*I r_k  =( s_k-1T*y_k-1)(/y_k-1 T *y_k-1) 
有了上面的方向计算算法后，LBFGS 算法用于解无约束优化问题，可以表示为如下算法： 
选择一个初始点 x_0，并选择收敛判断条件 ε> 0，以及常量 m(代表过去代数)一般为 6 
k  0  H_0  I，因此 r = H_0 *▽f(x_0) =▽f(x_0) 
while ||▽f(x_k)||  > ε 
计算从当前点 x_k 走到下一个点 x_k+1 的方向 
 p_k = –r 
采用 line search 策略计算步长 a_k 
x_k+1 = x_k + a_k * p_k 
if  k > m 
        删除 LBFGS 计算 H_k 时用不上的向量对(s_k-m, y_k-m) 
计算并保存 s_k = x_k+1 – x_k     y_k = ▽f(x_k+1) - ▽f(x_k) 
采用 LBFGS Hessian 矩阵近似算法计算 r 
k  k+1 
4．算法总结 用于解无约束优化算法的 Quasi-Newton Method 中的 LBFGS 算法到这里总算初步介绍完了，不过这里笔者要承认的是这篇文档省略了许多内容，包括算法收敛性的证明以及收敛速度证明等许多内容。因此读者若希望对这一块有一个更深入的认识可以参考以下两本书： 
Numerical Methods for Unconstrained Optimization and Nonlinear Equations 
J.E. Dennis Jr.   Robert B. Schnabel 
Numerical Optimization 
Jorge Nocedal   Stephen J. Wright 
同时我想引用一下侯捷老师的话： 
大哉问。学习需要明师。但明师可遇不可求，所以退而求其次你需要好书，并尽早建立自修的基础。迷时师渡，悟了自渡，寻好书看好书，就是你的自渡法门。切记，徒学不足以自行，计算机是实作性很强的一门科技，你一定要动手做，忌讳眼高手低。学而不思则罔，思而不学则殆，一定要思考、沉淀、整理。 因此这里附上一个我阅读后感觉代码还比较美观的 LBFGS 实现（C++）： 
http://www.chokkan.org/software/liblbfgs/ 
Naoaki Okazaki 
后我不得不承认这篇文档，对许多读者来说也许是几张废纸，甚至使读者昏昏欲睡，
或烦躁。因为该文档从头至尾并未涉及到一个例子。-_- || 
老天啊，饶恕我吧，毕竟愿望是美好的，但愿它能对你理解 LBFGS 算法有一点帮助，
哪怕是一点点… 
 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	       jianzhu  	 	 	 	 	 	 	 	 	 	 	 	 	 	 	jzhu.nlp@gmail.com 
                                                            2009-11-1-2009-12-1 
 注：这里以及后面所说的全局收敛，均指的是可以收敛到一个局部小值处  ↑
 注：Quasi-Newton Method 中要求 0< c_1< 0.5 ↑
 关于矩阵正定可以查阅 Linear Algebra and Its Applications 一书的后一章  ↑