当前位置：首页 > news >正文

数值线性代数: 共轭梯度法

news 2025/6/9 12:59:14

本文总结线性方程组求解的相关算法，特别是共轭梯度法的原理及流程。

零、预修

0.1 LU分解

设 $\boldsymbol{A}\in \mathbb{R}^{n\times n}$ ，若对于 $k\in \left [ 1,n-1 \right ]$ ，均有 $\left | \boldsymbol{A}\left ( 1:k,1:k \right ) \right |\neq 0$ ，则存在下三角矩阵 $\boldsymbol{L} \in\mathbb{R}^{n\times n}$ 和上三角矩阵 $\boldsymbol{U} \in\mathbb{R}^{n\times n}$ ，使得 $\boldsymbol{A}=\boldsymbol{L}\boldsymbol{U}$ 。

设 $\boldsymbol{A}\in \mathbb{R}^{n\times n}$ ，若对于 $k\in \left [ 1,n \right ]$ ，均有 $\left | \boldsymbol{A}\left ( 1:k,1:k \right ) \right |\neq 0$ ，则存在唯一的下三角矩阵 $\boldsymbol{L} \in\mathbb{R}^{n\times n}$ 和上三角矩阵 $\boldsymbol{U} \in\mathbb{R}^{n\times n}$ ，使得 $\boldsymbol{A}=\boldsymbol{L}\boldsymbol{U}$ ，并且 $\left |A \right |=U\left ( 1,1 \right )U\left ( 2,2 \right )\cdots U\left ( n,n \right )$ 。

0.2 Cholesky分解

若 $\boldsymbol{A}\in \mathbb{R}^{n\times n}$ 对称正定，则存在一个对角元均为正数的下三角矩阵 $\boldsymbol{L} \in\mathbb{R}^{n\times n}$ ，使得 $\boldsymbol{A}=\boldsymbol{L}\boldsymbol{L}^{T}$ 。

一、总论：迭代法求解线性方程组的一般思路

对于非奇异矩阵 $\boldsymbol{A}\in \mathbb{R}^{n\times n}$ ， $\boldsymbol{b}\in \mathbb{R}^{n}$ ，使用迭代法求解线性方程组 $\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$ 过程中，一般需要以下流程进行：

给定一个初始向量 $\boldsymbol{x}_{0}$ ；
构造一个递推公式 $\boldsymbol{x}_{k+1}=\boldsymbol{f}\left ( \boldsymbol{x}_{k},\boldsymbol{A},\mathbf{b} \right )$ ；
不断递推 $\boldsymbol{x}_{k+1}$ ，使其近似收敛于 $\boldsymbol{x}_{*}$ ；

下表列出了若干迭代算法的迭代公式。

方法	$\boldsymbol{A}$	迭代公式	备注
Jacobi迭代	非奇异	$\boldsymbol{A}=\boldsymbol{D}-\boldsymbol{L}-\boldsymbol{U} \\ \boldsymbol{x}_{k}=\boldsymbol{D}^{-1}\left ( \boldsymbol{L}+\boldsymbol{U} \right ) \boldsymbol{x}_{k-1}+\boldsymbol{D}^{-1}\boldsymbol{b}$
Gausss-Seidel迭代	非奇异	$\boldsymbol{A}=\boldsymbol{D}-\boldsymbol{L}-\boldsymbol{U} \\ \boldsymbol{x}_{k}=\left ( \boldsymbol{D}-\boldsymbol{L }\right )^{-1}\boldsymbol{U}\boldsymbol{x}_{k-1}+\left ( \boldsymbol{D}-\boldsymbol{L} \right )^{-1}b$
SOR迭代	非奇异	$\boldsymbol{A}=\boldsymbol{D}-\boldsymbol{L}-\boldsymbol{U} \\ \boldsymbol{L}_{\omega }=\left ( \boldsymbol{D}-\omega \boldsymbol{L}\right )^{-1} \left ( \left ( 1-\omega \right )\boldsymbol{D}+\omega \boldsymbol{U} \right )\\ \boldsymbol{x}_{k+1}= \boldsymbol{L}_{\omega }\boldsymbol{x}_{k}+\omega \left ( \boldsymbol{D}-\omega \boldsymbol{L} \right )^{-1}\boldsymbol{b}$
Steepest Descent	对称正定	$\boldsymbol{r}_{k}=\boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}\\ \boldsymbol{p}_{k}=\boldsymbol{r}_{k}\\ \alpha_{k}=\frac{\boldsymbol{r}_{k}^{T}\boldsymbol{p}_{k}}{\boldsymbol{p}_{k}^{T}\boldsymbol{A}\boldsymbol{p}_{k}}\\ \boldsymbol{x}_{k+1}=\boldsymbol{x}_{k}+\alpha _{k}\boldsymbol{p}_{k}$
Conjugate Gradient	对称正定	当 $k=1$ 时 $\boldsymbol{r}_{k}=\boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}\\ \boldsymbol{p}_{k}=\boldsymbol{r}_{k}\\ \alpha_{k}=\frac{\boldsymbol{r}_{k}^{T}\boldsymbol{p}_{k}}{\boldsymbol{p}_{k}^{T}\boldsymbol{A}\boldsymbol{p}_{k}}\\ \boldsymbol{x}_{k+1}=\boldsymbol{x}_{k}+\alpha _{k}\boldsymbol{p}_{k}$ 当 $k>1$ 时 $\alpha _{k}=\frac{\boldsymbol{r}_{k}^{T}\boldsymbol{r}_{k}}{\boldsymbol{p}_{k}^{T}\boldsymbol{A}\boldsymbol{p}_{k}}\\ \boldsymbol{x}_{k+1}=\boldsymbol{x}_{k}+\alpha \boldsymbol{p}_{k} \\ \boldsymbol{r}_{k+1}=\boldsymbol{r}_{k}-\alpha _{k}\boldsymbol{A}\boldsymbol{p}_{k} \\ \beta _{k}=\frac{\boldsymbol{r}_{k+1}^{T}\boldsymbol{r}_{k+1}}{\boldsymbol{r}_{k}^{T}\boldsymbol{r}_{k}}\\ \boldsymbol{p}_{k+1}=\boldsymbol{r}_{k+1}+\beta _{k}\boldsymbol{p}_{k}$

二、Projection Method

投影法将线性方程组求解问题转换成了最优值求解问题，是求解线性方程组的一大类方法。

在投影法中，令 $\boldsymbol{r}=\boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}$ ，构造列满秩矩阵 $\mathcal{K}\in \mathbb{R}^{n\times m}$ 与 $\mathcal{L}\in \mathbb{R}^{n\times m}$ ，寻找 $\boldsymbol{\tilde{x}}\in\mathcal{K}$ ，满足Petrov-Galerkin条件，即 $\forall \boldsymbol{y}\in \mathcal{L}$ ，均有 $\mathcal{L}^{T}\left ( \boldsymbol{b}-\boldsymbol{A}\boldsymbol{\tilde{x}} \right )=\boldsymbol{0}$ 。 $\mathcal{K}$ 称为搜索空间， $\mathcal{L}$ 称为约束空间。若 $\mathcal{L}=\mathcal{K}$ 时，称为正投影算法，否则称为斜投影算法。

三、Krylov Subspace Method

Krylov子空间法本质上也是一种投影法，其核心思想是在更小维度的Krylov子空间内寻找满足精度要求的近似解。即令 $\boldsymbol{r}_{0}=\boldsymbol{b}-\boldsymbol{A}\boldsymbol{x}_{0}$ ，构造了 $m$ 维Krylov子空间 $\mathcal{K}\left ( \boldsymbol{A},\boldsymbol{r}_{0} \right )=span\left ( \boldsymbol{r}_{0} , \boldsymbol{A}\boldsymbol{r}_{0}, \boldsymbol{A}^{2} \boldsymbol{r}_{0},\cdots ,\boldsymbol{A}^{m-1}\boldsymbol{r}_{0} \right )$ ，使得 $\mathcal{L}^{T}\left (\boldsymbol{b}-\boldsymbol{A}\boldsymbol{x} \right )=\boldsymbol{0}$ 。

选择不同的 $\mathcal{L}$ ，就对应不同的Krylov子空间法。

3.1 Steepest Descent Method

3.2 Hestenes-Stiefel Conjugate Gradient Method

3.3 Preconditioned Conjugate Gradient

参考书籍

Golub G H , Loan C F V .Matrix Computations.Johns Hopkins University Press,1996.

Ford W .Numerical Linear Algebra with Applications using MATLAB. 2014.

徐树方. 数值线性代数(第二版). 北京大学出版社, 2010.

参考文献

Hestenes M R , Stiefel E L .Methods of Conjugate Gradients for Solving Linear Systems. Journal of Research of the National Bureau of Standards (United States), 1952.