如果学习机器学习算法，你会发现，其实机器学习的过程大概就是定义一个模型的目标函数 J(θ) <script type="math/tex" id="MathJax-Element-1">J(\theta)</script>，然后通过优化算法从数据中求取 J(θ) <script type="math/tex" id="MathJax-Element-2">J(\theta)</script>取得极值时对应模型参数 θ <script type="math/tex" id="MathJax-Element-3">\theta</script>的过程，而学习到的参数就对应于机器学习到的知识。不管学习到的是好的还是无用的，我们知道这其中的动力引擎就是优化算法。在很多开源软件包中都有自己实现的一套优化算法包，比如stanford-nlp，希望通过本篇简要介绍之后，对于开源软件包里面的优化方法不至于太陌生。本文主要介绍三种方法，分别是梯度下降，共轭梯度法(Conjugate Gradient Method)和近似牛顿法(Quasi-Newton)。具体在stanford-nlp中都有对应的实现，由于前两种方法都涉及到梯度的概念，我们首先从介绍梯度开始。
梯度(Gradient)
什么是梯度，记忆中好像和高数里面的微积分有关。好，只要您也有这么一点印象就好办了，我们知道微积分的鼻祖是牛顿，人家是经典力学的奠基人，那么我们先来看看一道简单的物理问题：
一个小球在一个平面运动，沿着x轴的位移随时间的变化为： Sx=20−t2 <script type="math/tex" id="MathJax-Element-4">S_x=20-t^2</script>，沿着y轴的位移随时间的变化为： Sy=10+2t2 <script type="math/tex" id="MathJax-Element-5">S_y=10+2t^2</script>，现在求在 t0 <script type="math/tex" id="MathJax-Element-6">t_0</script>时刻小球的速度 v <script type="math/tex" id="MathJax-Element-7">v</script>？
大家都是为高考奋战过的人，这样的小题应该是送分题吧。牛老师告诉我们，只要通过求各个方向的分速度，然后再合成就可以求解得出。好，现在我们知道各个方向的位移关于时间的变化规律，我们来求各个方向的速度。如何求速度呢，牛老师说位移的变化率就是速度，那么我们来求在t0<script type="math/tex" id="MathJax-Element-8">t_0</script>时刻的变化率：
vx(t0)=limΔt−>0(20−(t0+Δt)2)−(20−t20)Δt=−2t0
<script type="math/tex; mode=display" id="MathJax-Element-9">v_x(t_0)=lim_{\Delta t->0}\frac{(20-(t_0+\Delta t)^2)-(20-t_0^2)}{\Delta t}=-2t_0</script>
vy(t0)=limΔt−>0(10+2(t0+Δt)2)−(10+2t20)Δt=4t0
<script type="math/tex; mode=display" id="MathJax-Element-10">v_y(t_0)=lim_{\Delta t->0}\frac{(10+2(t_0+\Delta t)^2)-(10+2t_0^2)}{\Delta t}=4t_0</script>，那么此时的合成速度
v <script type="math/tex" id="MathJax-Element-11">v</script>:
v=vx+vy=[−2t0,4t0]
<script type="math/tex; mode=display" id="MathJax-Element-12">v=v_x+v_y=[-2t_0,4t_0]</script>，此时的速度方向就是总位移变化最大的方向。搬到数学中，对于一个位移函数
S(x,y) <script type="math/tex" id="MathJax-Element-13">S(x,y)</script>，它各个维度的变化率就是其对应的偏导数
∂S(x,y)∂x,∂S(x,y)∂y
<script type="math/tex; mode=display" id="MathJax-Element-14">\frac{\partial S(x,y)}{\partial x},\frac{\partial S(x,y)}{\partial y}</script>，两者组合起来的向量就是该函数的梯度，所代表的含义上面已经说过，其方向代表函数变化最大的方向，模为变化率的大小。如果我们分别沿着
x,y <script type="math/tex" id="MathJax-Element-15">x,y</script>两个维度做微小的变化
Δx,Δy <script type="math/tex" id="MathJax-Element-16">\Delta x,\Delta y</script>,那么位移总体的变化将如下：
ΔS(x,y)≈∂S(x,y)∂x∗Δx+∂S(x,y)∂y∗Δy
<script type="math/tex; mode=display" id="MathJax-Element-17">\Delta S(x,y)\approx\frac{\partial S(x,y)}{\partial x}* \Delta x + \frac{\partial S(x,y)}{\partial y} * \Delta y</script>现在我们知道如何求取函数的梯度，而且如何利用梯度求取函数微小变化量了。
梯度下降法
做机器学习(监督学习)的时候，一般情况是这样的，有 N <script type="math/tex" id="MathJax-Element-18">N</script>条训练数据(X(i),y(i))<script type="math/tex" id="MathJax-Element-19">(X^{(i)},y^{(i)})</script>，我们的模型会根据 X <script type="math/tex" id="MathJax-Element-20">X</script>预测出对应的y<script type="math/tex" id="MathJax-Element-21">y</script>，也就是：
y(i)预测=f(X(i);θ)
<script type="math/tex; mode=display" id="MathJax-Element-22">y^{(i)}_{预测}=f(X^{(i)};\theta)</script>其中
θ=[θ1,θ2,θ3,...,θn] <script type="math/tex" id="MathJax-Element-23">\theta = [\theta _1,\theta_2,\theta_3,...,\theta_n]</script>是模型的参数。通常我们希望预测值和真实值是一致的，所以会引出一个惩罚函数:
cost(y(i)预测,y(i)真实)
<script type="math/tex; mode=display" id="MathJax-Element-24">cost(y^{(i)}_{预测},y^{(i)}_{真实})</script>而目标函数则是：
J(θ)=∑i=0Ncost(y(i)真实,y(i)预测)=∑i=0Ncost(y(i)真实,f(X(i);θ))
<script type="math/tex; mode=display" id="MathJax-Element-25">J(\theta)=\sum^N_{i=0}cost(y^{(i)}_{真实},y^{(i)}_{预测})=\sum^N_{i=0}cost(y^{(i)}_{真实},f(X^{(i)};\theta))</script>，我们目的是解决下面的优化问题:
argminθJ(θ)
<script type="math/tex; mode=display" id="MathJax-Element-26">argmin_{\theta} J(\theta)</script>，一般一组
θ <script type="math/tex" id="MathJax-Element-27">\theta</script>和
N <script type="math/tex" id="MathJax-Element-28">N</script>条数据会对应一个
J(θ)<script type="math/tex" id="MathJax-Element-29">J(\theta)</script>，也就是
N <script type="math/tex" id="MathJax-Element-30">N</script>维平面上的一个点，那么不同
θ<script type="math/tex" id="MathJax-Element-31">\theta</script>就可以得到一个
N <script type="math/tex" id="MathJax-Element-32">N</script>维的超平面(hyper plane)。特殊的假如
N=3<script type="math/tex" id="MathJax-Element-33">N=3</script>，我们可能的超平面就如下图所示：
如何找到最优的
θ <script type="math/tex" id="MathJax-Element-34">\theta</script>呢？一个想法是这样的：我们随机在超平面上取一个点，对应我们
θ <script type="math/tex" id="MathJax-Element-35">\theta</script>的初始值，然后每次改变一点
Δθ <script type="math/tex" id="MathJax-Element-36">\Delta \theta</script>，使
J(θ) <script type="math/tex" id="MathJax-Element-37">J(\theta)</script>也改变
ΔJ(θ) <script type="math/tex" id="MathJax-Element-38">\Delta J(\theta)</script>，只要能保证
ΔJ<0 <script type="math/tex" id="MathJax-Element-39">\Delta J < 0</script>就一直更新
θ <script type="math/tex" id="MathJax-Element-40">\theta</script>直到
J(θ) <script type="math/tex" id="MathJax-Element-41">J(\theta)</script>不再减少为止。具体如下：
- 随机初始化 θ <script type="math/tex" id="MathJax-Element-42">\theta</script>
- 对于每一个 θi <script type="math/tex" id="MathJax-Element-43">\theta_i</script>选择合适的 Δθi <script type="math/tex" id="MathJax-Element-44">\Delta \theta_i</script>，使得 J(θ+Δθ)−J(θ)<0 <script type="math/tex" id="MathJax-Element-45">J(\theta+\Delta \theta)-J(\theta)<0</script>,如果找不到这样的 Δθ <script type="math/tex" id="MathJax-Element-46">\Delta \theta</script>，则结束算法
- 对于每一个 θi <script type="math/tex" id="MathJax-Element-47">\theta_i</script>进行更新： θi=θi+Δθi <script type="math/tex" id="MathJax-Element-48">\theta_i = \theta_i + \Delta \theta_i</script>，回到第2步。
想法挺好的，那么如何找到所谓合适的 Δθ <script type="math/tex" id="MathJax-Element-49">\Delta \theta</script>呢？根据上一节中我们知道：
ΔJ(θ)≈∑i∂J(θ)∂θi∗Δθi
<script type="math/tex; mode=display" id="MathJax-Element-50">\Delta J(\theta) \approx \sum_i \frac{\partial J(\theta)}{\partial \theta_i}*\Delta \theta_i</script>,要如何保证
ΔJ(θ)<0 <script type="math/tex" id="MathJax-Element-51">\Delta J(\theta) <0</script>呢？我们知道，两个非0数相乘，要保证大于0，只要两个数一样即可，如果我们要保证
ΔJ(θ)>0 <script type="math/tex" id="MathJax-Element-52">\Delta J(\theta)>0</script>，只要另每一个
θi＝∂J(θ)∂θi <script type="math/tex" id="MathJax-Element-53">\theta_i＝\frac{\partial J(\theta)}{\partial \theta_i}</script>即可，此时
ΔJ(θ)≈∑i∂J(θ)∂θi∗∂J(θ)∂θi>0
<script type="math/tex; mode=display" id="MathJax-Element-54">\Delta J(\theta) \approx \sum_i \frac{\partial J(\theta)}{\partial \theta_i}*\frac{\partial J(\theta)}{\partial \theta_i} > 0</script>，有人疑问了，我们目标可是要使
ΔJ(θ)<0 <script type="math/tex" id="MathJax-Element-55">\Delta J(\theta)<0</script>，上面的做法刚好相反啊！反应快的人可能马上想到了，我们只需另每一个
θi＝−∂J(θ)∂θi <script type="math/tex" id="MathJax-Element-56">\theta_i＝-\frac{\partial J(\theta)}{\partial \theta_i}</script>不就好了，而这样的求取向量
θ <script type="math/tex" id="MathJax-Element-57">\theta</script>各个维度相对于
J(θ) <script type="math/tex" id="MathJax-Element-58">J(\theta)</script>的偏导数实际上就是求取
J(θ) <script type="math/tex" id="MathJax-Element-59">J(\theta)</script>的梯度！回忆上一节梯度的含义，表示
J(θ) <script type="math/tex" id="MathJax-Element-60">J(\theta)</script>变化最大的方向，想象一个球在上面的图中上方滚下来，而我们的做法是使他沿着最陡的方向滚。不错，我们找到了上述算法所说的合适的
Δθ <script type="math/tex" id="MathJax-Element-61">\Delta \theta</script>了！其实上述的算法就是我们本文的主角——梯度下降法(gradient descent)，完整算法如下：
- 随机初始化 θ <script type="math/tex" id="MathJax-Element-62">\theta</script>
- 求取 θ <script type="math/tex" id="MathJax-Element-63">\theta</script>的梯度 ∇θ <script type="math/tex" id="MathJax-Element-64">\nabla \theta</script>，也就是对于每个 θi <script type="math/tex" id="MathJax-Element-65">\theta_i</script>求取其偏导数 ∂J(θ)∂θi <script type="math/tex" id="MathJax-Element-66">\frac{\partial J(\theta)}{\partial \theta_i}</script>，并更新 θi=θi−η∗∂J(θ)∂θi(η>0并足够小) <script type="math/tex" id="MathJax-Element-67">\theta_i = \theta_i - \eta * \frac{\partial J(\theta)}{\partial \theta_i}(\eta > 0并足够小)</script>
- 判断 ∇θ <script type="math/tex" id="MathJax-Element-68">\nabla \theta</script>是否为0或者足够小，是就输出此时的 θ <script type="math/tex" id="MathJax-Element-69">\theta</script>，否则返回第2步
上述算法的第二步中多了一个未曾介绍的 η <script type="math/tex" id="MathJax-Element-70">\eta</script>，这是步伐大小，因为求取每一个维度的偏导，只是求取了该维度上的变化率，具体要变化多大就由 η <script type="math/tex" id="MathJax-Element-71">\eta</script>控制了， η <script type="math/tex" id="MathJax-Element-72">\eta</script>的选取更多考验的是你的工程能力，取太大是不可行的，这样导致算法无法收敛，取太小则会导致训练时间太长，有兴趣的可参考An overview of gradient descent optimization algorithms这篇文章中对 η <script type="math/tex" id="MathJax-Element-73">\eta</script>选取的一些算法。如何计算 ∂J(θ)∂θi <script type="math/tex" id="MathJax-Element-74">\frac{\partial J(\theta)}{\partial \theta_i}</script>呢？根据定义，可如下计算：
∂J(θ)∂θi=∑i=0N∂cost(y(i)真实,f(X(i);θ))∂θi
<script type="math/tex; mode=display" id="MathJax-Element-75">\frac{\partial J(\theta)}{\partial \theta_i}=\sum^N_{i=0}\frac{\partial cost(y^{(i)}_{真实},f(X^{(i)};\theta))}{\partial \theta_i}</script>，由于每次计算梯度都需要用到所有
N <script type="math/tex" id="MathJax-Element-76">N</script>条训练数据，所以这种算法也叫批量梯度下降法(Batch gradient descent)。在实际情况中，有时候我们的训练数据数以亿计，那么这样的批量计算消耗太大了，所以我们可以近似计算梯度，也就是只取
M(M<<N)<script type="math/tex" id="MathJax-Element-77">M(M<
共轭梯度法
上一节中，我们介绍了一般的梯度下降法，这是很多开源软件包里面都会提供的一种算法。现在我们来看看另外一种软件包也经常见到算法——共轭梯度法(Conjugate Gradient Method)，Jonathan在94年的时候写过《An Introduction to the Conjugate Gradient Method Without the Agonizing Pain》详细而直观地介绍了CG，确实文如其名。这里我只是简要介绍CG到底是一个什么样的东西，具体还需阅读原文，强烈推荐啊！
最陡下降法(Steepest Descent)
上一节中，我们介绍了如何反复利用 θi=θi−η∗∂J(θ)∂θi <script type="math/tex" id="MathJax-Element-78">\theta_i = \theta_i - \eta * \frac{\partial J(\theta)}{\partial \theta_i}</script>求得最优的 θ <script type="math/tex" id="MathJax-Element-79">\theta</script>，但是我们说选取 η <script type="math/tex" id="MathJax-Element-80">\eta</script>是一个艺术活，这里介绍一种 η <script type="math/tex" id="MathJax-Element-81">\eta</script>的选取方式。首先明确一点，我们希望每次改变 θ <script type="math/tex" id="MathJax-Element-82">\theta</script>，使得 J(θ) <script type="math/tex" id="MathJax-Element-83">J(\theta)</script>越来越小。在梯度确定的情况下，其实 ΔJ(θ) <script type="math/tex" id="MathJax-Element-84">\Delta J(\theta)</script>是关于 η <script type="math/tex" id="MathJax-Element-85">\eta</script>的一个函数：
ϕ(η)=J(θ−η∗∇θ)−J(θ)=ΔJ(θ)
<script type="math/tex; mode=display" id="MathJax-Element-86">\phi(\eta)=J(\theta-\eta*\nabla \theta) - J(\theta)=\Delta J(\theta)</script>，既然我们想让
J(θ) <script type="math/tex" id="MathJax-Element-87">J(\theta)</script>减小，那么干脆每一步都使得
|ΔJ(θ)| <script type="math/tex" id="MathJax-Element-88">|\Delta J(\theta)|</script>最大好了，理论上我们可以通过求导求极值，令：
ϕ′(η)=0
<script type="math/tex; mode=display" id="MathJax-Element-89">\phi'(\eta)=0</script>求得此时的
η <script type="math/tex" id="MathJax-Element-90">\eta</script>，这样每次改变
J(θ) <script type="math/tex" id="MathJax-Element-91">J(\theta)</script>是最大的，而实际操作中，我们一般采用line search的技术来求取
η <script type="math/tex" id="MathJax-Element-92">\eta</script>，也就是固定此时的梯度
∇θ <script type="math/tex" id="MathJax-Element-93">\nabla \theta</script>，也就是固定方向，尝试不同的
η <script type="math/tex" id="MathJax-Element-94">\eta</script>值，使得
J(θ−η∗∇θ)
<script type="math/tex; mode=display" id="MathJax-Element-95">J(\theta-\eta*\nabla \theta)</script>近似最小即可。这样固定方向的搜索和直线搜索没太大区别，也是名字的由来。表面看来，最陡梯度下降法应该是最优的啊，不但方向上是最陡的，而且走的步伐也是”最优”的，但是实际应用该方法的地方并不多，因为步伐的局部最优并不代表全局最优，导致其实际表现收敛速度比较慢(这却是优化算法的致命弱点！)。如果
J(θ) <script type="math/tex" id="MathJax-Element-96">J(\theta)</script>是一个二次函数也就是
J(θ)=θTAθ+bTθ+c <script type="math/tex" id="MathJax-Element-97">J(\theta)=\theta^TA\theta+b^T\theta+c</script>，通过运行算法，我们可以得到一个如下的轨迹：
我们可以发现，每一次走的步伐和上一次都是垂直的(事实上是可以证明的，在前面我推荐的文中有详细的证明:-))，这样必然有很多步伐是平行的，造成同一个方向要走好几次。研究最优化的人野心就来了，既然同一个方向要走好几次，能不能有什么办法，使得同一个方向只走一次就可以了呢？
共轭梯度法
Cornelius，Magnus和Eduard经过研究之后，便设计了这样的方法——共轭梯度法。
具体详细的原理还是强烈推荐《An Introduction to the Conjugate Gradient Method Without the Agonizing Pain》一文，这里我只是利用文章中的思路进行简要介绍。
何谓共轭(Conjugate)
查看维基的介绍，共轭梯度法(CG)最早的提出是为了解决大规模线性方程求解，比如下面形式:
Ax=b
<script type="math/tex; mode=display" id="MathJax-Element-98">Ax=b</script>，其中
A <script type="math/tex" id="MathJax-Element-99">A</script>是方形对称半正定的。如果
A<script type="math/tex" id="MathJax-Element-100">A</script>越大并且越稀疏，导致其他线性方程求解不可行的时候，CG方法就更奏效。
我们已经了解梯度为何物，现在就差修饰词共轭(Conjugate)了，那么何为共轭(conjugate)呢？对于两个非零向量
d(i),d(j) <script type="math/tex" id="MathJax-Element-101">d_{(i)},d_{(j)}</script>,如果满足
dT(i)Ad(j)=0
<script type="math/tex; mode=display" id="MathJax-Element-102">d_{(i)}^TAd_{(j)}=0</script>我们就称
d(i),d(j) <script type="math/tex" id="MathJax-Element-103">d_{(i)},d_{(j)}</script>是关于
A <script type="math/tex" id="MathJax-Element-104">A</script>共轭的，也就是说其实共轭不共轭是相对于
A<script type="math/tex" id="MathJax-Element-105">A</script>而言。我们知道，如果两个向量正交，他们的内积为0，所以如果两个向量关于
A <script type="math/tex" id="MathJax-Element-106">A</script>是共轭的，我们也可以称这两个向量关于
A<script type="math/tex" id="MathJax-Element-107">A</script>是正交的，也就是
A−orthogonal <script type="math/tex" id="MathJax-Element-108">A-orthogonal</script>。
共轭梯度法求解线性方程组
那么求解上面线性方程组的时候，假如我们已经找到 n <script type="math/tex" id="MathJax-Element-829">n</script>个两两共轭的方向{d(i)}<script type="math/tex" id="MathJax-Element-830">\{d_{(i)}\}</script>，如果将这些方向作为基，也就可以将 Ax=b <script type="math/tex" id="MathJax-Element-831">Ax=b</script>的解表示为 d(i) <script type="math/tex" id="MathJax-Element-832">d_{(i)}</script>的线性组合：
x∗=∑inαid(i)
<script type="math/tex; mode=display" id="MathJax-Element-833">x_*=\sum_i^n \alpha_i d_{(i)}</script>如果左右分别乘上
A <script type="math/tex" id="MathJax-Element-834">A</script>:
Ax∗=∑inαiAd(i)
<script type="math/tex; mode=display" id="MathJax-Element-835">Ax_*=\sum_i^n \alpha_i A d_{(i)}</script>
b=∑inαiAd(i)
<script type="math/tex; mode=display" id="MathJax-Element-836">b=\sum_i^n \alpha_i A d_{(i)}</script>
dT(k)b=∑inαidT(k)Ad(i)=αkdT(k)Ad(k)
<script type="math/tex; mode=display" id="MathJax-Element-837">d_{(k)}^Tb=\sum_i^n \alpha_i d^T_{(k)} A d_{(i)}\\=\alpha_k d^T_{(k)} A d_{(k)}</script>，也就是
αk=dT(k)bdT(k)Ad(k) <script type="math/tex" id="MathJax-Element-838">\alpha_k=\frac{d_{(k)}^T b}{d^T_{(k)} A d_{(k)}}</script>。现在问题就只在于如何找到这些神奇的共轭方向了，神奇之处在于这些共轭方向可以利用迭代方式求取，可以一次只求一个这样的方向向量
d(k) <script type="math/tex" id="MathJax-Element-839">d_{(k)}</script>，这也是该算法的核心之处。如果令
rk=b−Axk <script type="math/tex" id="MathJax-Element-840">r_k=b-Ax_k</script>,那么
d(k+1)=rk+βkd(k)
<script type="math/tex; mode=display" id="MathJax-Element-841">d_{(k+1)}=r_k+\beta_k d_{(k)}</script>其中
βk=rTk+1rk+1rTkrk <script type="math/tex" id="MathJax-Element-842">\beta_k = \frac{r^T_{k+1} r_{k+1}}{r_k^Tr_k}</script>
上一小节中利用Steepest Method的优化问题如果利用CG就变成了如下：
二维的情况下，可以保证只走两步就达到收敛(严格的证明请参靠推荐的论文)！
非线性共轭梯度法
机器学习算法中，我们碰到的大部分问题都是非线性的，上面我们只是讲解了线性共轭梯度法，那么它可以解决非线性优化问题吗？很遗憾，不行，但是经过修改，可以利用共轭梯度法求取局部最优解，下面展示非线性共轭梯度法的大致轮廓,对于一个非线性目标函数 J(θ) <script type="math/tex" id="MathJax-Element-122">J(\theta)</script>
- 随机初始化 θ0 <script type="math/tex" id="MathJax-Element-123">\theta_0</script>，并令 r0=d(0)=−J′(θ0) <script type="math/tex" id="MathJax-Element-124">r_0=d_{(0)} = -J'(\theta_0)</script>
- 对于k = 0,1,2….
- 利用line search找出使得 J(θk+αid(k)) <script type="math/tex" id="MathJax-Element-125">J(\theta_k+\alpha_i d_{(k)})</script>足够小的 αk <script type="math/tex" id="MathJax-Element-126">\alpha_k</script>
- θk+1=θk+αkd(k) <script type="math/tex" id="MathJax-Element-127">\theta_{k+1} = \theta_k + \alpha_k d_{(k)}</script>
- rk+1=−J′(θk+1) <script type="math/tex" id="MathJax-Element-128">r_{k+1} = -J'(\theta_{k+1})</script>
- d(k+1)=rk+1+βk+1d(k) <script type="math/tex" id="MathJax-Element-129">d_{(k+1)} = r_{k+1} + \beta_{k+1} d_{(k)}</script>
这里又出现了 βi <script type="math/tex" id="MathJax-Element-130">\beta_i</script>,对于 β <script type="math/tex" id="MathJax-Element-131">\beta</script>的研究，著名的方法有:Hestenes-Stiefel方法、Fletcher-Reeves方法、Polak-Ribiére-Polyak 方法和Dai-Yuan方法，分别对应于：
βHSk=(rk+1−rk)Trk+1dT(k)(rk+1−rk)
<script type="math/tex; mode=display" id="MathJax-Element-132">\beta_k^{HS}=\frac{(r_{k+1}-r_k)^Tr_{k+1}}{d_{(k)}^T(r_{k+1}-r_k)}</script>
βFRk=||rk+1||2||rk||2
<script type="math/tex; mode=display" id="MathJax-Element-133">\beta_k^{FR}=\frac{||r_{k+1}||^2}{||r_k||^2}</script>
βPRPk=(rk+1−rk)Trk+1||rk||2
<script type="math/tex; mode=display" id="MathJax-Element-134">\beta_k^{PRP}=\frac{(r_{k+1}-r_k)^Tr_{k+1}}{||r_k||^2}</script>
βDYk=||rk+1||2dT(k)(rk+1−rk)
<script type="math/tex; mode=display" id="MathJax-Element-135">\beta_k^{DY}=\frac{||r_{k+1}||^2}{d_{(k)}^T(r_{k+1}-r_k)}</script>
需要注意的是，非线性共轭梯度法并不能像解决线性系统那样，保证
n <script type="math/tex" id="MathJax-Element-136">n</script>步内收敛，一般我们迭代很多次直到
||rk||<ϵ||r0||<script type="math/tex" id="MathJax-Element-137">||r_{k}|| < \epsilon ||r_0||</script>。
像CG这样高效的方法，一般都有现成的工具库可以使用，只要我们提供目标函数的一次导函数和初始值，CG就能帮我们找到我们想要的了！再次推荐《
An Introduction to the Conjugate Gradient Method Without the Agonizing Pain》一文。
近似牛顿法(Quasi-Newton)
上一节中介绍开源软件包常见的方法Conjugate Gradient Method，这一节我们来介绍另一个常见的方法——Quasi-Newton Method。
牛顿法(Newton Method)
我们高中的时候数学课本上介绍过牛顿求根法，具体的做法是：对于一个连续可导的函数 f(x) <script type="math/tex" id="MathJax-Element-138">f(x)</script>，我们如何求取它的零点呢，看看维基百科是如何展示牛老师的方法:
如图所示，我们首先随机初始化 x0 <script type="math/tex" id="MathJax-Element-139">x_0</script>，然后每一次利用曲线在当前 xi <script type="math/tex" id="MathJax-Element-140">x_i</script>的切线与横轴的交点作为下一个尝试点 xi+1 <script type="math/tex" id="MathJax-Element-141">x_{i+1}</script>，具体更新公式：
xi+1=xi−f(xi)f′(xi)
<script type="math/tex; mode=display" id="MathJax-Element-142">x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)}</script>，直到
f(xi)≈0 <script type="math/tex" id="MathJax-Element-143">f(x_i)\approx0</script>为止。牛老师告诉我们一个求取函数0点的方法，那么对于我们本篇的优化问题有什么帮助呢，我们知道，函数的极值在于导数为0的点取得，那么我们可以利用牛老师的方法求得导数为0的点啊。我们目的是求取
J′(θ)=0 <script type="math/tex" id="MathJax-Element-144">J'(\theta)=0</script>对应的
θ <script type="math/tex" id="MathJax-Element-145">\theta</script>，那么我们可以依样画葫芦（假设
J(θ) <script type="math/tex" id="MathJax-Element-146">J(\theta)</script>是二阶可导的）按照如下更新：
- While J′(θ) <script type="math/tex" id="MathJax-Element-147">J'(\theta)</script>没有足够小：
- θ=θ−J′(θ)J′′(θ) <script type="math/tex" id="MathJax-Element-148">\theta = \theta - \frac{J'(\theta)}{J''(\theta)}</script>
其中 J′′(θ) <script type="math/tex" id="MathJax-Element-149">J''(\theta)</script>是一个矩阵：
J′′(θ)=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪∂J(θ)∂θ0∂θ0...∂J(θ)∂θn∂θ0∂J(θ)∂θ0∂θ1...∂J(θ)∂θn∂θ1.........∂J(θ)∂θ0∂θn∂J(θ)∂θn∂θn⎫⎭⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪=H(2)
<script type="math/tex; mode=display" id="MathJax-Element-150">J''(\theta)= \left\{ \begin{matrix} \frac{\partial J(\theta)}{\partial \theta_0\partial \theta_0} & \frac{\partial J(\theta)}{\partial \theta_0\partial \theta_1}&...& \frac{\partial J(\theta)}{\partial \theta_0\partial \theta_n} \\ . & . & . \\ . & . & . \\ . & . & . \\ \frac{\partial J(\theta)}{\partial \theta_n\partial \theta_0} & \frac{\partial J(\theta)}{\partial \theta_n\partial \theta_1}&...& \frac{\partial J(\theta)}{\partial \theta_n\partial \theta_n} \end{matrix} \right\} \tag{2}=H</script>也就是大名鼎鼎的Hession矩阵。而牛顿法更新中:
θ=θ−J′(θ)J′′(θ)=θ−ηH−1J′(θ)(实际应用我们会加入步伐η，可以用linesearch求得最优的η)
<script type="math/tex; mode=display" id="MathJax-Element-151">\theta = \theta - \frac{J'(\theta)}{J''(\theta)} = \theta - \eta H^{-1}J'(\theta)\\(实际应用我们会加入步伐\eta，可以用line search求得最优的\eta)</script>涉及到Hession矩阵的求逆过程，对于一些参数比较多的模型，这个矩阵将非常巨大，计算也极其耗时，所以这就为什么实际项目中很少直接使用牛老师的方法。不过之前我们介绍的方法都只利用了一阶信息，牛老师的方法启发了利用二阶信息优化方式。
L-BFGS算法
上一小节中，我们介绍了牛顿法，并且指出它一个严重的缺陷，就是计算Hession矩阵和求逆有时候内存和时间都不允许。那么有什么办法可以近似利用牛顿法呢，也就是有没有Quasi-Newton Method呢？答案是有的，BFGS算法就是一个比较著名的近似牛顿法，对于BFGS的介绍，另外有一篇博客有很好的介绍，具体参阅《Numerical Optimization: Understanding L-BFGS》，也是非常直观简洁的介绍，还附有Java和Scala源码，非常值得学习。
BFGS算法核心在于他利用迭代的方式(具体方式请参考上面推荐的博文，文章不长，可读性很强)近似求解Hession矩阵的逆，使得求解Hession矩阵的逆变得不再是神话。而迭代的过程步骤是无限制的，这也会导致内存不足问题，所以工程上利用有限步骤来近似BFGS求解Hession的逆，就成了Limit-BFGS算法。与很多算法一样，这个算法名字是取4位发明者的名字首字母命名的，所以单看名字是没有意义的:-)。
以上是几位大佬的尊荣。利用Quasi-Newton法，在处理数据规模不大的算法模型，比如Logistic Regression，可以很快收敛，是所有优化算法包不可或缺的利器。
参考引用
15
0
已为社区贡献3条内容
所有评论(0)