基于线搜索的下降算法基本思路
- 给定初始点\(x^0\),k=0;
- 判断\(x^k\)是否满足终止条件:是,则终止;
- 寻找\(x^k\)处的下降方向\(d^k\);
- 选择合适的步长\(α_k\)>0,使得;
- 令;令k=k+1;转第1步。
这里的终止条件一般为∇f(\(x^k\))=0,但是在实际计算的时候,我们不会直接判断为0,而是二范数小于一个非常小的数,如\(10^{-4}\)或者\(10^{-6}\)
\(||∇f(x^k)||_2\)≤ɛ,这个ɛ就是这个非常小的数。
如果f(x)是凸函数的时候,就已经找到了最优解;如果f(x)是非凸函数的时候,是一个困难问题,此时我们依然会使用梯度为0为终止条件,这样至少可以找到一个平稳点,它有可能是一个局部解。
终止条件也可以是如果\(||x^k-x^{k+N}||_2\)≤ɛ,说明又经过了N次迭代(N会根据实际的问题来取),并没有发生太大的改变,依然还是在一个非常小的邻域中,此时也可以终止。又有f(\(x^k\))-f(\(x^{k+N}\))≤ɛ,说明经过了N次迭代,函数值的改变是非常小的,此时也可以终止。
下降方向:可以使用负梯度方向,-∇f(\(x^k\)),这里称为最速下降法。它的收敛速度会比较慢。负梯度方向不是唯一选择,还有牛顿方向,共轭梯度方向等等。
步长:可以将α代入到函数中,f(\(x^k\)+α\(d^k\))=Φ(α),这是一个关于α的一元函数,这里就是求该一元函数最小
min Φ(α),α≥0
我们在确定这个\(α_k\)的过程就称为一维线搜索问题。
收敛性问题:点列首先得是收敛的,不能是发散的。不同的算法产生的点列,对于评价算法好坏会有一个“收敛速度”的概念。
线搜索方法
求解一元问题,其解记为\(α^*\).
- 如果f(x)为简单函数,如,如何线搜索?
这是一个二次函数,为二次项,H为正定的海塞矩阵;为线性项;b为常数。
由于H正定,则f(x)为一个凸函数,将f(x)代入 Φ(α)有
Φ(α)=\(1\over 2\)\((x^k+αd^k)^TH\)(\(x^k+αd^k\))+\(c^T\) (\(x^k+αd^k\))+b
=\(1\over 2\)\((d^k)^TH\)\(d^k\)\(α^2\)+(\((d^k)^TH\)\(x^k\)+\(c^T\)\(d^k\))α+b
进行线搜索,就是求
min Φ(α),由于α的二次项系数是大于0的,故这是一个开口向上的一元二次函数求最小,则直接导数为0即可。
Φ'(α)=\((d^k)^TH\)\(d^k\)α+ \((d^k)^TH\)\(x^k\)+\(c^T\)\(d^k\)=0
\(α^*\)=-\((d^k)^THx^k+c^Td^k\over (d^k)^THd^k\)
基于搜索区间的直接搜索法
我们的问题是min Φ(α),α>0
- 什么是搜索区间?包含\(α^*\);单谷(single basin);记为[\(a_0,b_0\)]
单谷是在一个区间中,Φ(α)函数图像形如上图,先单调递减,再单调递增。这样的一个区间 [\(a_0,b_0\)]称为搜索区间。
- 特征:选取λ,μ∈(\(a_0,b_0\))且λ<μ
- 若Φ(λ)<Φ(μ),则得新搜索区间为
- 若Φ(λ)≥Φ(μ),则得新搜索区间为
通过上图,我们可以发现,在Φ(λ)<Φ(μ)时,最小点一定不在[μ,\(b_0\)]区间,是肯定可以剔除的。新的搜索区间则为[\(a_0,μ\)]。
而在Φ(λ)≥Φ(μ)时, 最小点一定不在[\(a_0\),λ]区间,是肯定可以剔除的。新的搜索区间则为[\(λ,b_0\)]。
常用的直接搜索法
- 均匀搜索法(uniform search)
- 令先将区间[\(a_0,b_0\)]N等分;
- 比较相邻三个点对应的函数值,若对于某个i有,则
- 得新搜索区间[\(a_1,b_1\)]
这样我们可以将搜索区间变为原来的\(2\over N\),共计算了N-1个点的函数值。N的取值会根据具体问题的不同而不同。
- 黄金区间法(0.618法)
- 令
- 若,则,产生新搜索区间;
这里\(μ-a_0\over b_0-a_0\)=\(b_0-λ\over b_0-a_0\)=\(λ-a_0\over μ-a_0\)=0.618
第一次插入λ和μ,如果Φ(λ)<Φ(μ),则可以去掉[μ,\(b_0\)]区间,新搜索区间变为[\(a_0\),μ]。这里算了两个点的函数值。继续往下迭代,此时\(a_0\)变成了\(a_1\),μ变成了\(b_1\)。
新插入的两个点,我们会发现λ可以保持不变,只需要重新插入一个μ‘即可。这样就只需要再计算一个点的函数值即可,这是黄金区间法的好处。
- 基于导数信息的二分法
- 记区间中点,计算该点导数值Φ’(λ);
- 若Φ’(λ)=0,则\(α^*\)=λ;
- 若Φ’(λ)<0,则\(α^*\)∈[λ,\(b_0\)],此时λ处于函数值单调递减的部分;
- 若Φ’(λ)>0,则\(α^*\)∈[\(a_0\),λ],此时λ处于函数值单调递增的部分;
非精确线搜索(Inexact linear search)
- Armijo条件:
这是一个关于α的线性函数。我们知道Φ(0)+αΦ'(0)是Φ(α)在0这个点的切线。那么Φ(0)+0αΦ'(0)是一个常数,切线可以看成为Φ(0)+1αΦ'(0);\(c_1\)是一个介于0到1之间的数。
在上图中,平行于横轴α的虚线就是Φ(0)+0αΦ'(0)这个常数,Φ(α)在0上的切线虚线就是Φ(0)+1αΦ'(0);那么Φ(0)+\(c_1\)αΦ'(0)就是介于这两条虚线之间的经过Φ(0)的任意一条线,这条线与Φ(α)的非Φ(0)交点,横坐标即为上图中的实心点。那么α就在0到该实心点的区间范围内。我们在此区间内任选一点就可以作为\(α^*\)的选择。
但是这样的选取条件并不完美,这样选取的α(步长)是为了让目标函数有一个充分的下降,但是α可以选取非常小几乎接近于0的值,这样步长就非常的小,那么\(x^k\)和\(x^{k+1}\)并不会有明显的改变,所以我们希望将这种非常小的步长给排除掉,则有下面的法则。
- Goldstein法则:除Armijo条件外还要求μ满足:
该法则的意思是在满足Armijo条件的基础上,还需要大于等于 Φ(0)+\(c_2\)αΦ'(0)这条直线,即该直线与Φ(α)的非Φ(0)交点,它的横坐标即上图的△点,该点到实心点的范围就是α的范围,在此范围内任意选取一个α都可以作为\(α^*\)的选择。
收敛性:基于线搜索的迭代下降算法
- 为了保证全局收敛(在迭代下降算法的初始点的选取,初始点随便选,都能保证所产生的点列是满足收敛性的),步长\(α_k\)和下降方向\(d^k\)的选取均需满足一定的条件。
- 记当前点\(x^k\)处所选取下降方向\(d^k\)与负梯度方向-∇f(\(x^k\))的夹角为\(θ_k\):
表示两个向量的内积,一般来讲,我们选取的这两个向量的夹角是个锐角,则>0,故cos\(θ_k\)∈[0,\(π\over 2\))
- Zoutendijk定理说明:
- 在适当条件下,迭代点列{\(x^k\)}满足:这里是一个有界值,也是Zoutendijk条件。
- 当每次选取\(d^k\)时均保证\(θ_k\)满足cos\(θ_k\)≥δ(比较小的数),其中δ>0,则||∇f(\(x^k\))||->0。这里是能找到一个梯度为0的点,它是逐渐收敛为0的。我们将\(δ^2\)替代\(cos^2θ_k\),有\(δ^2\)\( \sum_{k=1}^{\infty}\)\({||∇f(x^k)||}^2\)<\(\infty\),如果单独求范数,则肯定有||∇f(\(x^k\))||->0。其实这里说的是在选择\(d^k\)的时候,不要让其与负梯度方向的夹角特别接近\(π\over 2\),即不要太接近垂直的位置,尽量让夹角小一点。
- Zoutendijk定理
假设函数f(x)有下界,且∇f(x)是Lipschitz连续的,即存在L>0使得
记迭代过程为,其中\(d^k\)为下降方向,步长\(α_k\)有1)精确线搜索或者2)Goldstein法则确定,则
证明:由,则
=\({(∇f(x^k)^Td^k)}^2\over {||d^k||}^2\)
已知
由1)精确线搜索,有
\(α_k\)=arg \(min_{α≥0}\) f(\(x^k+αd^k\)) \(α_k\)是所有α的最小值,则有
∀ α>0,均有f(\(x^{k+1}\))=f(\(x^k\)+\(α_kd^k\))≤f(\(x^k\)+α\(d^k\))
根据拉格朗日中值定理(可以参考高等数学整理 中的拉格朗日中值定理),有
f(\(x^k\)+α\(d^k\))=f(\(x^k\))+α∇\(f(ξ^k)^T\)\(d^k\),这里\(ξ^k\)=\(x^k\)+tα\(d^k\),t∈(0,1),是\(x^k\)到\(x^{k+1}\)上的一个点
= f(\(x^k\))+α∇\(f(x^k)^T\)\(d^k\)+α\((∇f(ξ^k)-∇f(x^k))^Td^k\)
≤f(\(x^k\))+α( ∇\(f(x^k)^T\)\(d^k\) +||\(∇f(ξ^k)-∇f(x^k)\)||⋅||\(d^k\)||)
= f(\(x^k\))+α||\(d^k\)||(\(∇f(x^k)^Td^k\over ||d^k||\)+ ||\(∇f(ξ^k)-∇f(x^k)\)||)
已知
令\(a_0\)=-\(∇f(x^k)^Td^k\over ||d^k||\),取\(α_0\)=\(1\over 2L\)⋅\(a_0\over ||d^k||\),可知
f(\(x^{k+1}\))≤f(\(x^k\))+\(α_0\)||\(d^k\)||⋅(-\(a_0\)+ ||\(∇f(ξ^k)-∇f(x^k)\)||)
≤f(\(x^k\))+\(1\over 2L\)\(a_0\)(-\(a_0\)+L||\(ξ^k-x^k\)||)
≤f(\(x^k\))+\(1\over 2L\)\(a_0\) ( -\(a_0\) +L\(α_0\)||\(d^k\)||)
= f(\(x^k\))+\(1\over 2L\)\(a_0\) ( -\(a_0\) +\(1\over 2\)\(a_0\))
=f(\(x^k\))-\(1\over 4L\)\(a_0^2\)= f(\(x^k\))-\(1\over 4L\)\(({-∇f(x^k)^Td^k\over ||d^k||⋅||∇f(x^k)||})^2\)\(||∇f(x^k)||^2\)
= f(\(x^k\))-\(1\over 4L\)\(cos^2θ_k||∇f(x^k)||^2\)
即 f(\(x^{k+1}\))≤ f(\(x^k\))-\(1\over 4L\)\(cos^2θ_k||∇f(x^k)||^2\)
得 f(\(x^{K+1}\))≤f(\(x^0\))-\(1\over 4L\)\(\sum_{k=0}^K\)\(cos^2θ_k||∇f(x^k)||^2\)
令K->\(\infty\),f(\(x^*\))≤ f(\(x^0\))-\(1\over 4L\)\(\sum_{k=0}^\infty\)\(cos^2θ_k||∇f(x^k)||^2\)
得得证
由2)Goldstein法则,有
则有
f(\(x^{k+1}\))≤f(\(x^k\))+\(c_1α_k∇f(x^k)^Td^k\) (1)
f(\(x^{k+1}\))≥f(\(x^k\))+\(c_2α_k∇f(x^k)^Td^k\) (2)
0<\(c_1\)<\(c_2\)<1
根据拉格朗日中值定理,有
f(\(x^{k+1}\))=f(\(x^k\))+\(α_k\)∇f\((ξ^k)^T\)\(d^k\) 这里\(ξ^k\)=\(x^k\)+t\(α_k\)\(d^k\),t∈(0,1)
= f(\(x^k\))+\(α_k\)∇\(f(x^k)^T\)\(d^k\)+\(α_k\)\((∇f(ξ^k)-∇f(x^k))^Td^k\)
由(2)有 ≥f(\(x^k\))+\(c_2α_k∇f(x^k)^Td^k\)
得 \((∇f(ξ^k)-∇f(x^k))^Td^k\) ≥(\(c_2\)-1)∇\(f(x^k)^T\)\(d^k\)
又 \((∇f(ξ^k)-∇f(x^k))^Td^k\) ≤||\(∇f(ξ^k)-∇f(x^k)\)||⋅||\(d^k\)||
≤L⋅||\(ξ^k-x^k\)||⋅||\(d^k\)||
≤L⋅\(α^k\)\(||d^k||^2\)
得 L⋅\(α^k\)\(||d^k||^2\)≥ (\(c_2\)-1)∇\(f(x^k)^T\)\(d^k\)
\(α^k\)≥\((c_2-1)\over L\)⋅\(∇f(x^k)^Td^k\over ||d^k||^2\)
由(1) f(\(x^{k+1}\))≤f(\(x^k\))+\(c_1α_k∇f(x^k)^Td^k\)
≤ f(\(x^k\))+\(c_1∇f(x^k)^Td^k\)⋅ \((c_2-1)\over L\)⋅\(∇f(x^k)^Td^k\over ||d^k||^2\)
=f(\(x^k\))+\(c_1(c_2-1)\over L\)⋅\((∇f(x^k)^Td^k)^2\over ||d^k||^2\)
=f(\(x^k\))+\(c_1(c_2-1)\over L\)\(cos^2θ_k||∇f(x^k)||^2\)
令c=|\(c_1(c_2-1)\over L\)| =f(\(x^k\))-c⋅\(cos^2θ_k||∇f(x^k)||^2\)
再根据精确搜索相同的方法,最终可得
常用求解方法
- 坐标轴交替下降法
- 最速下降法
- 牛顿法
- 拟牛顿法
这些方法的不同主要在于下降方向选择策略的不同。
- 收敛速度
设序列{\(x^k\)}收敛到\(x^*\),若存在极限
- 当0<β<1时,则称 {\(x^k\)}为线性收敛。
- 当β=0,则称 {\(x^k\)}为超线性收敛。
若存在某个p≥1,有
则称{\(x^k\)}为p阶收敛。
当p>1时,p阶收敛必为超线性收敛,但反之不一定成立。
- 朴素的算法:坐标轴交替下降法
基本思想:给定初始点\(x^0\),依次沿着坐标轴\(e_1,...,e_n\)进行搜索:
上图中,假定是一个二维空间,横坐标表示\(e_1\),纵坐标表示\(e_2\)。初始点\(x^0\)先沿着\(e_1\)的正方向走一点,再沿着\(e_2\)的正方向走一点,得到\(x^1\),此时函数值是下降的。再沿着\(e_1\)正方向走一点,再沿着\(e_2\)负方向走一点,得到\(x^2\),此时函数值是下降的,如此反复。
框架(Framework):
- 给定初始点\(x^0\),k=0,c>0;
- 判断是否满足;是,终止;
- 记,令其中
- 令转步1;
优点:不需成本即可获得搜索方向;当目标函数变量之间交叉程度较小(loosely coupled)时非常有效(极端情况如可分离函数);
f(x)=\(f_1(x_1)+f_2(x_2)+...+f_n(x_n)\)变量之间只属于各自的函数,各自的函数相加组成总函数为可分离函数。
缺点:对于一般问题所得点列未必收敛;
改进的方法:
任意取一点\(x^0\)沿着\(x_1\)和\(x_2\)方向找出一点,由于我们知道f()<f(\(x^0\)),故在\(x^0\)指向的方向上,目标函数值是下降的,沿着\(x^0\)指向的方向上进行线搜索,找出一个步长进行推进,得到新的点\(x^1\)。此时又回到坐标轴交替下降的过程,找出一点,再沿着\(x^1\)指向的方向进行推进,得到\(x^2\),以此往复。该方法在一些特殊的问题中可以明显的改进收敛的情况。
- 最速下降法(也称为梯度下降法)
基本思想:选择\(x^k\)处负梯度作为搜索方向,即
如果给定一个方向d,\(∇f(x^k)^Td\)<0,则d必然是下降方向,而\(∇f(x^k)^Td\)称为方向导数。
特点:简单直观;收敛;搜索方向只需计算∇f(\(x^k\));
缺点:
- 收敛速度慢(线性收敛);
- Zigzag现象;
- 不具备二次终止性(即在有限步内求得凸二次函数最优解);
原因:若迭代中步长\(α_k\)是的精确最小点,则,即
这意味着点列中新的点\(x^{k+1}\)的梯度方向与原始点\(x^k\)的梯度方向是垂直的。
例:设初始点为
这是一个凸二次函数,最优解为\((0,0)^T\),蓝色虚线是f(x)的等值线。\(x^0\)是最外圈的蓝色虚线与红线的交点。\(x^0\)的负梯度方向就是该红线指向内部第二圈蓝色虚线的方向。该红线与内部第二圈蓝色虚线的交点为\(x^1\)。\(x^1\)的负梯度方向就是该点红线指向内部第三圈蓝色虚线的方向,与第三圈蓝色虚线的交点为\(x^2\)。那么\(x^0x^1\)线段与\(x^1x^2\)线段的夹角是垂直\(90^\circ\)的。最速下降法会不断产生点列,不断接近中心最小点,但是不能够在有限步内取到。
- 牛顿法
基本思想:当前点\(x^k\)处选择方向的时候不仅考虑梯度的信息(一阶导数),还要考虑二阶导数,也就是海塞矩阵的信息可理解为:对\(x^k\)处的二次逼近函数(二阶泰勒展式)进行最小化;
上式中为常数项,为线性项,为二次项,其中为海塞矩阵。
二次函数求最小,一定是梯度为0为最小。即上式就变成求
∇f(\(x^k\))+\(∇^2f(x^k)(x-x^k)\)=0
x=\(x^k\)-\([∇^2f(x^k)]^{-1}∇f(x^k)\)
我们令x为\(x^{k+1}\),上式就为
\(x^{k+1}\)= \(x^k\)-\([∇^2f(x^k)]^{-1}∇f(x^k)\)
牛顿步(Newton step):
在\(x^k\)的位置,选择方向为,则步长\(α_k\)=1,此为纯牛顿法(Pure Newton method);
步骤:
- \(x^0\),ε很小,k=0
- 如果梯度||∇f(\(x^k\))||≤ε,终止,否则
- 计算牛顿方向
- \(x^{k+1}\)=\(x^k\)+\(d^k\),k=k+1,返回step 1
纯牛顿法只在一些前提下有效,即在\(x^k\)的位置,它的海塞矩阵正定\(∇^2f(x^k)\)0,则整个二次函数为凸函数,梯度为0为最小点,而\(d^k\)肯定是下降方向的。
如果 在\(x^k\)的位置,它的海塞矩阵\(∇^2f(x^k)\)非正定,有可能特征值存在0和负数。如果特征值λ存在0,则行列式不存在,逆矩阵无法求,牛顿方向无法计算。如果特征值λ<0,我们可以看一下梯度向量与\(d^k\)的内积(有关行列式和特征值的内容可以参考线性代数整理(三) 中的行列式和特征值和特征向量)
∇\(f(x^k)^Td^k\)=-∇\(f(x^k)^T\)\([∇^2f(x^k)]^{-1}\)∇f(\(x^k\))
此时我们无法判断上面的式子小于0,所以\(d^k\)未必是下降方向。
优点:当初始点\(x^0\)取的比较接近于收敛点\(x^*\),且\(∇^2f(x)\)满足较好性质(正定)时,二阶收敛(收敛速度快);二次终止性,牛顿法求解凸二次最小化问题,一步可得最优解(因为牛顿法本身就是二次近似函数,跟目标函数相同)。
缺点:计算量大(需计算Hesse矩阵);适用范围较窄;
例:设初始点为
对于凸二次目标函数,一步收敛:
- 修正牛顿法(Modified Newton method)
- 对于步长\(α_k\)的修正:首先判断\(α_k\)=1是否让目标函数充分下降;否,则采用线搜索方法重新确定\(α_k\);
- 对于方向(Hesse矩阵)的修正:选取若,(不需要修正)则选取否则,采取修正方法(多种):
- ,λ为适当正数保证\(B_k\)正定;也就是说当海塞矩阵非正定的时候,它的特征值会有0或负值的\(λ_1,...,λ_n\)中,\(λ_i\)+λ>0,i=1,...,n,即λ>\(max_{i=1...n}\){-\(λ_i\)}。这样就对对称矩阵的主对角线扰动,变成了正定矩阵。
- 考虑特征值分解其中令δ为一适当的正数。
在对海塞矩阵进行特征值分解的时候(有关特征值分解的内容可以参考线性代数整理(三) 的正交对角化),会有一系列的对角特征值,其中会有小于0或者等于0的数
我们将其中的负数使用一些正数δ进行替代,而本身的正数不变
- 拟牛顿法(Quasi-Newton method)
考虑f(x)在当前点\(x^k\)处的二次近似函数
其中
这里的意思是相对于纯牛顿法,这里在进行二次近似的时候是直接使用一个正定矩阵来代替海塞矩阵\(∇^2f(x^k)\)。对于\(B_k\)而言,我们希望它能体现一些二阶信息,同时它的获取的代价要小一点。
利用min \(m_k\)(x),梯度为0,得搜索方向
Quasi-Newton Framework:
- \(x^0\) ,ε很小,k=0,\(B_0\),可以对\(x^0\)的海塞矩阵进行修正(如果海塞矩阵非正定)得到\(B_0\)
- 如果梯度||∇f(\(x^k\))||≤ε,终止,否则
- 计算方向
- 使用线搜索法确定步长\(α_k\)
- \(x^{k+1}\)=\(x^k\)+\(α_k\)\(d^k\),确定\(B_{k+1}\),k=k+1,返回step 1
- 如何比较简便获得矩阵\(B_{k+1}\)?
拟牛顿方程(基本要求):
我们的目的是需要在\(x^{k+1}\)处确定一个\(B_{k+1}\)矩阵。根据二阶的拉格朗日中值定理有
\(∇f(x^{k+1})-∇f(x^k)\)=\(∇^2f(ξ)\)(\(x^{k+1}-x^k\)),其中ξ=λ\(x^k+(1-λ)x^{k+1}\),λ∈(0,1)
通过比较上式和拟牛顿方程,我们会发现
\(B_{k+1}\)就是为了代替这个海塞矩阵\(∇^2f(ξ)\)。在拟牛顿方程中,\(∇f(x^{k+1})\)、\(∇f(x^k)\)、\(x^{k+1}\)、\(x^k\)都是已知的,\(B_{k+1}\)是n*n阶的方阵,它里面有\(n(n+1)\over 2\)个元素待定,拟牛顿方程包含了n个方程组成方程组。
记拟牛顿方程简写为
若记拟牛顿方程也表示为
满足拟牛顿方程的矩阵很多!具体获取\(B_{k+1}\)方法如下
基于已有信息获取\(B_{k+1}\),或基于已有信息获取\(H_{k+1}\)
- 第一类方法:选择满足拟牛顿方程且与\(B_k\)近似的矩阵:当然也可以使用H矩阵min ||H-\(H_k\)||,s.t. H\(y_k\)=\(s_k\),H=\(H^T\)
- 第二类方法:对\(B_k\)(或者\(H_k\))进行校正;如,令
- rank-2校正,要求ΔB的秩为2:DFP方法,BFGS方法;(有关秩的内容可以参考线性代数整理(二) 中的矩阵的秩)
- rank-1校正,要求ΔB的秩为1:SR-1方法。
对于第二类方法,也可以对\(H_{k+1}=H_k+ΔH\),要求ΔH的秩为1或2。
- 拟牛顿之DFP(Davidon-Fletcher-Powell)方法
DFP方法可看作是对\(H_k\)进行rank-2校正,
秩1矩阵可以写成一个向量u乘以一个向量\(v^T\),即u\(v^T\),由于H是对称矩阵,所以这里为u\(u^T\),它是一个秩1的方阵;秩2的矩阵为两个秩1的矩阵相加,即u\(u^T\)+v\(v^T\)。
则秩2校正有
\(H_{k+1}\)=\(H_k\)+a u\(u^T\)+b v\(v^T\)
其中a,b为标量系数,u,v∈\(R^n\)
由拟牛顿方程有\(H_{k+1}⋅y_k=s_k\),代入得
\(H_ky_k+auu^Ty_k+bvv^Ty_k=s_k\)
\(H_ky_k+auu^Ty_k+bvv^Ty_k-s_k=0\)
令u=\(H_Ky_k\),得1+a\(u^Ty_k\)=0,即a=-\(1\over u^Ty_k\)=-\(1\over y_k^TH_ky_k\)
令v=\(s_k\),得b\(v^Ty_k\)-1=0,即b=\(1\over v^Ty_k\)=\(1\over s_k^Ty_k\)
将a,b,u,v代入秩2校正,则有
- 拟牛顿之BFGS(Broyden-Fletcher-Goldfarb-Shannon)方法
BFGS方法可看作对\(B_k\)进行rank-2校正,
上式可以通过\(B_{k+1}=B_k+auu^T+bvv^T\)以及拟牛顿方程\(B_{k+1}⋅s_k=y_k\)推导出来。
拟牛顿方向\(d^{k+1}=-B_{k+1}^{-1}∇f(x^{k+1})\)需要计算\(B_{k+1}^{-1}\),可利用Sherman-Morrison公式显示写出:
被认为是最有效的拟牛顿法;超线性收敛;
Broyden族:DFP与BFGS的线性组合!
{λ\(B_{k+1}^{DFP}+(1-λ)B_{k+1}^{BFGS}\)} λ∈[0,1]
- 拟牛顿之SR-1方法
SR-1方法是对\(B_k\)进行rank-1校正,
由于是秩1校正,有
\(B_{k+1}\)=\(B_k+auu^T\)
由拟牛顿方程有
\(B_{k+1}⋅s_k=y_k\)
代入,得
\(B_ks_k+auu^Ts_k=y_k\)
au(\(u^Ts_k\))=\(y_k-B_ks_y\)
令u=\(y_k-B_ks_y\),得a\(u^Ts_k\)=1,即a=\(1\over u^Ts_k\)=\(1\over (y_k-B_ks_k)^Ts_k\)
将a,u代入秩1校正,则有
SR-1方法迭代公式更简单,但不能保证\(B_{k+1}\)矩阵正定性;适当条件下能达到n步超线性收敛。
共轭梯度法
- 背景知识
在1950年代,该方法主要研究的是二次函数的最小值问题
min \(1\over 2\)\(x^T\)Qx+\(c^T\)x,Q是n阶的正定矩阵\(Q_n\)0,c是n维向量c∈\(R^n\)
因为Q正定,则目标函数是一个凸二次函数,求梯度为0即为该问题的解
∇f(x)=Qx+c=0
因为Q是矩阵,c是向量,故上式是一个方程组,当问题的规模n不大的时候可以使用高斯消元法来解(有关高斯消元法的内容可以参考线性代数整理 中的高斯消元法),如果规模n很大的话,使用高斯消元法的计算成本非常大。
于是有人提出了使用迭代的方法来求解,就是我们这里的共轭梯度法,在这里也称为线性共轭梯度法。
在1960年代,共轭梯度法被推广到了求一般性函数的最小值min f(x),不再局限于二次函数,该方法称为非线形共轭梯度法。
- 线性共轭梯度法
考虑问题
Q对角
这个问题中只有\(x_i\)的平方项,没有\(x_i\)的交叉项。当n=2的时候,f(x)的等值线如下
n=2的时候是一个椭圆,它的长短轴平行于坐标轴;当n≥3的时候是一个椭球面。使用之前说的坐标轴交替下降法,我们任取一个点,沿着横轴走,得到与椭圆等值线的一个交点,再沿着纵轴走即可得到目标函数的最小值。在2步之内就可以找到最小值,如果n≥3,沿着n个轴都走一遍就可以把最小值找出。
Q非对角
目标函数的等值线依然是一个椭球,只不过它的长短轴不再平行于坐标轴。我们可以对其进行正交对角化
Q=\(P^TDP\),D是以特征值组成的对角矩阵,P是特征向量组成的正交矩阵。(有关正交对角化的内容可以参考线性代数整理(三) 中的正交对角化)
我们对原目标函数进行正交替换,令=Px,有
此时使用Q对角的方式,沿着各个轴走,就可以找出的最小值。此时f(x)的最优解\(x^*\)=\(P^T\)
线性共轭梯度法被设计出来就是为了求解线性方程组的,如果要使用Q非对角的方法,就必须求解P正交矩阵,特征向量的获得就得解线性方程组。这个P正交矩阵是没有办法直接使用的,我们的问题就是找出n个n维向量(\(d^0,d^1,...,d^{n-1}\))组成的一个矩阵S。我们希望这个S能够像正交矩阵P一样,能够把Q进行对角化。这是共轭梯度法的一个基本的想法。这里的\(d^0,d^1,...,d^{n-1}\)之间的关系就称为共轭。
- 共轭方向
考虑正定矩阵Q及非零向量\(d^i,d^j\),若
则称\(d^i,d^j\)关于矩阵Q共轭;
向量组\(d^0,d^1,...,d^k\)关于矩阵Q共轭,指的是它们两两共轭;
共轭 vs 正交
- 如果\(d^0,d^1,...,d^k\)关于单位矩阵I共轭,则\(d^0,d^1,...,d^k\)正交;
- 如果\(d^0,d^1,...,d^k\)关于正定矩阵Q共轭,Q=\(P^2\),P正定,则\((d^i)^TQd^j=(d^i)^TP^TPd^j=(Pd^i)^T(Pd^j)\)=0,即如果原来的一个向量组关于Q矩阵共轭,那么使用P矩阵做变换以后的向量Pd为正交关系。
- 正交向量组是线性无关的,共轭向量组也是线性无关的。
- 共轭向量组最多只有n个向量\(d^0,d^1,...,d^{n-1}\)
- 共轭方向法
给定初始点\(x^0\)及一组关于Q共轭方向\(d^0,d^1,...,d^{n-1}\),令
其中
由于f(x)本身是一个凸二次函数,φ(α)的最小值即为求
φ‘(\(α_k\))=\(∇f(x^k+α_kd^k)^Td^k\)=\((Q(x^k+α_kd^k)+c)^Td^k\)=\((Qx^k+c)^Td^k+α_k(d^k)^TQd^k\)=0
计算得
Remark:这就是共轭方向法,我们需要给出这一组共轭方向\(d^0,d^1,...,d^{n-1}\),共轭方向法为一类方法,共轭梯度法是其中一种。
几何解释:
从\(x^0\)沿着\(d^0\)方向走\(α_0\)这么远,再沿着\(d^1\)方向走\(α_1\)这么远,直到最后沿着\(d^{n-1}\)方向走\(α_{n-1}\)这么远
\(x^0\)+\(α_0\)\(d^0\)+\(α_1\)\(d^1\)+...+\(α_{n-1}\)\(d^{n-1}\)
我们知道\(d^0,d^1,...,d^{n-1}\)是n个线性无关的n维向量,把它们组成一个矩阵S,则S为方阵,且可逆;
- \(S^TQS\)=,Q0;它的最终结果是一个n*n的方阵,如果i≠j的时候,由于\(d^i,d^j\)共轭,则\((d^i)^TQd^j\)=0;当i=j的时候,由于Q正定,\((d^i)^TQd^j\)>0;\(S^TQS\)整个是一个对角阵。
- \(S^{-1}\)(\(d^0,d^1,...,d^{n-1}\))=(\(S^{-1}\)\(d^0\),\(S^{-1}\)\(d^1\),...,\(S^{-1}\)\(d^{n-1}\))=I,由此可以看出\(S^{-1}\)\(d^i\)是一个单位坐标向量\(e^{i+1}\),i=0,...,n-1。
我们的目标函数为
令x=S或者写为=\(S^{-1}\)x
则目标函数可以转换成
其中\(S^TQS\)为对角矩阵,意味着二次项没有交叉项,等值面是与坐标轴平行的椭球面。
\(x^0\)+\(α_0\)\(d^0\)+\(α_1\)\(d^1\)+...+\(α_{n-1}\)\(d^{n-1}\)的迭代过程可以写为
\(S^{-1}\)\(x^0\)+\(α_0\)\(S^{-1}\)\(d^0\)+\(α_1\)\(S^{-1}\)\(d^1\)+...+\(α_{n-1}\)\(S^{-1}\)\(d^{n-1}\),则转化为
\(S^{-1}\)\(x^0\) +\(α_0\)\(e_1\)+\(α_1\)\(e_2\)+...+\(α_{n-1}\)\(e_n\)
在x空间中沿着共轭方向搜索的过程,其实对应到空间中,是沿着坐标轴去搜索的过程。
重要特征:
考虑问题
给定初始点\(x^0\)及一组关于Q共轭方向\(d^0,d^1,...,d^{n-1}\),令
点列{\(x^k\)}具有如下特征
- 在点列中的某一个点目标函数的梯度跟使用过的共轭方向满足垂直的关系。
- \(x^k\)其中,\(x^k=x^0+α_0d^0+...+α_{k-1}d^{k-1}\),它是目标函数的最优解,\(X^k\)是一个平移的子空间中求最小。当k=n的时候,\(x^n\)=arg min f(x),over \(X^n=R^n\),它说明共轭方向法能够在n步之内找到最优解。
证明1:当i=k-1,∇f\((x^k)^Td^{k-1}\)=0,k=1,...,n
我们知道步长\(α_{k-1}\)=arg min f(\(x^{k-1}+αd^{k-1}\))=φ(α)是对φ(α)求导为0
φ‘(\(α_{k-1}\))=∇f\((x^{k-1}+α_{k-1}d^{k-1})^Td^{k-1}\)= ∇f\((x^k)^Td^{k-1}\)=0
此时我们要看∇f\((x^k)^Td^i\),i=0,1,...,k-2
=\((Qx^k+c)^T\)\(d^i\)
=\((Q(x^{i+1}+α_{i+1}d^{i+1}+...+α_{k-1}d^{k-1})+c)^Td^i\) (\(d^i,d^{i+1},...,d^{k-1}\)共轭)
=\((Qx^{i+1}+c)^Td^i\)
=∇f\((x^{i+1})^Td^i\) (令i+1=k,i=k-1)
=0
1得证
证明2:记φ(\(a_0,a_1,...,a_{k-1}\))=\({1\over 2}(x^0+\sum_{i=0}^{k-1}a_id^i)^TQ(x^0+\sum_{i=0}^{k-1}a_id^i)+c^T(x^0+\sum_{i=0}^{k-1}a_id^i)\)
(\(α_0,α_1,...,α_{k-1}\))=arg min φ(\(a_0,a_1,...,a_{k-1}\))
求偏导
\({∂φ(a_1,a_2,...,a_{k-1})\over ∂a_i}|_{(α_0,α_1,...,α_{k-1})}\)=\((Q(x^0+α_0d^0+...+α_{k-1}d^{k-1})+c)^Td^i\)=\((Qx^k+c)^Td^i\)=∇f\((x^k)^Td^i\)=0,i=0,...,k-1
2得证
- 共轭梯度法(共轭方向法的一种)
在迭代下降过程中,借助当前点\(x^k\)的梯度信息构造共轭方向;
Main framework:
- 给定初始点\(x^0\),记
- 判断是否成立;是,则终止;
- 计算
- 令,并计算方向,保证\(d^{k+1}\)跟前面所有的方向都共轭。令k=k+1,转Step 1;
算法过程中,式具体为
上式中\(d^{k+1}\)与所有\(d^0,d^1,...,d^k\)均关于Q共轭。
这里我们可以很容易想到
\(d^{k+1}\)=-∇f(\(x^{k+1}\))+\(β_0d^0+β_1d^1+...+β_kd^k\)
这么做的目的是为了满足共轭\((d^{k+1})^TQd^i=0\),i=0...k
即 \((-∇f(x^{k+1})^T+β_0d^0+...+β_kd^k)^TQd^i=0\)
由共轭,有 -∇f\((x^{k+1})^TQd^i+β_i(d^i)^TQd^i\)=0
得 \(β_i\)=\(∇f(x^{k+1})^TQd^i\over (d^i)^TQd^i\),i=0...k
当i=0,1,...,k-1的时候,\(β_i\)的分子
\(∇f(x^{k+1})^TQd^i\)
\(d^i\)是方向,我们知道\(α_id^i=x^{i+1}-x^i\),则\(d^i={(x^{i+1}-x^i)}\)\(1\over α_i\),再两边同时乘以Q
\(Qd^i=Q{(x^{i+1}-x^i)}\)\(1\over α_i\)
=\((Qx^{i+1}-Qx^i)\)\(1\over α_i\)
=(∇f(\(x^{i+1}\))-∇f(\(x^i\)))\(1\over α_i\)
\(∇f(x^{k+1})^TQd^i\)=\(∇f(x^{k+1})^T\) (∇f(\(x^{i+1}\))-∇f(\(x^i\)))\(1\over α_i\)
由于
\(d^{i+1}=-∇f(x^{i+1})+"d^0,d^1,...,d^i"\),这里\("d^0,d^1,...,d^i"\)表示之前方向的一个组合
\(∇f(x^{i+1})\)=-\(d^{i+1}+"d^0,d^1,...,d^i"\)
在 \(∇f(x^{k+1})^T\) (∇f(\(x^{i+1}\))-∇f(\(x^i\)))\(1\over α_i\)中
\(∇f(x^{k+1})^T∇f(x^{i+1})\)=\(∇f(x^{k+1})^T\)(\(-d^{i+1}+"d^0,d^1,...,d^i"\))
根据共轭方向法的特征1, \(∇f(x^{k+1})^T\)(\(-d^{i+1}+"d^0,d^1,...,d^i"\))=0
同样\(∇f(x^{k+1})^T∇f(x^i)\)=0
故 \(β_i\)的分子\(∇f(x^{k+1})^TQd^i\)=0,故最终只剩下
- 公式简化
共轭梯度法的步长公式
可简化为
这里简化的目的是为了把公式中的表达式跟目标函数的梯度建立起比较直观的联系。
步长公式的分子
∇f\((x^k)^Td^k\)
由于\(d^k=-∇f(x^k)+β_{k-1}d^{k-1}\),则
∇f\((x^k)^Td^k\)=∇f\((x^k)^T\)(\(-∇f(x^k)+β_{k-1}d^{k-1}\))=- ∇f\((x^k)^T\) ∇f(\(x^k\))+\(β_{k-1}\) ∇f\((x^k)^T\) \(d^{k-1}\)
根据特征1,有 \(β_{k-1}\) ∇f\((x^k)^T\) \(d^{k-1}\)=0
故 ∇f\((x^k)^Td^k\)= - ∇f\((x^k)^T\) ∇f(\(x^k\)),故最终\(α_k\)可简化为
------------------------------------------------------------------------------
共轭梯度法步长公式中的系数
可简化为
\(β_k=\)\(∇f(x^{k+1})^T(∇f(x^{k+1})-∇f(x^k))\over ∇f(x^k)^T∇f(x^k)\) (1)
=\(∇f(x^{k+1})^T∇f(x^{k+1})\over ∇f(x^k)^T∇f(x^k)\) (2)
\(β_k\)的分子
\(∇f(x^{k+1})^TQd^k\)=\(∇f(x^{k+1})^T(∇f(x^{k+1})-∇f(x^k))\)\(1\over α_k\),则
=\(1\over α_k\)\(∇f(x^{k+1})^T(∇f(x^{k+1})-∇f(x^k))\over (d^k)^TQd^k\)
将代入上式,就有了
\(β_k=\)\(∇f(x^{k+1})^T(∇f(x^{k+1})-∇f(x^k))\over ∇f(x^k)^T∇f(x^k)\)
由于\(∇f(x^{k+1})^T∇f(x^k)=0\),故上式可化简为
\(β_k=\)\(∇f(x^{k+1})^T∇f(x^{k+1})\over ∇f(x^k)^T∇f(x^k)\)
我们这里仍然在讨论一个二次函数,但是如果把线性共轭梯度法推广到非线性,要求解的是一般性的问题min f(x),我们会发现在简化的\(β_k\)中没有出现二次目标函数中的Q、c,在(1)(2)式中出现的都只有目标函数的梯度,跟目标函数具体的样子没关系,仅要求目标函数可微,梯度可求。当目标函数是一般函数的时候,(1)(2)式的值未必一样,只有在目标函数是二次函数的时候是相等的。
- 非线形共轭梯度法:(FR方法,PRP方法)
- 给定初始点\(x^0\),记
- 判断是否成立;是,则终止;
- 利用线性搜索计算步长\(α_k\);(这里跟线性共轭梯度法不同)
- 令并计算方向其中或者(这两个公式只能二选一,并保持一致性)
- 令k=k+1,转Step 1。
一些说明
- 在实践中,为保证每次产生的方向\(d^{k+1}\)为下降方向,可能对\(β_k\)进行调整;
- 具有二次终止性;
- 实现过程中常采用n步重启策略,可达到n步二阶收敛\(\lim_{k->∞}{||x^{k+n}-x^*||\over ||x^k-x^*||^2}=γ\);
这里说的是当迭代的次数比较多的时候,离当前点比较远的点\(x^0,x^1\)的信息已经不起作用了,于是在当前点重新设置为负梯度方向,再依次来一遍。另外一种情况是当当前点落入了类似二次函数的阶段时,此时进行重启方向,就相当于是线性共轭梯度法在起作用,它能够在n步之内找到最小点,此时能够发挥线性共轭梯度法的优势。
跟其他的方法相比,它的收敛速度是明显比最速下降法快的;跟牛顿法和拟牛顿法比,牛顿法是需要计算海塞矩阵的,拟牛顿法是需要每一步计算一个B矩阵,共轭梯度法不牵涉任何矩阵的情况,它只是梯度做一些运算,存储量较小,比较适合于当n比较大的大规模问题。