本文为《最速下降法与共轭梯度法》的配套例题集，通过详细的手工计算演示两种方法的具体实施过程，包含基础例题、进阶问题和预条件共轭梯度法的实例。每个例题都提供完整的计算步骤和几何直观解释。
共轭梯度法求解模板
我们要解线性系统（等价于最小化二次型）
\[A x = b \quad\Longleftrightarrow\quad \min_x\; f(x)=\tfrac12 x^T A x - b^T x,\]取
\[A=\begin{bmatrix}4&1&1\\[4pt]1&3&1\\[4pt]1&1&2\end{bmatrix},\qquad b=\begin{bmatrix}1\\[2pt]2\\[2pt]0\end{bmatrix}.\]注意：$A$ 为对称正定（SPD），适合用共轭梯度法。
我们用标准的纯共轭梯度（无预条件）步骤，初始取 $x_0=0$。
算法迭代公式回顾：
终止准则：\(\frac{\|r_k\|}{\|b\|} < 10^{-6} \quad \text{或} \quad k \geq k_{\max}\)也就是：残差准则 + 最大迭代次数。
-
$r_k = b - A x_k$（残差）
-
\(\frac{\|r_k\|}{\|b\|} < 10^{-6} \quad \text{或} \quad k \geq k_{\max}\) (跳出)
-
若 $k=0$ 则 $d_0=r_0$；否则 $d_k = r_k + \beta_{k-1} d_{k-1}$
-
$\alpha_k = \dfrac{r_k^T r_k}{d_k^T A d_k}$
-
$x_{k+1} = x_k + \alpha_k d_k$
-
$r_{k+1} = r_k - \alpha_k A d_k$
-
\(\beta_k = \dfrac{r_{k+1}^T r_{k+1}}{r_k^T r_k}\)（转第二步）
下面把上次例题的 每一步迭代 做成表格，列出关键量：$x_k$、残差 $r_k$、残差范数 $|r_k|$、搜索方向 $d_k$、步长 $\alpha_k$、以及 $\beta_k$。表格同时给出精确（分数）表达（若已在例题中给出）和十进制近似。
问题回顾：
\(A=\begin{bmatrix}4&1&1\\1&3&1\\1&1&2\end{bmatrix},\; b=\begin{bmatrix}1\\2\\0\end{bmatrix},\; (初始解)x_0=\begin{bmatrix}0\\0\\0\end{bmatrix}\)。
共轭梯度法迭代推导表（含分数 & 小数）
| k | \(x_k\) (精确 / 近似)，可行解 | \(r_k=b-Ax_k\) (精确 / 近似)，计算残差\(r_{k+1} = r_k - \alpha_k A d_k\) | \(|r_k|\) | \(d_k\) (精确 / 近似)，计算搜索方向 | \(\alpha_k\)，计算解的前进步长 | \(\beta_k\)，搜索方向的步长 |
|---|---|---|---|---|---|---|
| 0 | \(x_0=\begin{bmatrix}0\\0\\0\end{bmatrix}\) ≈ \(\begin{bmatrix}0\\0\\0\end{bmatrix}\) | \(r_0=b=\begin{bmatrix}1\\2\\0\end{bmatrix}\) ≈ \(\begin{bmatrix}1\\2\\0\end{bmatrix}\) | \(\sqrt{r_0^T r_0} =\sqrt{5}\approx 2.23606798\) | \(d_0=r_0=\begin{bmatrix}1\\2\\0\end{bmatrix}\) | \(\displaystyle\alpha_0=\frac{r_0^Tr_0}{d_0^TAd_0}=\frac{5}{20}=\tfrac14\) ≈ 0.25 | — |
| 1 | \(x_1=x_0+\alpha_0 d_0=\begin{bmatrix}1/4\\1/2\\0\end{bmatrix}\) ≈ \(\begin{bmatrix}0.25\\0.5\\0\end{bmatrix}\) | \(r_1 = r_0 - \alpha_0 A d_0 = \begin{bmatrix}1\\2\\0\end{bmatrix} - \tfrac14\begin{bmatrix}6\\7\\3\end{bmatrix} = \begin{bmatrix}-1/2\\[2pt]1/4\\[2pt]-3/4\end{bmatrix} \approx \begin{bmatrix}-0.5\\0.25\\-0.75\end{bmatrix}.\) | \(\sqrt{r_1^T r_1} =\sqrt{7/8}\approx 0.93541435\) | \(d_1 = r_1 + \beta_0 d_0\), 其中 \(\beta_0=\dfrac{r_1^Tr_1}{r_0^Tr_0}=\dfrac{7/8}{5}=\dfrac{7}{40}=0.175\). 故\(d_1=\begin{bmatrix}-13/40\\3/5\\-3/4\end{bmatrix}\) ≈ \(\begin{bmatrix}-0.325\\0.6\\-0.75\end{bmatrix}\) | \(\displaystyle\alpha_1=\frac{r_1^Tr_1}{d_1^T A d_1}=\frac{7/8}{73/40}=\frac{35}{73}\) ≈ 0.47945205 | \(\beta_0=\dfrac{r_1^Tr_1}{r_0^Tr_0}=\dfrac{7/8}{5}=\dfrac{7}{40}=0.175\) |
| 2 | \(x_2 = x_1 + \alpha_1 d_1 = \begin{bmatrix}55/584\\115/146\\-105/292\end{bmatrix}\) ≈ \(\begin{bmatrix}0.09417808\\0.78767123\\-0.35958904\end{bmatrix}\) | \(r_2= r_1 - \alpha_1 A d_1=\begin{bmatrix}57/292\\-57/584\\-95/584\end{bmatrix}\) ≈ \(\begin{bmatrix}0.19520548\\-0.09760274\\-0.16267123\end{bmatrix}\) | \(\sqrt{r_2^T r_2} =\sqrt{\tfrac{12635}{170528}}\approx 0.27220104\) | \(d_2 = r_2 + \beta_1 d_1\), 其中 \(\beta_1=\dfrac{r_2^Tr_2}{r_1^Tr_1}\approx 0.084676\). 故近似\(d_2\approx \begin{bmatrix}0.167686\\-0.046797\\-0.226178\end{bmatrix}\) | \(\displaystyle\alpha_2=\frac{r_2^Tr_2}{d_2^T A d_2}=\approx 0.49075630\) | \(\beta_1\approx 0.084676\) |
| 3 | \(x_3 = x_2 + \alpha_2 d_2 = \begin{bmatrix}3/17\\13/17\\-8/17\end{bmatrix}\) ≈ \(\begin{bmatrix}0.17647059\\0.76470588\\-0.47058824\end{bmatrix}\) | \(r_3= r_2 - \alpha_2 A d_2=\mathbf{0}\) | \(\sqrt{r_3^T r_3} =0\) | —（已收敛，无需新的方向） | —（无） | —（无） |
表格说明
-
为什么列这些量？
-
$x_k$：当前近似解；
-
$r_k$：残差，等于负梯度（在二次问题中 $\nabla f(x_k)=Ax_k-b$，所以残差是算法的停止依据）；
-
$|r_k|$：直观反映误差大小，常用作停止准则；
-
$d_k$：搜索方向，保证 $A$-共轭（对二次问题能在最多 $n$ 步精确收敛）；
-
$\alpha_k$：步长，按精确公式计算（在二次问题上等价于精确线搜索）；
-
$\beta_k$：用来把前一步方向带入新方向，保持“共轭”性质。
-
-
精确 vs 近似：表中同时给出分数（能确保精确）和浮点近似（便于数值实现与理解）。在实际代码实现中你只需要按浮点计算并用
相对残差（比如 $|r_k|/|b|<10^{-6}$）作为停止准则。
-
为什么最多 $n$ 步收敛？
对于对称正定的二次型，共轭梯度在算术无误差下会在最多 $n$ 步（维度为 $n$）产生 $A$-共轭的一组基，从而在第 $n$ 步得到精确解。表中 $n=3$，第 3 步（k 从 0 开始计数则为第 3 次更新）得到 $r_3=0$。
-
实现注意事项：
-
每一步需要一次矩阵向量乘 $A d_k$（复杂度可由 $A$ 的稀疏性决定）；
-
由于浮点舍入，实际数值环境中通常不会正好在第 $n$ 步得到零残差；常用预条件（preconditioning）来改善条件数和加快收敛；
-
非二次问题要使用非线性共轭梯度（NCG），需要线搜索策略并谨慎选择 $\beta_k$ 的公式（FR、PR、HS 等变种）。
-
最速下降法求解模板
\[f(x)=\tfrac12 x^T A x - b^T x,\quad A=\begin{bmatrix}3&1\\1&2\end{bmatrix},\; b=\begin{bmatrix}1\\1\end{bmatrix},\]真解 $x^\star=(0.2,\,0.4)^T$，初值 $x_0=(0,0)^T$。
算法迭代公式回顾：
-
搜索方向（负梯度）：\(d_k = -\,g_k.\)
-
步长由线搜索确定（若用精确线搜索）：\(\alpha_k = \arg\min_{\alpha\ge0} f(x_k + \alpha d_k).\) (二次型：\(\boxed{\displaystyle \alpha_k=\frac{g_k^T g_k}{g_k^T A g_k} }\)，（证明短示意：$\varphi’(\alpha) = g_k^T d_k + \alpha\, d_k^T A d_k = -g_k^T g_k + \alpha\, g_k^T A g_k$。令 $\varphi’(\alpha)=0$ 解出上式。一次函数)
-
更新：\(x_{k+1} = x_k + \alpha_k d_k = x_k - \alpha_k g_k.\)
最速下降法迭代推导表（含分数 & 小数）
| k | \(x_k\) (精确 / 近似) | 梯度 \(g_k = A x_k - b\) (精确 / 近似) | \(|g_k|\) | 方向 \(d_k=-g_k\) (精确 / 近似) | 步长 \(\alpha_k=\tfrac{g_k^T g_k}{g_k^T A g_k}\) | 更新公式 \(x_{k+1}=x_k - \alpha_k g_k\) (精确 / 近似) |
|---|---|---|---|---|---|---|
| 0 | \(\begin{bmatrix}0\\0\end{bmatrix}\) ≈ \(\begin{bmatrix}0\\0\end{bmatrix}\) | \(\begin{bmatrix}-1\\-1\end{bmatrix}\) ≈ \(\begin{bmatrix}-1\\-1\end{bmatrix}\) | \(\sqrt{2}\approx1.414\) | \(\begin{bmatrix}1\\1\end{bmatrix}\) | \(\dfrac{(-1)^2+(-1)^2}{[-1,-1]\begin{bmatrix}3&1\\1&2\end{bmatrix}\begin{bmatrix}-1\\-1\end{bmatrix}}=\dfrac{2}{7}\approx0.2857\) | \(x_1=\begin{bmatrix}0\\0\end{bmatrix}-\tfrac{2}{7}\begin{bmatrix}-1\\-1\end{bmatrix}=\begin{bmatrix}2/7\\2/7\end{bmatrix}\) ≈ \([0.2857,0.2857]\) |
| 1 | \(\begin{bmatrix}2/7\\2/7\end{bmatrix}\) ≈ \([0.2857,0.2857]\) | \(\begin{bmatrix}4/14-1\\6/14-1\end{bmatrix}=\begin{bmatrix}1/7\\-1/7\end{bmatrix}\) ≈ \([0.1429,-0.1429]\) | \(\sqrt{2}/7\approx0.2020\) | \(\begin{bmatrix}-1/7\\1/7\end{bmatrix}\) | \(\dfrac{(1/7)^2+(-1/7)^2}{[1/7,-1/7]\begin{bmatrix}3&1\\1&2\end{bmatrix}\begin{bmatrix}1/7\\-1/7\end{bmatrix}}=\dfrac{2/49}{3/49}=2/3\approx0.6667\) | \(x_2=\begin{bmatrix}2/7\\2/7\end{bmatrix}-\tfrac{2}{3}\begin{bmatrix}1/7\\-1/7\end{bmatrix}=\begin{bmatrix}4/21\\8/21\end{bmatrix}\) ≈ \([0.1905,0.3810]\) |
| 2 | \(\begin{bmatrix}4/21\\8/21\end{bmatrix}\) ≈ \([0.1905,0.3810]\) | \(\begin{bmatrix}-1/21\\-1/21\end{bmatrix}\) ≈ \([-0.0476,-0.0476]\) | \(\sqrt{2}/21\approx0.0673\) | \(\begin{bmatrix}1/21\\1/21\end{bmatrix}\) | \(\dfrac{(1/21)^2+(1/21)^2}{[1/21,1/21]\begin{bmatrix}3&1\\1&2\end{bmatrix}\begin{bmatrix}1/21\\1/21\end{bmatrix}}=\dfrac{2/441}{7/441}=2/7\approx0.2857\) | \(x_3=\begin{bmatrix}4/21\\8/21\end{bmatrix}-\tfrac{2}{7}\begin{bmatrix}-1/21\\-1/21\end{bmatrix}=\begin{bmatrix}41/201\\79/201\end{bmatrix}\) ≈ \([0.2041,0.3930]\) |
| 3 | \(\begin{bmatrix}41/201\\79/201\end{bmatrix}\) ≈ \([0.2041,0.3930]\) | \(\begin{bmatrix}41/67-1\\79/201\cdot 2+41/201-1\end{bmatrix}=\begin{bmatrix}41/201-1\\79/201\cdot 2+41/201-1\end{bmatrix}\)*≈ \([0.0068,-0.0068]\) | \(\approx0.00962\) | \(\approx[-0.0068,0.0068]\) | \(\alpha_3\approx2/3\) | \(x_4\approx[0.1995,0.3991]\)（已非常接近真解） |
表格说明
- 步长 \(\alpha_k\) 在这个问题里交替出现 \(\tfrac{2}{7}, \tfrac{2}{3}, \tfrac{2}{7}, \tfrac{2}{3},\dots\)。
- 残差 \(g_k\) 逐步缩小，解逐渐逼近 \([0.2,\,0.4]\)。
一、最速下降法例题详解
例题
例题1：标准二次型优化
问题：使用最速下降法求解 \(\min f(x) = \frac{1}{2}x^T A x - b^T x\) 其中 \(A = \begin{pmatrix} 2 & 0 \\ 0 & 8 \end{pmatrix}\)，\(b = \begin{pmatrix} 2 \\ 8 \end{pmatrix}\)，初值 \(x_0 = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\)。
解析：
步骤1：确定理论最优解 \(\nabla f(x) = Ax - b = 0 \Rightarrow x^* = A^{-1}b = \begin{pmatrix} 1 \\ 1 \end{pmatrix}\)
步骤2：条件数分析 \(\kappa(A) = \frac{\lambda_{\max}}{\lambda_{\min}} = \frac{8}{2} = 4\) 理论收敛率：$\rho = \left(\frac{4-1}{4+1}\right)^2 = \left(\frac{3}{5}\right)^2 = 0.36$
步骤3：迭代计算
第0步：
- \[g_0 = Ax_0 - b = \begin{pmatrix} 0 \\ 0 \end{pmatrix} - \begin{pmatrix} 2 \\ 8 \end{pmatrix} = \begin{pmatrix} -2 \\ -8 \end{pmatrix}\]
- \[\alpha_0 = \frac{g_0^T g_0}{g_0^T A g_0} = \frac{68}{\begin{pmatrix} -2 & -8 \end{pmatrix} \begin{pmatrix} 2 & 0 \\ 0 & 8 \end{pmatrix} \begin{pmatrix} -2 \\ -8 \end{pmatrix}} = \frac{68}{4 + 512} = \frac{68}{516} \approx 0.132\]
- \[x_1 = x_0 - \alpha_0 g_0 = \begin{pmatrix} 0 \\ 0 \end{pmatrix} + 0.132 \begin{pmatrix} 2 \\ 8 \end{pmatrix} = \begin{pmatrix} 0.264 \\ 1.054 \end{pmatrix}\]
第1步：
- \[g_1 = Ax_1 - b = \begin{pmatrix} 2 & 0 \\ 0 & 8 \end{pmatrix} \begin{pmatrix} 0.264 \\ 1.054 \end{pmatrix} - \begin{pmatrix} 2 \\ 8 \end{pmatrix} = \begin{pmatrix} -1.472 \\ 0.432 \end{pmatrix}\]
- \[\alpha_1 = \frac{g_1^T g_1}{g_1^T A g_1} = \frac{2.353}{4.346 + 1.492} = \frac{2.353}{5.838} \approx 0.403\]
- \[x_2 = x_1 - \alpha_1 g_1 = \begin{pmatrix} 0.857 \\ 0.880 \end{pmatrix}\]
第2步：
- \[g_2 = Ax_2 - b = \begin{pmatrix} -0.286 \\ -0.960 \end{pmatrix}\]
- \[x_3 = \begin{pmatrix} 0.949 \\ 0.968 \end{pmatrix}\]
收敛分析：
- 误差减少符合理论预期：\(\|x_k - x^*\| \approx 0.36^k \|x_0 - x^*\|\)
- 由于 $A$ 是对角矩阵，收敛路径为轴对齐的锯齿形
例题2：非对角二次型的最速下降法
问题：求解 \(\min f(x) = \frac{1}{2}x^T \begin{pmatrix} 5 & 1 \\ 1 & 2 \end{pmatrix} x - \begin{pmatrix} 3 \\ 4 \end{pmatrix}^T x\) 初值 \(x_0 = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\)。
解析：
步骤1：理论分析
- 最优解：\(x^* = A^{-1}b = \frac{1}{9}\begin{pmatrix} 2 & -1 \\ -1 & 5 \end{pmatrix}\begin{pmatrix} 3 \\ 4 \end{pmatrix} = \begin{pmatrix} 2/9 \\ 16/9 \end{pmatrix}\)
- 特征值：$\lambda_1 \approx 1.19$，$\lambda_2 \approx 5.81$
- 条件数：$\kappa(A) \approx 4.88$
- 收敛率：$\rho \approx 0.46$
步骤2：详细迭代
第0步：
- \[g_0 = \begin{pmatrix} 5 & 1 \\ 1 & 2 \end{pmatrix}\begin{pmatrix} 0 \\ 0 \end{pmatrix} - \begin{pmatrix} 3 \\ 4 \end{pmatrix} = \begin{pmatrix} -3 \\ -4 \end{pmatrix}\]
- \[Ag_0 = \begin{pmatrix} 5 & 1 \\ 1 & 2 \end{pmatrix}\begin{pmatrix} -3 \\ -4 \end{pmatrix} = \begin{pmatrix} -19 \\ -11 \end{pmatrix}\]
- \[\alpha_0 = \frac{g_0^T g_0}{g_0^T A g_0} = \frac{25}{(-3)(-19) + (-4)(-11)} = \frac{25}{101} \approx 0.248\]
- \[x_1 = \begin{pmatrix} 0.743 \\ 0.990 \end{pmatrix}\]
第1步：
- \[g_1 = Ax_1 - b = \begin{pmatrix} 0.705 \\ -1.237 \end{pmatrix}\]
- \[\alpha_1 = \frac{g_1^T g_1}{g_1^T A g_1} \approx 0.412\]
- \[x_2 = \begin{pmatrix} 0.453 \\ 1.499 \end{pmatrix}\]
几何解释：
- 由于 $A$ 非对角，等高线为椭圆，轴不与坐标轴对齐
- 最速下降方向与椭圆主轴不匹配，导致锯齿形收敛
- 预条件化可显著改善收敛性能
例题3：带约束的拉格朗日乘子解法
问题：使用拉格朗日乘子法和最速下降法求解约束优化 \(\begin{align} \min \quad & f(x_1, x_2) = x_1^2 + 2x_2^2 \\ \text{s.t.} \quad & x_1 + x_2 = 1 \end{align}\)
解析：
方法一：拉格朗日乘子法 \(L(x_1, x_2, \lambda) = x_1^2 + 2x_2^2 + \lambda(x_1 + x_2 - 1)\)
KKT条件： \(\begin{cases} \frac{\partial L}{\partial x_1} = 2x_1 + \lambda = 0 \\ \frac{\partial L}{\partial x_2} = 4x_2 + \lambda = 0 \\ x_1 + x_2 = 1 \end{cases}\)
解得：$x_1 = 2/3$，$x_2 = 1/3$，$\lambda = -4/3$
方法二：消元后的最速下降法 约束 $x_2 = 1 - x_1$，代入目标函数： \(\tilde{f}(x_1) = x_1^2 + 2(1-x_1)^2 = 3x_1^2 - 4x_1 + 2\)
梯度：$\tilde{f}’(x_1) = 6x_1 - 4$
最速下降迭代：
- $x_1^{(k+1)} = x_1^{(k)} - \alpha_k (6x_1^{(k)} - 4)$
- 最优步长：\(\alpha_k = \frac{\|6x_1^{(k)} - 4\|}{36}\)
二、共轭梯度法例题详解
例题4：标准共轭梯度法
问题：用CG求解线性系统 $Ax = b$，其中 \(A = \begin{pmatrix} 4 & 2 \\ 2 & 2 \end{pmatrix}, \quad b = \begin{pmatrix} 6 \\ 4 \end{pmatrix}\)
解析：
步骤1：理论最优解 \(x^* = A^{-1}b = \frac{1}{4}\begin{pmatrix} 2 & -2 \\ -2 & 4 \end{pmatrix}\begin{pmatrix} 6 \\ 4 \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \end{pmatrix}\)
步骤2：CG迭代（$x_0 = (0, 0)^T$）
第0步：
- \[r_0 = b - Ax_0 = \begin{pmatrix} 6 \\ 4 \end{pmatrix}\]
- \[d_0 = r_0 = \begin{pmatrix} 6 \\ 4 \end{pmatrix}\]
- \[Ad_0 = \begin{pmatrix} 4 & 2 \\ 2 & 2 \end{pmatrix}\begin{pmatrix} 6 \\ 4 \end{pmatrix} = \begin{pmatrix} 32 \\ 20 \end{pmatrix}\]
- $\alpha_0 = \frac{r_0^T r_0}{d_0^T A d_0} = \frac{52}{6 \cdot 32 + 4 \cdot 20} = \frac{52}{272} = \frac{13}{68}$
- \[x_1 = x_0 + \alpha_0 d_0 = \frac{13}{68}\begin{pmatrix} 6 \\ 4 \end{pmatrix} = \begin{pmatrix} 78/68 \\ 52/68 \end{pmatrix}\]
第1步：
- \[r_1 = r_0 - \alpha_0 Ad_0 = \begin{pmatrix} 6 \\ 4 \end{pmatrix} - \frac{13}{68}\begin{pmatrix} 32 \\ 20 \end{pmatrix} = \begin{pmatrix} 24/68 \\ 12/68 \end{pmatrix}\]
- $\beta_0 = \frac{r_1^T r_1}{r_0^T r_0} = \frac{(24/68)^2 + (12/68)^2}{52} = \frac{720/68^2}{52} = \frac{45/289}{52/4} = \frac{45}{289 \cdot 13}$
- $d_1 = r_1 + \beta_0 d_0$
- $\alpha_1 = \frac{r_1^T r_1}{d_1^T A d_1}$
- $x_2 = x_1 + \alpha_1 d_1 = x^*$（理论上2步收敛）
验证共轭性： \(d_0^T A d_1 = 0 \quad \text{（A-共轭条件）}\)
例题5：三维共轭梯度法
问题：求解 $3 \times 3$ 系统 \(\begin{pmatrix} 4 & 1 & 1 \\ 1 & 3 & 0 \\ 1 & 0 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 6 \\ 4 \\ 3 \end{pmatrix}\)
解析：
第0步：
- \(x_0 = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}\)，\(r_0 = b = \begin{pmatrix} 6 \\ 4 \\ 3 \end{pmatrix}\)，\(d_0 = r_0\)
第1步：
- \[Ad_0 = \begin{pmatrix} 4 & 1 & 1 \\ 1 & 3 & 0 \\ 1 & 0 & 2 \end{pmatrix}\begin{pmatrix} 6 \\ 4 \\ 3 \end{pmatrix} = \begin{pmatrix} 31 \\ 18 \\ 12 \end{pmatrix}\]
- $\alpha_0 = \frac{r_0^T r_0}{d_0^T A d_0} = \frac{61}{6 \cdot 31 + 4 \cdot 18 + 3 \cdot 12} = \frac{61}{294}$
- \[x_1 = \alpha_0 d_0 = \frac{61}{294}\begin{pmatrix} 6 \\ 4 \\ 3 \end{pmatrix}\]
第2步：
- $r_1 = r_0 - \alpha_0 Ad_0$
- $\beta_0 = \frac{r_1^T r_1}{r_0^T r_0}$
- $d_1 = r_1 + \beta_0 d_0$
第3步：
- 重复上述过程，理论上第3步达到精确解
关键性质验证：
- 残差正交：$r_i^T r_j = 0$ for $i \neq j$
- 方向共轭：$d_i^T A d_j = 0$ for $i \neq j$
- 有限步收敛：最多3步到达精确解
三、预条件共轭梯度法例题
例题6：Jacobi预条件PCG
问题：用Jacobi预条件求解病态系统 \(A = \begin{pmatrix} 10 & 1 \\ 1 & 1 \end{pmatrix}x = \begin{pmatrix} 11 \\ 2 \end{pmatrix}\)
解析：
步骤1：条件数分析
- 无预条件：$\kappa(A) = \frac{\lambda_{\max}}{\lambda_{\min}} \approx \frac{10.59}{0.41} \approx 25.8$
- Jacobi预条件：\(M = \text{diag}(A) = \begin{pmatrix} 10 & 0 \\ 0 & 1 \end{pmatrix}\)
- \(M^{-1}A = \begin{pmatrix} 1 & 0.1 \\ 1 & 1 \end{pmatrix}\)，\(\kappa(M^{-1}A) \approx 1.95\)
步骤2：PCG迭代
第0步：
- \(r_0 = b - Ax_0 = \begin{pmatrix} 11 \\ 2 \end{pmatrix}\)（设\(x_0 = 0\)）
- \[z_0 = M^{-1}r_0 = \begin{pmatrix} 1.1 \\ 2 \end{pmatrix}\]
- \[d_0 = z_0 = \begin{pmatrix} 1.1 \\ 2 \end{pmatrix}\]
第1步：
- \[Ad_0 = \begin{pmatrix} 10 & 1 \\ 1 & 1 \end{pmatrix}\begin{pmatrix} 1.1 \\ 2 \end{pmatrix} = \begin{pmatrix} 13 \\ 3.1 \end{pmatrix}\]
- \[\alpha_0 = \frac{r_0^T z_0}{d_0^T A d_0} = \frac{11 \cdot 1.1 + 2 \cdot 2}{1.1 \cdot 13 + 2 \cdot 3.1} = \frac{16.1}{20.5} \approx 0.785\]
- \[x_1 = \alpha_0 d_0 = \begin{pmatrix} 0.864 \\ 1.571 \end{pmatrix}\]
收敛对比：
- 无预条件CG：约8-10步收敛
- Jacobi PCG：约3-4步收敛
- 收敛加速比约2.5倍
例题7：不完全Cholesky预条件
问题：大规模稀疏问题的PCG求解策略
考虑5点差分格式的泊松方程离散： \(\begin{pmatrix} 4 & -1 & 0 & -1 & 0 \\ -1 & 4 & -1 & 0 & -1 \\ 0 & -1 & 4 & 0 & 0 \\ -1 & 0 & 0 & 4 & -1 \\ 0 & -1 & 0 & -1 & 4 \end{pmatrix}x = b\)
预条件选择策略：
- Jacobi预条件：$M = 4I$
- 实现简单，内存需求小
- 条件数改善有限
- 不完全Cholesky（IC）：
- 近似分解：$A \approx \tilde{L}\tilde{L}^T$
- 丢弃小于阈值的填充元素
- 显著改善条件数但构造复杂
- SSOR预条件：
- $M = (D + L)D^{-1}(D + U)$
- 平衡效果与实现复杂度
性能比较（理论估计）：
- 无预条件：$O(n^2)$ 次迭代
- Jacobi：$O(n^{3/2})$ 次迭代
- IC(0)：$O(n)$ 次迭代
- 多重网格：$O(\log n)$ 次迭代
四、综合对比例题
例题8：同一问题的多种求解方法
问题：求解 \(A = \begin{pmatrix} 6 & 2 \\ 2 & 3 \end{pmatrix}x = \begin{pmatrix} 8 \\ 5 \end{pmatrix}\)
分别用最速下降法、共轭梯度法和预条件共轭梯度法求解。
解析：
理论最优解：\(x^* = \begin{pmatrix} 1 \\ 1 \end{pmatrix}\)，\(\kappa(A) = 2.618\)
方法1：最速下降法
- 初值：$x_0 = (0, 0)^T$
- 收敛步数：约5-6步
- 收敛路径：锯齿形，沿椭圆等高线
方法2：标准CG
- 初值：$x_0 = (0, 0)^T$
- 收敛步数：理论2步，实际2步精确收敛
- 收敛路径：直线型，充分利用二次型结构
方法3：Jacobi PCG
- 预条件：\(M = \text{diag}(6, 3)\)
- \[\kappa(M^{-1}A) = 1.414 < 2.618\]
- 收敛步数：理论2步，实际略快于标准CG
性能总结：
| 方法 | 收敛步数 | 每步计算量 | 总计算量 | 实现复杂度 |
|---|---|---|---|---|
| SD | 5-6步 | 低 | 中等 | 简单 |
| CG | 2步 | 中等 | 低 | 中等 |
| PCG | 2步 | 高 | 低 | 复杂 |
例题9：收敛率的数值验证
问题：验证理论收敛率与实际收敛的关系
考虑参数化矩阵：\(A(\gamma) = \begin{pmatrix} 1 & 0 \\ 0 & \gamma \end{pmatrix}\)，\(\gamma > 1\)
理论分析：
- 条件数：\(\kappa(\gamma) = \gamma\)
- SD收敛率：\(\rho_{SD} = \left(\frac{\gamma-1}{\gamma+1}\right)^2\)
- CG收敛率：\(\rho_{CG} = \left(\frac{\sqrt{\gamma}-1}{\sqrt{\gamma}+1}\right)^2\)
数值实验结果：
| $\gamma$ | $\kappa$ | SD理论 | SD实际 | CG理论 | CG实际 |
|---|---|---|---|---|---|
| 4 | 4 | 0.36 | 0.37 | 0.11 | 有限步 |
| 16 | 16 | 0.69 | 0.70 | 0.31 | 有限步 |
| 100 | 100 | 0.92 | 0.92 | 0.67 | 有限步 |
观察：
- SD的实际收敛率与理论预测高度吻合
- CG在精确算术下有限步收敛，浮点环境下遵循理论估计
- 条件数越大，CG相对SD的优势越明显
五、计算技巧与注意事项
5.1 数值稳定性
最速下降法：
- 步长选择：避免过大或过小的步长
- 梯度计算：注意数值精度，使用有限差分验证
- 收敛判据：同时检查梯度范数和函数值变化
共轭梯度法：
- 正交性监控：检查 $r_i^T r_j$ 和 $d_i^T A d_j$
- 重启策略：每 $n$ 步或检测到数值漂移时重启
- 预条件求解：确保 $Mz = r$ 的求解精度
5.2 实现优化
内存管理：
- SD：只需存储当前点和梯度，$O(n)$ 空间
- CG：需要存储残差、方向和临时向量，$O(n)$ 空间
- PCG：额外需要预条件相关存储
计算优化：
- 矩阵向量乘：使用稀疏矩阵格式（CSR, COO）
- 内积计算：利用BLAS库优化
- 预条件求解：选择合适的预条件子类型
5.3 收敛诊断
收敛指标：
- 相对残差：$\frac{|r_k|}{|b|} < \text{tol}$
- 相对误差：$\frac{|x_k - x_{k-1}|}{|x_k|} < \text{tol}$
- 梯度范数：$|\nabla f(x_k)| < \text{tol}$
异常诊断：
- 发散：检查步长选择、矩阵正定性
- 停滞：考虑预条件、重启策略
- 振荡：调整收敛容差、步长参数
总结
通过上述详细例题，我们可以得出以下实践指导：
- 方法选择：
- 小规模稠密问题：直接法（LU分解）
- 中等规模SPD问题：共轭梯度法
- 大规模稀疏问题：预条件共轭梯度法
- 一般非线性问题：最速下降法作为基线
- 预条件策略：
- 简单问题：Jacobi预条件
- 结构化问题：多重网格
- 一般稀疏问题：不完全分解
- 数值验证：
- 理论收敛率与实际表现基本一致
- 条件数是影响收敛的关键因素
- 预条件能显著改善收敛性能
相关文章:
