共轭梯度法
最速下降法
考虑线性方程组 \(Ax=b\) ,其中 \(A\) 是对称正定阵,定义二次泛函
\[\varphi(x) = x^TAx - 2b^Tx \]定理 5.1.1 设 \(A\) 对称正定,则 \(Ax=b\) 的解等价于二次泛函 \(\varphi(x)\) 的极小值点.
证明 事实上,直接计算梯度有
\[\mathrm{grad}\ \varphi(x) = 2(Ax-b) = -2r \]则 \(Ax=b\) 的解是 \(\varphi(x)\) 的极小值点.
通常类似于盲人爬山法,给定任意一个初始向量 \(x_0\) ,确定一个下山方向 \(p_0\) ,沿着经过 \(x_0\) 方向为 \(p_0\) 的直线 \(x=x_0+\alpha p_0\) 找一个点 \(x_1 = x_0+\alpha_0 p_0\) ,使得该点是这条直线上 \(\varphi(x)\) 的最小值,以此类推,得到一系列 \(\alpha_k,\ p_k,\ x_k\) ;
我们先考虑如何确定步长 \(\alpha_k\) ,令
\[f(\alpha) = \varphi(x_k+\alpha p_k) = \alpha^2p_k^TAp_k - 2\alpha r_k^T p_k + \varphi(x_k) \]只需要找到导数为零的位置即可,则
\[\alpha_k = \dfrac{r_k^Tp_k}{p_k^TAp_k} \]这样就确定了步长,容易验证这恰好就是 \(f\) 的极小值点,因此我们得到 \(x_{k+1}\) ;
接着,由于我们知道梯度 \(\mathrm{grad}\ \varphi(x) = -2r\) ,梯度的负方向就是下降最快的方向,因此 \(p_k = r_k\) .
定理 5.1.2 设 \(A\) 的特征值为 \(0<\lambda_1\le \cdots\le \lambda_n\) ,则上述方法得到的序列 \(\{x_k\}\) 满足
\[\|x_{k} - x_*\|_A = \left(\dfrac{\lambda_n-\lambda_1}{\lambda_n+\lambda_1}\right)^k\|x_0-x_*\|_A \]其中 \(x_* = A^{-1}b,\ \|x\|_A = \sqrt{x^TAx}\) .
引理 5.1.1 设 \(A\) 的特征值为 \(0<\lambda_1\le \cdots\le \lambda_n\) , \(P(t)\) 是一个 \(t\) 的多项式,则
\[\|P(A)x\|_A \le \max_{1\le i\le n}|P(\lambda_i)|\|x\|_A,\quad x\in\mathbb{R}^n \]此引理说明了 \(A\) 的多项式与其特征值的多项式之间的关系.
共轭梯度法
最速下降法虽然保证了每一步都选取最优的下降方向和步长,但是并不一定是整体最优,因此我们需要进行改进:
给定 \(x_0\) 在第一步仍选取 \(p_0 = r_0\) ,则有
\[\alpha_0 = \dfrac{r_0^Tr_0}{p_0^TAp_0},\quad x_1 = x_0 + \alpha_0p_0,\quad r_1 = b-Ax_1 \]在之后第 \(k+1\) 步,我们在 \(r_k,\ p_{k-1}\) 张成的平面
\[\pi_2 = \{x=x_k+\xi r_k+\eta p_{k-1}: \xi,\eta\in\mathbb{R}\} \]内找到新的下山方向,我们可以考虑二元函数
\[\psi(\xi,\eta) = \varphi(x_k+\xi r_k+\eta p_{k-1}) \]为了找出极值点,可以令 \(\phi\) 的偏导为零
\[\begin{aligned} 0 = \dfrac{\partial \psi}{\partial \xi} &= 2\left(\xi r_k^TA r_k + \eta r_k^TA p_{k-1} - r_k^Tr_k\right)\\ 0 = \dfrac{\partial \psi}{\partial \eta} &= 2\left(\xi r_k^TA p_{k-1} + \eta p_{k-1}^TA p_{k-1}\right) \end{aligned} \]其中第二个式子用到了 \(r_k,\ p_{k-1}\) 正交,这是因为每当选定了下山方向 \(p_{k-1}\) ,都会按照最速下降法中的方法取
\[\alpha_{k-1} = \dfrac{r_{k-1}^Tp_{k-1}}{p_{k-1}^TAp_{k-1}},\quad x_k = x_{k-1} + \alpha_{k-1}p_{k-1},\quad r_k = b - Ax_k \]这样我们就有
\[\begin{aligned} r_k^Tp_{k-1} &= b^Tp_{k-1} - x_{k}^TAp_{k-1}\\ &= b^Tp_{k-1} - (x_{k-1}^T + \alpha_{k-1}p_{k-1}^T)Ap_{k-1}\\ &= r_{k-1}^Tp_{k-1} - \alpha_{k-1}p_{k-1}^TAp_{k-1}\\ &= 0 \end{aligned} \]设唯一极小值点
\[\widetilde{x} = x_k+\xi_0r_k + \eta_0p_{k-1} \]其中 \(\xi_0,\ \eta_0\) 满足
\[\left\{ \begin{aligned} &\xi_0r_k^TAr_k + \eta_0r_k^TAp_{k-1} = r_k^Tr_k\\ &\xi_0r_k^TAp_{k-1} + \eta_0p_{k-1}^TAp_{k-1} = 0 \end{aligned} \right. \]注意到上式蕴含着 \(r_k\neq0\) 必有 \(\xi_0\neq 0\) ,事实上,若 \(\xi_0=0\) ,则
\[\left\{ \begin{aligned} &\eta_0r_k^TAp_{k-1} = r_k^Tr_k\\ &\eta_0p_{k-1}^TAp_{k-1} = 0 \end{aligned} \right.\ \Rightarrow\ \eta_0 \neq 0,\ p_{k-1} = 0\ \Rightarrow\ r_k = 0 \]因此根据极值点位置,我们可以取
\[p_k = \dfrac{1}{\xi_0}(\widetilde{x}-x_k) = r_k+\dfrac{\eta_0}{\xi_0}p_{k-1} \]作为新的下山方向,令 \(\beta_{k-1} = \eta_0/\xi_0\) ,代入第二个方程得到
\[p_k = r_k+\beta_{k-1}p_{k-1},\quad \beta_{k-1} = -\dfrac{r_k^TAp_{k-1}}{p_{k-1}^TAp_{k-1}} \]注意这样确定的 \(p_k\) 满足 \(p_k^TAp_{k-1} = 0\) ,即 \(p_k,\ p_{k-1}\) 相互共轭,于是有
\[\begin{aligned} \alpha_k &= \dfrac{r_k^Tp_k}{p_k^TAp_k},\quad x_{k+1} = x_k + \alpha_kp_k,\quad r_{k+1} = b-Ax_{k+1}\\ \beta_k &= -\dfrac{r_{k+1}^TAp_k}{p_k^TAp_k},\quad p_{k+1} = r_{k+1} + \beta_kp_k \end{aligned} \]为了简化计算,注意到
\[r_k^Tr_{k+1} = r_k^Tp_{k-1} = r_{k+1}^Tp_k = 0,\quad k=1,2,\cdots \]这些关系在下面的定理中得到证明。从上式可以导出
\[\begin{aligned} r_{k+1} &= b-Ax_{k+1} = b - A(x_k+\alpha_kp_k) = r_{k-1} - \alpha_kAp_k\\ r_{k+1}^TAp_k &= \dfrac{1}{\alpha_k}r_{k+1}^T\alpha_kAp_k = \dfrac{1}{\alpha_k}r_{k+1}^T(r_k-r_{k+1}) = -\dfrac{1}{\alpha_k}r_{k+1}^Tr_{k+1}\\ p_k^TAp_k &= \dfrac{1}{\alpha_k}p_k^T\alpha_kAp_k = \dfrac{1}{\alpha_k}p_k^T(r_k-r_{k+1}) = \dfrac{1}{\alpha_k}p_k^Tr_k\\ & = \dfrac{1}{\alpha_k}(r_{k}^T+\beta_{k-1}p_{k-1}^T)r_k = \dfrac{1}{\alpha_k}r_k^Tr_k \end{aligned} \]从而我们得到
\[\alpha_k = \dfrac{r_k^Tr_k}{p_k^TAp_k},\quad \beta_k = \dfrac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k} \]这意味着每一步迭代只需使用一次 \(A\) 计算 \(p_k^TAp_k\) .
定理 5.2.1 共轭梯度法得到的向量组 \(\{r_i\},\ \{p_i\}\) 满足
- \(p_i^Tr_j = 0,\ 0\le i
- \(r_i^Tr_j = 0,\ i\neq j,\ 0\le i,j\le k\)
- \(p_i^TAp_j = 0,\ i\neq j,\ 0\le i,j\le k\)
- \(\mathrm{span}\{r_0,\cdots,r_k\} = \mathrm{span}\{p_0,\cdots,p_k\} = \mathcal{K}(A,r_0,k+1)\) ,其中
通常称之为 Krylov 子空间.
证明 应用归纳法,可以验证它们对于 \(k=1\) 都成立,其中 (1) 式前面已经证明;对于 (2) 式,由归纳假设
\[\mathrm{span}\{r_0,\cdots,r_k\} = \mathrm{span}\{p_0,\cdots,p_k\} \]而由 (1) 有 \(r_{k+1}\) 与 \(\mathrm{span}\{p_0,\cdots,p_k\}\) 正交,即证; (3) 式利用
\[p_{k+1} = r_{k+1} + \beta_kp_k,\quad Ap_i = \dfrac{1}{\alpha_i}(r_i-r_{i+1}) \]通过定义可证;由于
\[r_{k+1} = r_k - \alpha_kAp_k\in\mathcal{K}(A,r_0,k+2)\\ p_{k+1} = r_{k+1} + \beta_kp_k\in\mathcal{K}(A,r_0,k+2) \]并且由之前的结论,两个向量组是线性无关的,因此即证.
由于每一步得到的两个向量组正交,因此至多 \(n\) 步就可以得到精确的解.
定理 5.2.2 共轭梯度法得到的近似解 \(x_k\) 满足
\[\varphi(x) = \min\{\varphi(x):x\in x_0+\mathcal{K}(A,r_0,k)\} \]或者
\[\|x_k-x_*\|_A = \min\{\|x-x_*\|_A:x\in x_0+\mathcal{K}(A,r_0,k)\} \]其中 \(\|x\|_A=\sqrt{x^TAx}\) , \(x_*\) 是方程组 \(Ax=b\) 的解, \(\mathcal{K}(A,r_0,k)\) 是 Krylov 子空间.
定理 5.3.1 若 \(A = I + B\) ,且 \(\mathrm{rank}(B) = r\) ,则共轭梯度法至多迭代 \(r+1\) 步即可得到精确解.
证明 注意到 \(\mathrm{rank}(B) = r\) 蕴含着
\[\mathrm{span}\{r_0,Ar_0,\cdots,A^kr_0\} = \mathrm{span}\{r_0,Br_0,\cdots,B^kr_0\} \]事实上 \(A^k = (I+B)^k\) ,即 \(A\) 的幂总是可以表示为 \(B\) 的幂的线性组合,而右边子空间维数不超过 \(r+1\) ,即证.
定理 5.3.2 共轭梯度法得到的近似解 \(x_k\) 满足
\[\|x_k-x_*\|_A\le 2\left(\dfrac{\sqrt{\kappa_2}-1}{\sqrt{\kappa_2}+1}\right)^k\|x_0-x_*\|_A \]其中 \(\kappa_2=\kappa_2(A)=\|A\|_2\left\|A^{-1}\right\|_2\) .