共轭梯度法及其在LR中的应用（一）

共轭方向

无约束最优化方法的核心问题就是选择搜索方向

定义：设$A_{n\times n}$为对称正定矩阵，若$R^{n}$中的两个方向$d^{(1)}$和$d^{(2)}$满足$d^{(1)T}Ad^{(2)}=0$，则称这两个方向关于$A$共轭，或称它们关于$A$正交。

若$d^{(1)}$,$d^{(2)}$,...,$d^{(k)}$为$R^{n}$中k个方向关于$A$共轭，满足$d^{(i)T}Ad^{(j)}=0,i\neq j(i,j=1,...,k)$，也即$d^{(i)}$与$Ad^{(j)}$正交，称这组方向是A共轭的。

若A为单位阵，则这两个方向关于A共轭等价于两个方向正交。

二次函数  $f(x)=\frac{1}{2}(x-\bar{x})^{T}A(x-\bar{x})$ ，$\bar{x}$为一个定点，

$f(x)$的等值面$\frac{1}{2}(x-\bar{x})^{T}A(x-\bar{x})=c$是以$\bar{x}$为中心的椭球面。

由于$\bigtriangledown f(\bar{x})=A(\bar{x}-\bar{x})=0$，A正定，所以$\bar{x}$为$f(x)$的极小点。

设$x^{(1)}$为等值面上一点，该等值面在点$x^{(1)}$处的法向量$\bigtriangledown f(x^{(1)})=A(x^{(1)}-\bar{x})$ （注：二元函数一条等值线上一点的法向量即梯度）

设$d^{(1)}$为等值面在$x^{(1)}$处的一个切向量，显然$d^{(1)}$与$\bigtriangledown f(x^{(1)})$正交。同时，令$d^{(2)}=\bar{x}-x^{(1)}$，所以

$d^{(1)T}Ad^{(2)}=0$

结论：等值面上一点处的切向量与由这一点指向极小点的向量关于A共轭。

所以，极小化这个二次函数，若一次沿着 $d^{(1)}$和$d^{(2)}$进行一维搜索，必将到达最小点。

 

---

所谓最速真的是最速吗？

回顾最速下降法，用最速下降法极小化目标函数时，相邻两个搜索方向是正交的。

因为梯度是函数局部性质，从局部看，最速下降方向确实是函数值下降最快的方向，然而从整体看，由于锯齿现象的影响，则走过了许多弯路，使收敛速度大为减慢。最速下降法并不是收敛最快的方向。

对于二次型，从初始点$x_{(0)}$出发，最速下降法可由如下描述

下图中极小点位于$x_{(*)}=(0,0)^{T}$，令$e_{(i)}=x_{(*)}-x_{(i)}$为每次迭代的误差向量，显然这个向量是无法计算的因为我们并不知道$x_{(*)}$。

但是我们可以计算

$Ae_{(i)}=A(x_{(*)}-x_{(i)})=Ax_{(*)}-Ax_{(i)}=b-Ax_{(i)}$

记残差向量$r_{(i)}=b-Ax_{(i)}=Ae_{(i)}$，可以发现$r_{(i)}=-\bigtriangledown f(x_{(i)}$。

最速下降法即通过贪婪地减少误差项，在下一个点$x_{(i+1)}$时$r_{(i+1)}\perp r_{(i)}$(can only make square turns)。这个性质对于对称正定矩阵Hessian矩阵A的（二维）二次型来说，就意味着搜索方向一定会大于2次。为了确保迭代次数少于两次，就需要选择一个起始点$x_{(0)}$，$\bigtriangledown f(x_{(0)})\parallel e_{(0)}$或者$\bigtriangledown f(x_{(0)})\perp e_{(0)}$。当然我们无法计算误差向量。

共轭梯度法的基本思想：利用已知点处的梯度构造一组共轭方向，并沿这组方向进行搜索，求出目标函数值得极小点。根据共轭方向的基本性质，这种方法具有二次终止性。

对于二次型的Hessian矩阵$A$为常值，所以可以构造一组搜索方向来避免redundant caculations。

假设有一组n个独立线性无关的方向向量$d_{(0)},....,d_{(n-1)}$以及对应的尺度$\alpha_{(0)},....,\alpha_{(n-1)}$

$e_{(0)}=\sum_{j=0}^{n-1}\alpha_{(j)}d_{(j)}$

我们要做的就是将误差$e_{(0)}$迭代减小为0.假定我们知道$e_{(0)}$并求解上式，则时间复杂度为$O(n^{3}$。但是，如果方向向量$d_{(0)},....,d_{(n-1)}$是相互正交的，则可以在$O(n)$求解得$\alpha_{(j)}$。但是问题是我们是不知道$e_{(0)}$的。

 又是先前的那个trick，我们计算$Ae_{(0)}$

$Ae_{(0)}=\sum_{j=0}^{n-1}\alpha_{(j)}Ad_{(j)}$

求解这个问题仍然是$O(n^{3}$，相比$d_{(0)},....,d_{(n-1)}$是相互正交的，我们让$d^{(i)T}Ad^{(j)}=0,i\neq j$。求解$\alpha_{(k)}$可通过左乘$d_{(k)}^{T}$，

$d_{(k)}^{T}Ae_{(0)}=\sum_{j=0}^{n-1}\alpha_{(j)}d_{(k)}^{T}Ad_{(j)}=\alpha_{(k)}d_{(k)}^{T}Ad_{(k)}$

那么我们可以求得

$\alpha_{(k)}=\frac{d_{(k)}^{T}Ae_{(0)}}{d_{(k)}^{T}Ad_{(k)}}=-\frac{d_{(k)}^{T}r_{(0)}}{d_{(k)}^{T}Ad_{(k)}}$

其中$Ad_{(k)}\neq 0$

共轭梯度法每一步迭代仅依赖于矩阵A和先前的迭代项，除了$\alpha_{(k)}$，实际上，$\alpha_{(k)}$也可以写成先前的迭代项形式，因此我们只需要存储一些向量和矩阵A。

当A_{n\times n}，CG方法可以在不超过n次迭代内得到精确解$x_{(*)}$。但是当在计算机上执行舍入运算时，这个性质无法保证。因为计算误差会减少方向向量之间的共轭性，从而使得算法可能会超过n次迭代。

CG的收敛率取决于

其中$\kappa=A/a$为条件数，Hesse矩阵$A$的最小特征值$a>0$，最大特征值为$A$。CG终止条件为$\lVert e_{i}\rVert_{A}/\lVert e_{i}\rVert_{A}<\epsilon$。这种最坏情况下CG的时间复杂度为$O(m\kappa$，$m$为A中的非零迹项。显然，最坏情况发生在A的特征值为均匀分布，而较好的收敛率则出现在有多重特征值时。

所以，选择迭代终止条件也是一个很有讲究的问题。

---

 用于二次函数的共轭梯度法

对于二次型$f:R^{n}\rightarrow R$

$f(x)=\frac{1}{2}x^{T}Ax-bx+c$

$A_{n \times n}$为对称正定矩阵，$x$和$b$为长度为$n$的向量，$c$为尺度。

假定$A$  SPD，即$x^{T}Ax\ge 0,\vee x\in R^{n}$

那么必有global minimum，最小点在point、line or plane上取得，不会有local minimum。

令$f^{'}(x)=Ax-b=0,f^{''}(x)=A$，则$Ax=b$

算法步骤：

1）选取初始点$x_{(0)}$，计算出目标函数在该点的梯度，若$\lVert g_{(0)}\rVert=0$，则停止计算；否则，令第一个搜索方向$d_{(0)}=-\bigtriangledown f(x_{(0)})=-g_{(0)}$，并记为此时方程的残差$r_{(0)}=b-Ax_{(0)}$

2)沿方向$d_{(0)}$搜索，得到点$x_{(1)}$，计算在点$x_{(1)}$处的梯度，若$\lVert g_{(1)}\rVert \neq0$，则利用$-g_{(1)}$和$d_{(0)}$构造第2个搜索方向$d_{(1)}$。先暂时假定我们已知该搜索方向。

一般地，若已知点$x_{(k-1)}$和搜索方向$d_{(k-1)}$，则从$x_{(k-1)}$出发，沿$d_{(k-1)}$进行搜索，得到

$x_{(k)}=x_{(k-1)}+\lambda_{(k-1)}(d_{(k-1)}$

其中步长$\lambda_{(k-1)}$满足$f(x_{(k-1)}+\lambda_{(k-1)}(d_{(k-1)})=min\quad x_{(k-1)}+\lambda d_{(k-1)}$。

令$\varphi(\lambda)=x_{(k-1)}+\lambda d_{(k-1)}$，求令$\varphi(\lambda)$的极小点。

令$\varphi^{'}(\lambda)=\bigtriangledown f(x_{(k)})^{T} d_{(k-1)}=0$即$\bigtriangledown f(x_{(k-1)}+\lambda d_{(k-1)})^{T} d_{(k-1)}$

根据二次函数的梯度表达式，可得

$(Ax_{(k)}-b)^{T}d_{(k-1)}=0$

$(A(x_{(k-1)}+\lambda_{(k-1)}d_{(k-1)})-b)^{T}d_{(k-1)}=0$

$(g_{(k-1)}+\lambda_{(k-1)}Ad_{(k-1)})^{T}d_{(k-1)}=0$

所以根据上式可得

$\lambda_{(k-1)}=-\frac{g_{(k-1)}d_{(k-1)}}{d_{(k-1)}^{T}Ad_{(k-1)}}$

3)计算$f(x)$在点$x_{(k)}$处的梯度，若$\lVert g_{(k)}\rVert \neq0$，则停止计算；否则，利用$-g_{(k)}$和$d_{(k-1)}$构造下一个搜索方向$d_{(k)}$，并使$d_{(k)}$和$d_{(k-1)}$关于$A$共轭

$d_{(k)}=-g_{k}+\beta_{(k)}d_{(k-1)}$

上式两端左乘$d_{(k-1)}^{T}A$，并令

$d_{(k-1)}^{T}Ad_{(k)}=-d_{(k-1)}^{T}Ag_{k}+\beta_{(k)}d_{(k-1)}^{T}Ad_{(k-1)}=0$

所以可得

$\beta_{(k)}=\frac{d_{(k-1)}^{T}Ag_{(k)}}{d_{(k-1)}^{T}Ad_{(k-1)}}$

综上分析，在第1个搜索方向取负梯度的前提下，伴随计算点的增加，构造出一组搜索方向，可以证明这组方向是关于A共轭的，因此，共轭梯度法具有二次终止性。（证明请参见《最优化方法》 陈宝林）

 

 

---

共轭梯度法：旨在求解线性方程Ax=b，A为正定矩阵。仅需一阶导数信息，但克服了最速下降法收敛慢的特点，又避免了牛顿法需要存储和计算Hessian矩阵并求逆的缺点。所以适合求解大型系统。当A is well-conditioned时，其远远快于其他方法如高斯消除法（时间复杂度$O(n^{3})$，显然不适用于求解大型线性方程)。通过迭代方法去寻找$x$的近似解。

Well-conditioned情况下：

Not Well-conditioned情况下：

 

参考资料：

1、陈宝林   《最优化理论与算法》  第二版

2、PaulKomarek       “LogisticRegressionforDataMiningand High-DimensionalClassification”

3、http://blog.163.com/yuyang_tech/blog/static/216050083201210294100593/

 

补充知识：
矩阵条件数的定义：

矩阵条件数的性质：

Matlab中计算矩阵条件数的函数：