标签归档:共轭梯度

Notes on Conjugate Gradient Method

版权声明:原创作品,欢迎转载,但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明,否则将追究法律责任。

这是一个关于共轭梯度法的笔记。请大家注意的是,这是个笔记,并不是一个教程,因此少不了跳跃和欠解释的地方。对CG方法了解不多的同学请移步这里

线性方程组和极小化问题

一个关于对称矩阵A的线性方程组Ax=b等价于求解如下极小值问题:

    \[ f(x) = \frac{1}{2} \transp{x} A x - \transp{b} x \]

这很容易说明,我们微分目标函数得

(1)   \begin{eqnarray*} \td f &=& \frac{1}{2}(\td\transp{x}Ax+\transp{x}A\td x)-\transp{b} \td x\\ &=& \transp{(Ax - b)} \td x \end{eqnarray*}

所以\td f=0意味着Ax=b.

x^*为问题的解,e为偏离极小值点的位移,即

(2)   \begin{eqnarray*} Ax^*&=&b\\ e&=&\Delta x=x-x^* \end{eqnarray*}

我们定义残差或负梯度r=b-Ax,容易看出,r正比于偏离梯度零点x^*的位移e:

    \[ r=-Ae \]

最速下降法和共轭梯度法简述

最速下降法:搜索方向为本轮迭代初始点的梯度方向,搜索到梯度与搜索方向正交的位置开始下一轮迭代。

共轭梯度法:搜索方向为本轮迭代初始点的梯度方向(残差方向)用A-内积做正交化,扣除掉之前所有搜索方向的分量给出,搜索到梯度与搜索方向正交的位置开始下一轮迭代。

共轭梯度法的理解

首先考虑一个简单情形,即A=I时。此时f退化,且x的不同分量解耦,f的极小值问题分解为各分量上的极小值问题。我们只需要在各分量方向上找到极小值,问题就得解了。具体的说,我们从任意一个初始点出发,在方向d_0上寻找极小值,然后从这个极小值出发继续在d_1寻找d_1方向的极小,以此类推最后到达全局极小值。

A为正定对称矩阵的一般情形,其实只是上述情形中的x表述在一个新坐标系下的结果(下文详细解释)。为了求解这种情形,我们考察上述方法的求解轨迹在新坐标系里的对应。

首先,A=I时的各个正交搜索方向在新坐标系里对应着一组A-正交的方向,因此我们需要设法找出这样一组正交方向,然后各自求解这些方向上的极小值。

求解特定方向上极小值的方法不变,极小点就是这个方向上梯度正交于搜索方向时的点(原因很简单,这说明搜索方向上梯度分量为零)。

最后一点,如何找到这样一组A-正交方向?

一般方法是找一组现成的完备基,进行A-正交化。问题是计算量比较大。

记迭代到第i步时,所搜索过的这样一组A-正交的搜索方向单位矢量为\{d_0,d_1,...,d_{i-1}\},它们撑起的子空间

(3)   \begin{eqnarray*} D_i=\mathrm{span}\{d_0,d_1,...,d_{i-1}\} \end{eqnarray*}

即第i步前已经搜索过的子空间。

共轭梯度法的一个关键之处在于它可以从每次迭代的r_i中快速构造出新的搜索方向d_i,而不需要对r_i进行完整的正交化手续(即不需依次扣除r_id_0,d_1d_{i-1}方向的分量)。原因在于r_i自然的和D_{i-1}A-正交,因此构造d_i只需要从r_{i}中扣除掉d_{i-1}分量即可。

因此,共轭梯度算法有两处关键,一处为通过利用A-正交关系来把搜索问题转换到一个特殊坐标系,使得各搜索方向上的极小值问题解耦;一是利用r_iD_{i-1}是A-正交的这一规律来充分简化A-正交搜索方向的构造。

接下来我们详细论证这个算法。

第一个关键点,坐标变换和子问题解耦

(4)   \begin{eqnarray*} A&=&\transp{T}T\\ y&=&T x\\ b^{\prime}&=&T^{-T} b\\ \end{eqnarray*}

其中T^{-T}表示矩阵T求逆后做转置,很拗口,不过这个不重要。

于是我们有

(5)   \begin{eqnarray*} f(x) &=& \frac{1}{2} \transp{y} y - b^{\prime T} y\\ &=& \frac{1}{2}\sum_{n}y_n^2 - b^{\prime}_n y_n \equiv \frac{1}{2}\sum_{n} P_n(y_n) \end{eqnarray*}

可见f已经被分解为一系列独立的子问题P_n(y_n)T正是将问题解耦所用的坐标变换。

在新的\{y_n\}坐标系中,整个极小值问题可以通过在一组正交方向上依次寻找极小点来解决。但求解T需要对A做分解,并非一个计算上经济的解决方法。这一个变换的真正意义在于,我们可以通过它来找到求解轨迹在\{x_n\}坐标系中的对应。

T是线性变换,因此\{y_n\}中的一组直线仍被变换为\{x_n\}中的一组直线。\{y_n\}坐标系中,我们依次沿着一组正交方向求解极小值,因此只需要找到这组方向在\{x_n\}坐标系中的对应(一组未必正交的方向),就可以等价的直接在这组方向之上寻找极小值,而不需要经过坐标变换。

考虑这组正交方向单位矢\{y_n\}

(6)   \begin{eqnarray*} &&\transp{y_i} y_j = \delta_{ij}\\ &\Leftrightarrow& \transp{(T x_i)} T x_j=\delta_{ij}\\ &\Leftrightarrow& \transp{x_i} \transp{T}T x_j = \delta_{ij}\\ &\Leftrightarrow& \transp{x_i} A x_j = \delta_{ij} \end{eqnarray*}

可见,\{y_n\}这一组正交的方向在变换到\{x_n\}坐标系中后,虽然不是\{x_n\}中通常意义的正交矢量组,却可以在A-内积\langle \bullet,\bullet \rangle_A意义下成为一组A-正交矢量。这里的A-内积,指的是以对称矩阵A定义的内积运算

(7)   \begin{eqnarray*} \langle a,b \rangle_A \equiv \transp{a} A b \end{eqnarray*}

需要注意的是,上述结论是可逆的,如果\langle x_i,x_j \rangle_A = 0,同样也可以得出,x_i,x_j对应着\{y_n\}坐标系中两个正交的方向。因此,可以直接在\{x_n\}坐标系下,寻找一组A-正交的方向,然后在这一组方向上依次极小化目标函数,就可以找到整个问题的极小值。

第二个关键点,D_i子空间结构和正交关系

算法的更新规则为

(8)   \begin{eqnarray*} e_{i+1}&=&e_i + \alpha_i d_i\\ r_{i+1}&=&-A(e_i+\alpha_i d_i)=r_i-\alpha_i A d_i\\ d_{i}&=&r_i + \sum_{j<i} \beta_{ij} d_{j}\\ \beta_{ij}&=& - \frac{\transp{r_i}d_j}{\transp{d_j}Ad_j}\\ &=& - \frac{\langle r_i,d_j \rangle }{\langle d_j,d_j \rangle _A} \end{eqnarray*}

其中\alpha_i为第i步搜索方向上的位移长度,它的大小可以通过r_{i+1}应该和d_i正交这一要求定出

(9)   \begin{eqnarray*} &&\transp{r_{i+1}}d_i=0\\ &\Rightarrow&\transp{(e_i+\alpha_i d_i)} A d_i =0\\ &\Rightarrow&\alpha_i = - \frac{\langle d_i,e_i \rangle_A }{\langle d_i,d_i \rangle_A}= \frac{\langle d_i,r_i \rangle}{\langle d_i,d_i \rangle_A} \end{eqnarray*}

算法将参数空间划分成了一个层级结构,下面我们来看看这个结构的性质。

根据上述更新规则,可以发现D_i可以由另外几个矢量组撑起。首先,由d_i的更新规则可以发现,d_ir_id_{j<i}线性组合而成,另一方面,我们有初始条件d_0=r_0,归纳而知,d_i可以由\{r_{j\leq i}\}线性组合出来。于是我们有

(10)   \begin{eqnarray*} D_i=\mathrm{span}\{r_0,r_1,...,r_{i-1}\} \end{eqnarray*}

接着我们考虑r_i的更新规则。可以看到,D_{i+1}=\mathrm{span}\{D_{i},Ad_i\}.于是我们可以从D_1=\{d_0\}D_1=\{r_0\}构造出D_i

(11)   \begin{eqnarray*} D_i&=&\mathrm{span}\{d_0,Ad_0,A^2d_0,...,A^{i-1}d_0\} = \mathrm{span}\{d_0,A D_{i-1}\}\\ &=&\mathrm{span}\{r_0,Ar_0,A^2r_0,...,A^{i-1}r_0\} = \mathrm{span}\{r_0,A D_{i-1}\} \end{eqnarray*}

接下来我们再考虑各子空间和向量之间的正交关系。首先,e_iD_i的补空间里。这点是明显的,因为每一个优化步都会优化掉e_j在相应d_j方向上的分量,因而剩余的e_i就只位于D_i的补空间里。于是我们有

(12)   \begin{eqnarray*} e_i = \sum_{j \geq i} \delta_j d_j \end{eqnarray*}

这样ed之间的A-正交关系就很明显了。由于\langle d_i,d_j \rangle_A = \delta_{ij},而e_i只包含d_j>i分量,因而

(13)   \begin{eqnarray*} \langle e_i,d_{j<i} \rangle_A =0 \end{eqnarray*}

e_i和子空间D_i也是A-正交的。这意味着r_i正交于D_i

(14)   \begin{eqnarray*} &&\transp{e_i} A d_{j<i}=0\\ &\Rightarrow&\transp{(A e_i)} d_{j<i}=0\\ &\Rightarrow&\langle r_i,d_{j<i}\rangle =0 \end{eqnarray*}

到这里就可以解释为什么我们有\langle r_i,d_{j<i-1} \rangle_A =0了:

(15)   \begin{eqnarray*} &&\langle r_i,D_i \rangle =0\\ &\Rightarrow&\langle r_i,\mathrm{span}\{d_0,A D_{i-1}\} \rangle =0\\ &\Rightarrow&\langle r_i,A D_{i-1} \rangle =0\\ &\Rightarrow&\langle r_i,D_{i-1} \rangle_A =0 \end{eqnarray*}

至此我们完成了第二个关键点的论证,且更清楚的了解了参数空间,尤其是D_i子空间的结构。