概述
优化问题就是在给定限制条件下寻找目标函数\(f(\mathbf{x}),\mathbf{x}\in\mathbf{R}^{\mathbf{n}}\)的极值点。极值可以分为整体极值或局部极值,整体极值即函数的最大/最小值,局部极值就是函数在有限邻域内的最大/最小值。通常都希望能求得函数的整体极值,但总体来说这是很困难的,目前有一些启发式算法可以在某种程度上处理全局极值的计算问题,但是并不能保证每次都能求出整体极值。另外,有一类“凸优化”问题可以比较简单地计算得到全局极值。
本文的目的是对一些常用的优化算法做一个整理和归纳,希望能够尽量理清各类算法的基本思路。算法中涉及到的公式推导和定理证明一般都不太复杂,为了避免过多的公式冲淡本文的主题(公式一个个敲出来太麻烦),所以采取能省则省的原则,如有需要可以查阅本文列出的参考文献。另外,文献[3]中还有部分算法的代码实现(该书被导师誉为数值计算领域的圣经)。为了简单起见,下文介绍的方法都以光滑函数的无约束优化问题为例。另外,由于极大值和极小值是等价的,下文默认待求的是函数极小值。
首先,根据自变量\(\mathbf{x}\)的维度\(\mathbf{n}\)的不同,可以将问题划分为一维和多维两种情况。对于一维优化方法,由于只有一个维度,所以就不存在选择“搜索方向”的问题,只需要在该方向上迭代计算极值即可。而对于多维优化方法,情况就复杂得多,一般的处理方法主要有直线搜索法和置信域法两种。对于我而言,看过的大部分书上都重点介绍的是直线搜索法,置信域法出镜率相对较低。
本文参考文献:
[1] Press, W.H., 2007. Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press.
[2] Boyd, S. and Vandenberghe, L., 2004. Convex optimization. Cambridge university press.
[3] Nocedal, J. and Wright, S., 2006. Numerical optimization. Springer Science & Business Media.
一维优化方法
黄金分割法
黄金分割法求函数极值的思想和求根算法中的二分法类似。在二分法中,首先需要确定一个根所在的区间\((a,b)\),其中\(a\)和\(b\)满足\(ab<0\)以确保至少区间中至少有有一个根。然后逐步缩小该区间直到达到预设精度。在黄金分割法中,同样也需要确定极值所在的范围。这个可以用一个三元组\((a,b,c)\)表示,其中\(a<b<c\)且\(f(b)<min\{f(a),f(c)\}\)。这样,在区间\((a,c)\)内存在一极值。
接下来,需要逐步设法缩小区间\((a,c)\)。在区间\((a,c)\)内选择新的一点\(x\),则\(x\)落在区间\((a,b)\)或\((b,c)\)内。以\(x\in(b,c)\)为例,计算\(f(x)\)的值并与\(f(b)\)作比较。如果\(f(b)<f(x)\),那么取新的三元组为\((a,b,x)\),如果\(f(b)>f(x)\),那么将新的三元组取为\((b,x,c)\)。这样一来,总能保证新的三元组确定的区间长度小于原来的区间。重复这一过程直到区间范围小于预设的精度即可得到一个极值点。
接下来讨论\(x\)的选取方法。令\(b\)分区间\((a,c)\)的比例为\(w\),即
\(\frac{b-a}{c-a}=w\),\(\frac{c-b}{c-a}=1-w\)
而下一实验点\(x\)与\(a,b,c\)的关系为:
\(\frac{x-b}{c-a}=z\)
那么下一区间长度是原区间长度的\(w+z\)或\(1-w\)倍。为了最小化最坏情况出现的概率,令
\(w+z=1-w\)
即\(z=1-2w\),可见\(x\)就是\(b\)点关于区间中心的对称点,所以\(x\)处于\((a,b)\)和\((b,c)\)中较长的哪个区间内。
关于\(w\)的选取,可以采用scale similarity的思想,即每一步的区间划分比例是相等的,如果\(z\)是最佳比例,那么\(w\)也是最佳,所以有
\(\frac{z}{1-w}=w\)
再利用\(z=1-2w\)可以解出\(w\approx 0.38197\),\(1-w\approx 0.61803\),即黄金分割比例。而该算法也因此得名。
下面将黄金分割算法做一个归纳:首先给定一个包含极值点的初始三元组\((a,b,c)\),然后下一个试探点\(x\)选在较长区间段的0.38197分割点。初始三元组并不需要满足黄金分割比例,算法会在迭代过程中自动达到一个稳定的重复比例。显然黄金分割算法是一种线性收敛算法,每一步计算都能使区间范围缩小到原来的0.618。
另外,关于“预设精度”的选取问题,有一个原则就是,不能设置精度为数值的机器精度\(\epsilon\),即区间范围为\((\bar{x}-\epsilon \bar{x},\bar{x}+\epsilon \bar{x} )\),其中\(\bar{x}\)为极值点,这一点和求根算法不同。为了说明该原则,考虑函数\(f(x)\)在极值点的泰勒展开:
\(f(x)\approx f(\bar{x})+\frac{1}{2}f\'\'(\bar{x})(x-\bar{x})^2\) (1)
当
\(|x-\bar{x}|<\sqrt{\epsilon}|\bar{x}|\sqrt{\frac{2|f(\bar{x})|}{\bar{x}^2f\'\'(\bar{x})}}\) (2)
时,式(1)右端可以写成\((1+\pm\epsilon)f(\bar{x})\),其中正负号由\(f(\bar{x})\)的符号决定。而式(2)中最后的平方根项对于绝大多数函数都是单位阶大小的(即可以认为是1)。所以\(|x-\bar{x}|\)最多可以到达\(sqrt{\epsilon}|\bar{x}|\)量级。更高的精度已经没有意义,因为\(f(x)\)和\(f(\bar{x})\)将不能分辨。所以“it is hopeless to ask for a bracketing interval of width less than \(\epsilon\) times its central value.”
抛物线插值和Brent算法
黄金分割法的优点是稳定可靠,缺点是收敛速度很慢。为了加快收敛速度,可以考虑采用抛物线插值的方法在计算极值。其思路是利用光滑函数在极值点的抛物线性质,在极值点附近用3个点确定的抛物线来近似函数\(f(x)\),然后再求抛物线的极值点\(\tilde{x}\)。用\(\tilde{x}\)替代上文的黄金分割点进行迭代更新。这种方法在很靠近极值点的时候相对黄金分割法有更快的收敛速度,但是在远离极值点也有可能不工作。所以,可以将黄金分割法和抛物线插值法综合起来应用,在抛物线插值法可以成功的时候,用该方法进行一次迭代,如果抛物线插值法不适用,则退回到黄金分割法。同样地,可以类比求根算法中二分法和牛顿法的结合算法。Brent算法就是一种结合两种方法的优化算法,具体可以参见文献[1]。
多维优化方法
概述
如前所述,多维优化问题可以分解成为在多个方向上的一维优化问题。具体来说,多维优化的每一步需要解决两个问题:
(a). 确定搜索方向\(\mathbf{p}\)
(b). 确定沿搜索方向前进的步长\(\delta\mathbf{x}\)
解决以上问题一般有两个途径:直线搜索(line search)法和置信域(trust region)法。其中直线搜索方法解决问题的思路是依次解决问题(a)和(b),即先确定搜索方向\(\mathbf{p}\),然后在该搜索方向上执行一维优化,确定\(\delta\mathbf{x}\)。而置信域方法则可以同时得到\(\mathbf{p}\)和\(\delta\mathbf{x}\)。
直线搜索方法
直线搜索方法的迭代公式可以写为:
\(\mathbf{x}_{k+1}=\mathbf{x}_k+\alpha_k\mathbf{p}_k\)
其中\(\mathbf{p}_k\)和\(\alpha_k\)分别表示搜索方向和步长。直线搜索方法的就是要确定\(\mathbf{p}_k\)和\(\alpha_k\),下面分别介绍。
步长选择
假设搜索方向\(\mathbf{p}_k\)已知,一个最直观的确定\(\alpha_k\)的思路就是令\(\alpha_k\)为函数\(\phi(\alpha)=f(\mathbf{x}+\alpha\mathbf{p}_k)\)的极小值(全局的或局部的)。这是个一维最优化问题,可以利用上一章介绍的方法进行求解。这种方法称作“精确直线搜索”,因为它试图计算的是沿着搜索方向上目标函数精确的极小值。但是,在每一步都求解一个一维优化问题的代价通常会比较大,所以通常采用的是一种称作“非精确直线搜索”的方法。该方法在每一步中寻找一个不是“特别坏”的步长,相对精确直线搜索的计算量更少。
Wolfe条件
Wolfe条件是一种用于非精确直线搜索中的判据,它包含两部分:
\(f(\mathbf{x}_k+\alpha_k\mathbf{p}_k)\leq f(\mathbf{x}_k)+c_1\alpha_k\nabla f(\mathbf{x}_k)^T\mathbf{p}_k\)
\(\nabla f(\mathbf{x}_k+\alpha_k\mathbf{p}_k)^T\mathbf{p}_k\ge c_2\nabla f(\mathbf{x}_k)^T\mathbf{p}_k\)
其中\(0<c_1<c_2<1\)。上式的第一个称作充分下降(sufficient decrease)条件,第二个条件称作曲率条件(curvature condition)。
在充分下降条件中,\(\alpha_k\)的取值保证\(\phi(\alpha)\)低于其与纵轴的交点\(\phi(0)=f(\mathbf{x}_k)\)。只要满足这个条件的\(\alpha_k\)均被认为是合适的。而曲率条件是要求\(\alpha_k\)的取值使\(\phi(\alpha)\)不能出现剧烈的下降,只能是较为平坦的下降、不下降或上升。因为剧烈的下降意味着\(\alpha_k\)的值还可以再更新以得到更合适的值,平坦的下降或者不下降则说明对\(\alpha_k\)进行更新的效果不明显,所以无需更新。注意曲率条件并不要求\(\phi(\alpha)\)的取值小于\(\phi(0)\),如果没有充分下降条件约束,仅凭曲率条件完全可以选择出\(\alpha_k\)使\(\phi(\alpha)\)比原始值更大。
根据上述Wolfe条件中的曲率条件可知,选取\(\alpha_k\)使\(\phi\'(\alpha)\)为一个很大的正数也是被允许的,然而这意味着在该点附近函数值还会发生很大的改变,即该点仍然距离极值点或驻点较远。为了强制要求直线搜索得到的下一步迭代点在极值点或驻点附近,可以采用所谓的强Wolfe条件:
\(f(\mathbf{x}_k+\alpha_k\mathbf{p}_k)\leq f(\mathbf{x}_k)+c_1\alpha_k\nabla f(\mathbf{x}_k)^T\mathbf{p}_k\)
\(|\nabla f(\mathbf{x}_k+\alpha_k\mathbf{p}_k)^T\mathbf{p}_k|\le c_2|\nabla f(\mathbf{x}_k)^T\mathbf{p}_k|\)
可见强Wolfe条件相对Wolfe条件的改变仅仅是要求\(\phi\'(\alpha)\)不能取太大的正数。
可以证明对于所有的连续可微的有界函数,都能找到满足Wolfe条件的步长。证明可以参考文献[3]的引理3.1。
关于采用Wolfe条件的非精确直线搜索算法的实现,可以参考文献[3]的3.5节。
回溯
在Wolfe条件中,充分下降条件总可以靠设置足够小的步长\(\alpha_k\)得到满足,所以只有充分下降条件一般是不能选出较好的直线搜索步长的。但是可以利用一种叫回溯(backtracing)的方法,使得仅仅依靠充分下降条件也能选出合适的步长。回溯算法的伪代码如下:
choose \(\bar{\alpha}>0\),\(\rho\in (0,1)\),\(c\in (0,1)\);
set \(\alpha=\bar{\alpha}\);
repeat until \(f(\mathbf{x}_k+\alpha\mathbf{p}_k)\leq f(\mathbf{x}_k)+c\alpha\nabla f(\mathbf{x}_k)^T\mathbf{p}_k\)
set \(\alpha=\rho\alpha\);
end
\(\alpha_k=\alpha\)
回溯算法是将步长从一个起始值\(\bar{\alpha}\)开始,按照比例\(\rho\)不断缩小,直到满足充分下降条件时即可停止,因此可以避免步长选择过小的问题。关于参数\(\bar{\alpha}>0\),\(\rho\in (0,1)\),\(c\in (0,1)\)的选取,文献[2]给出的范围是\(\bar{\alpha}=1\),\(\rho\in (0.1,0.8)\),\(c\in (0.01,0.3)\)。
回溯算法比较适合于牛顿法,但是对于拟牛顿法和共轭梯度法表现不好[3]。
收敛性分析
文献[3]中比较详细地讨论了Wolfe条件下非精确直线搜索的收敛性。其中Zoutendijk定理说明,只要保证搜索方向不是垂直于当前梯度方向,那么梯度会随着迭代而最终趋近于0,即迭代会收敛到驻点。对于梯度下降法,搜索方向沿着负梯度方向,所以总是收敛的。
收敛速率
[2-3]都讨论了\(f(\mathbf{x})\)为强凸函数(\(\nabla^2f-m\mathbf{I}\)为正定矩阵)情况下采用精确直线搜索梯度下降方法的收敛性。收敛速度都是线性的,而且收敛速度与\(f(\mathbf{x})\)在最优解处的黑塞矩阵的条件数(或者\(f(\mathbf{x})\)下水平集的条件数,因为\(f(\mathbf{x})\)的下水平集在极值点处可以用黑塞矩阵表示的椭球做近似)关系很大。黑塞矩阵条件数越接近1,收敛速度越快。如果条件数远大于1,梯度下降法的收敛速度会变得非常慢。而牛顿法和拟牛顿法的思路就是通过线性变换使黑塞矩阵条件数近似等于1,据此得到更快的收敛速度。
搜索方向选择
最速下降法
对目标函数\(f(\mathbf{x})\)在\(\mathbf{x}\)处做一阶泰勒展开,有:
\(f(\mathbf{x}+\mathbf{v})\approx f(\mathbf{x})+\nabla f(\mathbf{x})^T\mathbf{v}\)
其中上式右端第二项代表目标函数在\(\mathbf{x}\)沿方向\(\mathbf{v}\)的方向导数。如果该方向导数为负数,则为下降方向。
所谓的最速下降法,就是寻找\(\mathbf{v}\)使目标函数在沿着该方向的方向导数最小(注意下降方法要求方向导数为负值,方向导数最小即绝对值最大)。但是显然可以发现,只要\(\mathbf{v}\)是下降方向,增加\(\mathbf{v}\)的长度便可以使方向导数任意小,所以为了让问题有意义,需要限制\(\mathbf{v}\)的大小在一定范围内。
根据以上思路,定义规范化的最速下降方向如下:
\(\mathbf{v}=\textrm{argmin}\{\nabla f(\mathbf{x})^T\mathbf{v}\},\|\mathbf{v}\|=1\)
其中\(\|\|\)表示\(\mathbf{R}^{\mathbf{n}}\)上的任意范数。
当\(\|\|\)取做1-范数时,对应的最速下降法称作“坐标下降法”,此时\(\mathbf{v}\)的方向沿着梯度\(\nabla f(\mathbf{x})\)绝对值最大的分量;当\(\|\|\)取做欧几里得范数(即单位矩阵\(I\)定义的2-范数)时,就得到所谓的“梯度下降法”,因为显然\(\mathbf{v}\)沿着负梯度方向的方向导数最小;当\(\|\|\)取做\(f(\mathbf{x})\)的二阶导数\(\nabla^2 f(\mathbf{x})=\mathbf{H}\)定义的2-范数时,得到的就是牛顿法。而各种拟牛顿法也可以认为是取\(\|\|\)为\(2-\mathbf{P}^{-1}\)范数的最速下降法,其中\(\mathbf{P}\)为对称正定(symmetric positive definite, SPD)矩阵,通常用来模拟真实的黑塞矩阵\(\mathbf{H}\)
坐标下降法
坐标下降法的下降方向为:
\(\mathbf{v}=\textrm{argmin}\{\nabla f(\mathbf{x})^T\mathbf{v}\},\|\mathbf{v}\|_1=1\)
令\(i\)满足\(\|\nabla f(\mathbf{x})\|_{\infty}=|(\nabla f(\mathbf{x}))_i|\),即\(i\)代表梯度向量中最大的那个分量,则易知下降方向为:
\(\mathbf{v}=-\textrm{sign}[(\nabla f(\mathbf{x}))_i]\mathbf{e}_i\)
对于上式一个直观的解释如下:由于\(\|\mathbf{v}\|_1=1\),即\(\mathbf{v}\)的坐标分量绝对值之和为1,可以把\(\nabla f(\mathbf{x})^T\mathbf{v}\)看做是对梯度向量\(\nabla f(\mathbf{x})^T\)的各个分量以\(\mathbf{v}\)为权重进行加权求和。所以只有当权重都集中在梯度向量最大的分量上才会得到最大值(对应负数的最小值)。
坐标下降法每次只对变量\(\mathbf{x}\)的一个分量进行更改,该算法可以极大地简化或者省略步长选择步骤。
梯度下降法
当\(\|\|\)取做欧几里得范数时,可见\(\mathbf{v}\)就是沿着负梯度方向,所以该方法称作梯度下降法。在采用精确直线搜索方法的情况下,梯度下降发每一步搜索的方向都是沿着当前位置的负梯度方向\(\nabla f(\mathbf{x}_{k})\),而搜索步长的选择都使得在该方向上目标函数达到最小值,即沿着搜索方向上的方向导数为0,有\(\nabla f(\mathbf{x}_{k+1})^T\mathbf{v}_k=0\)。所以下一步的搜索方向\(\mathbf{v}_{k+1}=\nabla f(\mathbf{x}_{k+1})\)与当前搜索方向\(\mathbf{v}_k\)垂直。在二维情况下,采用精确直线搜索的梯度下降法的搜索路径会呈现出台阶形状,参考[2]的图9.2和[3]的图3.7。当采取非精确直线搜索时,连续两次搜索方向不一定垂直,但还是会呈现出锯齿状。文献[2]的图9.3给出了回溯直线搜索法的搜索路径。出现这种锯齿状路径的原因是每次搜索方向不能保证互相“独立”,沿着\(\mathbf{v}_{k+1}\)方向上的极小值不一定是(基本上都不会是)沿着\(\mathbf{v}_{k}\)方向上的极小值,所以下一次搜索时又要重复进行计算。
梯度下降法因为只用到了目标函数的一阶导数,涉及到的计算较为简单,而且梯度下降方法选择的搜索方向可以保证是收敛的。其缺点就是收敛速度非常的慢,尤其是目标函数下水平集的条件数很大的时候。
更新:梯度下降的一个解释可以参见MLAPP 13.4节。
牛顿法
从收敛性分析可知,最速下降法的收敛速度与待优化函数\(f(\mathbf{x})\)在极值点的黑塞矩阵\(\mathbf{H}\)条件数有很大关系,\(\mathbf{H}\)越接近单位矩阵\(\mathbf{I}\),收敛性越好。牛顿法相比梯度下降法,采用\(2-\mathbf{H}\)范数来替代梯度下降法中的欧几里得范数。这相当于对原始变量\(\mathbf{x}\)做了个线性变换\(\bar{\mathbf{x}}=\mathbf{H}^{\frac{1}{2}}\mathbf{x}\)。在此变换下,可以发现黑塞矩阵\(\bar{\mathbf{H}}=\mathbf{H}^{-1}\mathbf{H}=\mathbf{I}\),因此可以期望牛顿法具有更好的收敛性(当距离极值点较近的时候,牛顿法具有二次收敛速度)。
牛顿法的另一个解释是在极小值附近用二次曲面来近似目标函数\(f(\mathbf{x})\),然后求二次曲面的极小值作为步径。为了看出这一点,考虑\(f(\mathbf{x})\)的二阶泰勒展开:
\(f(\mathbf{x}+\mathbf{v})\approx f(\mathbf{x})+\nabla f(\mathbf{x})^T\mathbf{v}+\frac{1}{2}\mathbf{v}^T\nabla^2 f(\mathbf{x})\mathbf{v}\)
上式是关于\(\mathbf{v}\)的二次函数,易知其极小值为\(\mathbf{v}=-[\nabla^2 f(\mathbf{x})]^{-1}\nabla f(\mathbf{x})\),刚好是牛顿法的步径。相比之下,梯度下降法则是用平面来近似(参见最速下降法的推导)。
牛顿法相比梯度下降法的优点有:
- 在极值点附近具有二次收敛性。一旦进入二次收敛阶段,最多再经过6次的迭代便可以得到很高精度的解。
- 具有仿射不变性,对坐标选择或者目标函数下水平集不敏感(牛顿法本身相当于做了一次线性变换以保证目标函数下水平集近似为\(\mathbf{I}\))。
- 和问题规模有很好的比例关系,求解\(\mathbf{R}^{1000}\)中问题的性能和\(\mathbf{R}^{10}\)中问题的性能类似。
- 性能不依赖于参数的选择。
牛顿法的缺点是需要计算黑塞矩阵和其逆矩阵,这使得存储和计算的开销相比梯度下降法要大很多。为了解决这个问题,人们提出了所谓的拟牛顿法。
另外,牛顿法还有个问题就是在远离极值点的时候,目标函数的黑塞矩阵可能并不是正定的,所以牛顿法给出的搜索方向并不能保证是下降方向。为了解决这个问题,可以对黑塞矩阵做一些更改,比如加上一个矩阵\(\mathbf{E}\)使其变得正定。有关这方面的详细讨论可参见文献[3]。牛顿法的C++代码可以参考文献[2]。
拟牛顿法
拟牛顿法与牛顿法类似,只是用一个SPD矩阵代替精确的黑塞矩阵。
共轭方向法
一组矢量\(\{\mathbf{v}_1,\mathbf{v}_2,\cdots,\mathbf{v}_n\}\)关于SPD矩阵\(\mathbf{A}\)共轭是指对于任意的\(i\neq j\),都有:
\(\mathbf{v}_i^T\mathbf{A}\mathbf{v}_j=0\)
由于\(\mathbf{A}\)是SPD矩阵,所以上式的意思实际上就是\(\mathbf{v}_i\)和\(\mathbf{v}_j\)关于\(\mathbf{A}\)定义的内积是正交的。
共轭方向法就是在问题空间\(\mathbf{R}^n\)中选择一组互相共轭的基向量作为搜索方向。相比于最速下降法,它的一个显著优点是每次搜索是“互不干扰”的,所以最多需要\(N\)次便可以保证求出目标函数的一个极小值,其中\(N\)是问题的维度。为了看出这一点,仍然考虑目标函数的二阶泰勒展开式:
\(f(\mathbf{x}+\mathbf{v})\approx f(\mathbf{x})+\nabla f(\mathbf{x})^T\mathbf{v}+\frac{1}{2}\mathbf{v}^T\nabla^2 f(\mathbf{x})\mathbf{v}\)
令\(\mathbf{b}=-\nabla f(\mathbf{x})\),\(\mathbf{A}=\nabla^2 f(\mathbf{x})\),则有:
\(f(\mathbf{x}+\mathbf{v})\approx f(\mathbf{x})-\mathbf{b}^T\mathbf{v}+\frac{1}{2}\mathbf{v}^T\mathbf{A}\mathbf{v}\)
易知\(f(\mathbf{v})\)的梯度可以近似等于\(\mathbf{Ax-b}\)。另外值得一提的是,求解线性方程组\(\mathbf{Ax=b}\)可以认为是关于二次型\(f(\mathbf{x})=c-\mathbf{b}^T\mathbf{x}+\frac{1}{2}\mathbf{x}^T\mathbf{A}\mathbf{x}\)的优化问题,即寻找令\(f(\mathbf{x})\)梯度为0的\(\mathbf{x}\)。
当\(\mathbf{x}\)沿某一方向移动时,梯度的变化可以表示为\(\delta\nabla f=\mathbf{A}(\delta\mathbf{x})\)。假设已经沿着方向\(\mathbf{v}_i\)达到了极小值,现在要沿着\(\mathbf{v}_j\)进行搜索。可以发现沿\(\mathbf{v}_j\)的搜索并不会破坏之前沿\(\mathbf{v}_i\)的搜索结果,因为根据共轭性可知,沿\(\mathbf{v}_j\)搜索造成的梯度的变化\(\mathbf{A}\mathbf{v}_j\)垂直于之前的搜索方向\(\mathbf{v}_i\),即\(\mathbf{v}_i^T\mathbf{A}\mathbf{v}_j=0\)。所以更新后的梯度仍然垂直于\(\mathbf{v}_i\),也就是说沿着\(\mathbf{v}_i\)的方向导数仍然为0,新的搜索不会破坏沿之前方向上已经找出的极小值。所以沿着一组共轭方向集依次做搜索时,不需要做重复搜索。文献[3]给出了更正式的证明,每次搜索得到的\(\mathbf{x}_k\)都是\(f(\mathbf{x})\)在子空间\(\mathbf{x}_0+\textrm{span}\{\mathbf{v}_1,\mathbf{v}_2,\cdots,\mathbf{v}_{k-1}\}\)上的极小值,其中\(\mathbf{x}_0\)是初始试探解。
共轭方向法也可以从几何上比较直观地进行理解。考虑关于\(\mathbf{v}\)的二次函数\(\varphi(\mathbf{v})=a+\mathbf{b}^T\mathbf{v}+\frac{1}{2}\mathbf{v}^T\mathbf{A}\mathbf{v}\),如果\(\mathbf{A}\)在标准正交基\(E=\{\mathbf{e}_1,\mathbf{e}_2,\cdots,\mathbf{e}_n\}\)下是对角阵,那么\(\varphi(\mathbf{v})\)对应的图形就是主轴沿笛卡尔坐标系坐标轴的椭球。此时显然可以直接沿着各坐标轴的方向逐次求极小值,而且在每个坐标轴方向上求极小值都不会对其他坐标产生影响,最多需要遍历\(E\)便可求出\(\varphi(\mathbf{v})\)的极小值。
如果\(\mathbf{A}\)不是对角的,\(\varphi(\mathbf{v})\)所代表的椭球的主轴方向将发生倾斜,不再沿着坐标轴,这时再沿着坐标轴的方向逐次求极小值的时候,就会发现每次新的方向求极小都可能影响到之前已经搜索得到的极小值点。为了解决这个问题,可以考虑把该组基做线性变换,得到新的基\(P=\{\mathbf{p}_1,\mathbf{p}_2,\cdots,\mathbf{p}_n\}\)。在基\(P\)下,\(\mathbf{A}\)成为对角的,椭球的主轴也沿着\(P\)的方向,这样就回到了之前的情况。等价地在原空间中看,现在新的搜索方向就是沿着一组关于\(\mathbf{A}\)共轭的基向量\(P\)来进行。
共轭梯度法
共轭梯度法(Conjucate Gradient,CG)是一类特殊的共轭方向法,它的第一个共轭矢量\(\mathbf{p}_0\)选择为初始点的负梯度方向,接下来的共轭矢量则是前一个共轭矢量和当前点负梯度方向的线性组合的形式。
共轭梯度法的一个特点是在计算当前共轭矢量的时候只需要之前的一个矢量即可。也就是说,计算\(\mathbf{p}_k\)时只需要\(\mathbf{p}_{k-1}\),计算得到的\(\mathbf{p}_k\)自动与之前的共轭方向\(\mathbf{p}_1,\mathbf{p}_2,\cdots,\mathbf{p}_{k-2}\)正交。它的另一个特点是并不需要计算目标函数\(f(\mathbf{x})\)的二阶导数,整个迭代过程只需要计算目标函数的一阶导数以及直线搜索。
文献[3]对共轭梯度法的理论基础做了比较详细的介绍,包括基本原理,收敛性,算法伪代码等;文献[1]从实现的角度进行了介绍,并给出了代码。 关于共轭梯度法和拟牛顿法的优劣,文献[1]中说道**On the other hand, there is not, as far as we know, any overwhelming advantage that the quasi-Newton methods hold over the conjugate gradient techniques, except perhaps a historical one. **
此外,在矩阵计算里面,共轭梯度法可以归结到克雷洛夫子空间(Krylov Subspace)方法中去。作为克雷洛夫子空间方法中的一个判据,共轭梯度法寻找的是最小化\(\|\mathbf{r}_k\|_{A^{-1}}\)的解\(\mathbf{x}_k\)。其中\(\mathbf{r}_k=A\mathbf{x}_k-b\)是上文中的梯度,同时也是近似解\(\mathbf{x}_k\)对应的残量,\(\|\dot\|_{A^{-1}}\)表示矩阵\(A^{-1}\)定义的2-范数。其他的判据还有包括MINRES、GMRES、SYMMLQ等等。
置信域法
置信域法的思路是在一个有限区域内把目标函数近似地用一个二次模型函数\(m_k(\mathbf{p})=f(\mathbf{x}_k)+\nabla f(\mathbf{x}_k)^T\mathbf{p}+\frac{1}{2}\mathbf{p}^T\mathbf{B}_k\mathbf{p}\)替代,在该区域内寻找\(m_k(\mathbf{p})\)的极小值\(\mathbf{p}\),并令\(\mathbf{x}_{k+1}=\mathbf{x}_{k}+\mathbf{p}\)。如果\(m_k\)中的\(\mathbf{B}_k\)是目标函数的二阶导数\(\nabla^2 f(\mathbf{x}_k)\),就称作置信域-牛顿法,可见此时置信域法和直线搜索牛顿法的几何解释是一致的。