线性方程组直接解法的不足：

  • 准确度不可调整：只有精确解法
  • 复杂度高，对稀疏矩阵没有较好的策略

迭代法基本概念

类似于不动点迭代法，对方程

Ax=b\mathbf {Ax}=\mathbf b

构造与原方程等价的方程

x=g(x)\mathbf x = \mathbf g(\mathbf x)

从而迭代计算

x(k+1)=g(x(k))\mathbf x^{(k+1)}=\mathbf g(\mathbf x^{(k)})

为了使每步计算量尽可能小，要求g\mathbf g为线性函数

x(k+1)=Bx(k+1)+f\mathbf x^{(k+1)}=\mathbf {Bx}^{(k+1)}+\mathbf f

其中B\mathbf B为常矩阵，f\mathbf f为常向量

上面的方法称为一阶定常迭代法

一阶定常迭代法的构造方法：矩阵分裂法

令

A=M−N\mathbf A=\mathbf M-\mathbf N

则

Ax=b⇔Mx−Nx=b⇔x=M−1Nx+M−1b\mathbf {Ax}=\mathbf b\Leftrightarrow \mathbf{Mx}-\mathbf{Nx}=\mathbf b\Leftrightarrow \mathbf x=\mathbf M^{-1}\mathbf{Nx}+\mathbf M^{-1}\mathbf b

对M\mathbf M的选取需要考虑以下问题：

  • M\mathbf M非奇异
  • 迭代法易收敛、收敛快
  • 求解以M\mathbf M为系数矩阵的方程计算量小

迭代法判停准则

类似非线性方程求根的迭代法：

  • 残差判据：∣∣b−Ax(k)∣∣≤ε1||\mathbf b-\mathbf{Ax}^{(k)}||\leq \varepsilon_1
  • 误差判据：∣∣x(k)−x(k−1)∣∣≤ε2||\mathbf x^{(k)}-\mathbf x^{(k-1)}||\leq\varepsilon_2

一阶定常迭代法的收敛性可以如下考察：

考虑误差

e(k)=x(k)−x∗\mathbf e^{(k)}=\mathbf x^{(k)}-\mathbf x^*

由于

x∗=Bx∗+f\mathbf x^*=\mathbf {Bx}^*+\mathbf f

x(k+1)=Bx(k)+f\mathbf x^{(k+1)}=\mathbf {Bx}^{(k)}+\mathbf f

可知

e(k+1)=Be(k)\mathbf e^{(k+1)}=\mathbf {Be}^{(k)}

全局收敛等价于

lim⁡k→+∞Bk=O\lim_{k\to+\infty}\mathbf B^k=\mathbf O

为此，我们介绍下面一些补充知识

矩阵收敛基本理论

矩阵收敛定义为每个元素收敛

可以证明，所有矩阵算子范数等价

进而可以证明，矩阵序列收敛等价于误差的算子范数收敛到0

矩阵收敛是比向量收敛更强的结论：

若

lim⁡k→+∞A(k)=A\lim\limits_{k\to+\infty}\mathbf A^{(k)}=\mathbf A

则

∀x,  lim⁡k→+∞A(k)x=Ax\forall \mathbf x,\ \ \lim_{k\to+\infty}\mathbf A^{(k)}\mathbf x=\mathbf {Ax}

这可以通过算子范数的定义证明

定义（谱半径） 矩阵的谱半径ρ(A)\rho(\mathbf A)为其特征值绝对值的最大值

定理（谱半径和算子范数）

  • 对任意矩阵，ρ(A)≤∣∣A∣∣\rho(\mathbf A)\leq ||\mathbf A||
  • 对实对称矩阵，∣∣A∣∣2=ρ(A)||\mathbf A||_2=\rho(\mathbf A)

证明是容易的

定理（矩阵收敛与谱半径）

lim⁡k→+∞Bk=O⇔ρ(B)<1\lim_{k\to+\infty}\mathbf B^k=\mathbf O\Leftrightarrow \rho(\mathbf B)<1

证明

对矩阵做Jordan分解：

B=XJX−1\mathbf B=\mathbf{XJX}^{-1}

则

Bk=XJkX−1\mathbf B^k=\mathbf{XJ}^k\mathbf X^{-1}

各个Jordan块的乘方独立，只需考虑一个Jordan块

Jλ,m=[λ1λ1⋱⋱λ1λ]m×m\mathbf J_{\lambda,m}=\begin{bmatrix} \lambda&1\\ &\lambda&1\\ &&\ddots&\ddots\\ &&&\lambda&1\\ &&&&\lambda \end{bmatrix}_{m\times m}

则

Jλ,mk=[λk(k1)λk−1(k2)λk−2⋯(km−1)λk−m+1λk(k1)λk−1⋯(km−2)λk−m+2λk⋯(km−3)λk−m+3⋱⋮λk]\mathbf J_{\lambda,m}^k=\begin{bmatrix} \lambda^k&\tbinom{k}{1}\lambda^{k-1}&\tbinom{k}{2}\lambda^{k-2}&\cdots&\tbinom{k}{m-1}\lambda^{k-m+1}\\ &\lambda^k&\tbinom{k}{1}\lambda^{k-1}&\cdots&\tbinom{k}{m-2}\lambda^{k-m+2}\\ &&\lambda^k&\cdots&\tbinom{k}{m-3}\lambda^{k-m+3}\\ &&&\ddots&\vdots\\ &&&&\lambda^k \end{bmatrix}

由于

lim⁡k→+∞(kl)λk−l=lim⁡k→+∞1l!λlk!λk(k−l)!\lim_{k\to+\infty}\dbinom{k}{l}\lambda^{k-l}=\lim_{k\to+\infty}\frac{1}{l!\lambda^l}\frac{k!\lambda^k}{(k-l)!}

∣1l!λlk!λk(k−l)!∣≤∣1l!λlklλk∣\left|\frac{1}{l!\lambda^l}\frac{k!\lambda^k}{(k-l)!}\right|\leq\left|\frac{1}{l!\lambda^l}k^l\lambda^k\right|

显然此极限为00当且仅当∣λ∣<1|\lambda|<1

迭代法的收敛性

ρ(B)<1\rho(\mathbf B)<1是全局收敛的充分条件，只需对一种算子范数成立即可，且能保证对任意初值都收敛（即全局收敛）

可以类似定义迭代法的收敛阶：pp阶收敛定义为

lim⁡k→+∞∣∣e(k+1)∣∣∣∣e(k)∣∣p=c≠0\lim_{k\to+\infty}\frac{||\mathbf e^{(k+1)}||}{||\mathbf e^{(k)}||^p}=c\neq 0

可以证明一阶定常迭代法1阶收敛，且c=ρ(B)c=\rho(\mathbf B)

常见的迭代法

Jacobi迭代法

分量描述

考虑原方程

{a11x1+a12x2+a13x3=b1a21x1+a22x2+a23x3=b2a31x1+a32x2+a33x3=b3\begin{cases} a_{11}x_1+a_{12}x_2+a_{13}x_3=b_1\\ a_{21}x_1+a_{22}x_2+a_{23}x_3=b_2\\ a_{31}x_1+a_{32}x_2+a_{33}x_3=b_3 \end{cases}

对第ii个方程，在左端只保留第ii个未知量

{x1=−1a11(a12x2+a13x3)+b1a11x2=−1a22(a21x1+a23x3)+b2a22x3=−1a33(a31x1+a32x2)+b3a33\begin{cases} x_1=-\frac 1 {a_{11}}(a_{12}x_2+a_{13}x_3)+\frac{b_1}{a_{11}}\\ x_2=-\frac 1 {a_{22}}(a_{21}x_1+a_{23}x_3)+\frac{b_2}{a_{22}}\\ x_3=-\frac 1 {a_{33}}(a_{31}x_1+a_{32}x_2)+\frac{b_3}{a_{33}} \end{cases}

以此为基础构造迭代法，即

{x1(k+1)=−1a11(a12x2(k)+a13x3(k))+b1a11x2(k+1)=−1a22(a21x1(k)+a23x3(k))+b2a22x3(k+1)=−1a33(a31x1(k)+a32x2(k))+b3a33\begin{cases} x_1^{(k+1)}=-\frac 1 {a_{11}}(a_{12}x_2^{(k)}+a_{13}x_3^{(k)})+\frac{b_1}{a_{11}}\\ x_2^{(k+1)}=-\frac 1 {a_{22}}(a_{21}x_1^{(k)}+a_{23}x_3^{(k)})+\frac{b_2}{a_{22}}\\ x_3^{(k+1)}=-\frac 1 {a_{33}}(a_{31}x_1^{(k)}+a_{32}x_2^{(k)})+\frac{b_3}{a_{33}} \end{cases}

矩阵描述

将A\mathbf A的对角元素和非对角元素拆开：

A=D+N,D=[a11a22a33]\mathbf A=\mathbf D+\mathbf N,\quad \mathbf D=\begin{bmatrix}a_{11}\\&a_{22}\\&&a_{33}\end{bmatrix}

于是

x=−D−1Nx+D−1b\mathbf x=-\mathbf D^{-1}\mathbf N\mathbf x+\mathbf D^{-1}\mathbf b

Jacobi迭代每步迭代计算量相当于一次矩阵乘向量，如果是系数矩阵，计算量为Nnz(A)N_{\rm nz}(\mathbf A)次乘法

计算过程不改变矩阵A\mathbf A，因此所有矩阵元素的计算均可只做一次

Gauss-Seidel迭代法

分量描述

回忆Jacobi迭代：

{x1(k+1)=−1a11(a12x2(k)+a13x3(k))+b1a11x2(k+1)=−1a22(a21x1(k)+a23x3(k))+b2a22x3(k+1)=−1a33(a31x1(k)+a32x2(k))+b3a33\begin{cases} x_1^{(k+1)}=-\frac 1 {a_{11}}(a_{12}x_2^{(k)}+a_{13}x_3^{(k)})+\frac{b_1}{a_{11}}\\ x_2^{(k+1)}=-\frac 1 {a_{22}}(a_{21}x_1^{(k)}+a_{23}x_3^{(k)})+\frac{b_2}{a_{22}}\\ x_3^{(k+1)}=-\frac 1 {a_{33}}(a_{31}x_1^{(k)}+a_{32}x_2^{(k)})+\frac{b_3}{a_{33}} \end{cases}

这样的迭代不方便原位计算，如果我们在计算后面的分量的(k+1)(k+1)次迭代时，使用前面的分量已计算完的(k+1)(k+1)次迭代，则得到

{x1(k+1)=−1a11(a12x2(k)+a13x3(k))+b1a11x2(k+1)=−1a22(a21x1(k+1)+a23x3(k))+b2a22x3(k+1)=−1a33(a31x1(k+1)+a32x2(k+1))+b3a33\begin{cases} x_1^{(k+1)}=-\frac 1 {a_{11}}(a_{12}x_2^{(k)}+a_{13}x_3^{(k)})+\frac{b_1}{a_{11}}\\ x_2^{(k+1)}=-\frac 1 {a_{22}}(a_{21}x_1^{(k+1)}+a_{23}x_3^{(k)})+\frac{b_2}{a_{22}}\\ x_3^{(k+1)}=-\frac 1 {a_{33}}(a_{31}x_1^{(k+1)}+a_{32}x_2^{(k+1)})+\frac{b_3}{a_{33}} \end{cases}

显然这样的迭代法仍然满足不动点性质

矩阵描述

上面的过程本质上相当于把原来的矩阵分成了对角、严格上三角、严格下三角三部分：

A=D+L+U\mathbf A=\mathbf D+{\mathbf L}+{\mathbf U}

然后对

Dx+Lx+Ux=b\mathbf {Dx}+{\mathbf L}\mathbf x+{\mathbf U}\mathbf x=\mathbf b

构造迭代

(D+L)x(k+1)=−Ux(k)+b(\mathbf {D}+{\mathbf L})\mathbf x^{(k+1)}=-{\mathbf U}\mathbf x^{(k)}+\mathbf b

即

x(k+1)=−(D+L)−1Ux(k)+(D+L)−1b\mathbf x^{(k+1)}=-(\mathbf {D}+{\mathbf L})^{-1}{\mathbf U}\mathbf x^{(k)}+(\mathbf {D}+{\mathbf L})^{-1}\mathbf b

这对应矩阵分裂法中的M\mathbf M为A\mathbf A的（不严格的）下三角部分D+L\mathbf D+\mathbf L

SOR迭代法

SOR迭代法全称为逐次超松弛（Successive Over Relaxation）

在G-S迭代法的基础上引入松弛因子ω\omega，本质上为G-S法算得的迭代结果与迭代前的结果进行加权平均

具体如下：

{x1(k+1)=(1−ω)x1(k)+ω[−1a11(a12x2(k)+a13x3(k))+b1a11]x2(k+1)=(1−ω)x1(k)+ω[−1a22(a21x1(k+1)+a23x3(k))+b2a22]x3(k+1)=(1−ω)x1(k)+ω[−1a33(a31x1(k+1)+a32x2(k+1))+b3a33]\begin{cases} x_1^{(k+1)}=(1-\omega)x_1^{(k)}+\omega\left[-\frac 1 {a_{11}}(a_{12}x_2^{(k)}+a_{13}x_3^{(k)})+\frac{b_1}{a_{11}}\right]\\ x_2^{(k+1)}=(1-\omega)x_1^{(k)}+\omega\left[-\frac 1 {a_{22}}(a_{21}x_1^{(k+1)}+a_{23}x_3^{(k)})+\frac{b_2}{a_{22}}\right]\\ x_3^{(k+1)}=(1-\omega)x_1^{(k)}+\omega\left[-\frac 1 {a_{33}}(a_{31}x_1^{(k+1)}+a_{32}x_2^{(k+1)})+\frac{b_3}{a_{33}}\right] \end{cases}

矩阵表示为

x(k+1)=(1−ω)x(k)+ωD−1(−Lx(k+1)−Ux(k)+b)\mathbf x^{(k+1)}=(1-\omega)\mathbf x^{(k)}+\omega\mathbf D^{-1}\left(-\mathbf {Lx}^{(k+1)}-\mathbf{Ux}^{(k)}+\mathbf b\right)

(D+ωL)x(k+1)=[(1−ω)D−ωU]x(k)+ωb(\mathbf D+\omega\mathbf L)\mathbf x^{(k+1)}=\left[(1-\omega)\mathbf D-\omega\mathbf U\right]\mathbf x^{(k)}+\omega\mathbf b

x(k+1)=(D+ωL)−1[(1−ω)D−ωU]x(k)+ω(D+ωL)(−1)b\mathbf x^{(k+1)}=(\mathbf D+\omega\mathbf L)^{-1}\left[(1-\omega)\mathbf D-\omega\mathbf U\right]\mathbf x^{(k)}+\omega(\mathbf D+\omega\mathbf L)^{(-1)}\mathbf b

三种迭代法的收敛性

使用迭代矩阵判断

考虑三种迭代法的迭代矩阵：

  • Jacobi：BJ=−D−1(A−D)\mathbf B_J=-\mathbf D^{-1}(\mathbf A-\mathbf D)
  • G-S：BG=−(D+L)−1U\mathbf B_G=-(\mathbf {D}+{\mathbf L})^{-1}{\mathbf U}
  • SOR：BS=(D+ωL)−1[(1−ω)D−ωU]\mathbf B_S=(\mathbf D+\omega\mathbf L)^{-1}\left[(1-\omega)\mathbf D-\omega\mathbf U\right]

注 这里的定义均为A=L+D+U\mathbf A=\mathbf L+\mathbf D+\mathbf U

只有Jacobi迭代的迭代矩阵是可以简单求出的

定理（Jacobi充分条件）

若∣∣BJ∣∣1<1||\mathbf B_J||_1<1或∣∣BJ∣∣∞<1||\mathbf B_J||_\infty<1，则Jacobi迭代法收敛

定理（G-S充分条件）

若∣∣BJ∣∣1<1||\mathbf B_J||_1<1或∣∣BJ∣∣∞<1||\mathbf B_J||_\infty<1，则G-S迭代法收敛

注 这里的条件与Jacobi迭代相同，不需要使用G-S迭代的迭代矩阵

使用系数矩阵判断

定理（Jacobi对称正矩阵的充分必要条件）

若A\mathbf A为对角元均正的对称矩阵，则Jacobi迭代法全局收敛等价于A\mathbf A与2D−A2\mathbf D-\mathbf A都正定

定义（可约矩阵）

设A∈Rn×n, n≥2\mathbf A\in\mathbb R^{n\times n},\ n\geq2，若存在排列阵P\mathbf P，使

PTAP=[A11A12OA22]\mathbf P^T\mathbf A\mathbf P=\begin{bmatrix}\mathbf A_{11}&\mathbf A_{12}\\\mathbf O&\mathbf A_{22}\end{bmatrix}

其中A11,A22\mathbf A_{11},\mathbf A_{22}均为方阵，则称A\mathbf A为可约矩阵，反之为不可约矩阵

定理（对角占优定理）

严格对角占优矩阵或不可约的弱对角占优矩阵非奇异

定理（对角占优与Jacobi/G-S迭代法收敛性）

对满足对角占优定理的A\mathbf A，求解Ax=b\mathbf {Ax}=\mathbf b的Jacobi和G-S迭代法都收敛

定理（对角占优与SOR迭代法收敛性）

对满足对角占优定理的A\mathbf A，若松弛因子0<ω≤10<\omega\leq 1，则求解Ax=b\mathbf {Ax}=\mathbf b的SOR迭代法收敛

定理（对称正定矩阵与迭代法收敛性）

系数矩阵为对称正定矩阵时：

  • G-S迭代法收敛
  • 0<ω<20<\omega<2时的SOR迭代法收敛

定理（SOR法收敛的必要条件）

SOR法收敛的必要条件是0<ω<20<\omega<2

共轭梯度法

最速下降法理论基础

最速下降法是除了一阶定常迭代法以外的一种新的迭代法构造思路

可以解决一阶定常迭代法对大规模矩阵收敛慢的问题

将解线性方程组化为搜索极值点

此方法只对对称正（负）定矩阵有效

设A\mathbf A对称，定义函数

ϕ(x)=12xTAx−bTx\phi(\mathbf x)=\frac12\mathbf x^T\mathbf{Ax}-\mathbf b^T\mathbf x

此函数为Rn\mathbb R^n上的标量场，其梯度为

∇ϕ(x)=[∂ϕ∂x1⋯∂ϕ∂xn]=(Ax−b)T\nabla\phi(\mathbf x)=\begin{bmatrix}\frac{\partial\phi}{\partial x_1}&\cdots&\frac{\partial\phi}{\partial x_n}\end{bmatrix}=(\mathbf{Ax}-\mathbf b)^T

Hassian矩阵为

Hϕ(x)=[∂2ϕ∂x1∂x1⋯∂2ϕ∂x1∂xn⋮⋱⋮∂2ϕ∂xn∂x1⋯∂2ϕ∂xn∂xn]=A\mathbf H_\phi(\mathbf x)=\begin{bmatrix}\frac{\partial^2\phi}{\partial x_1\partial x_1}&\cdots&\frac{\partial^2\phi}{\partial x_1\partial x_n}\\\vdots&\ddots&\vdots\\\frac{\partial^2\phi}{\partial x_n\partial x_1}&\cdots&\frac{\partial^2\phi}{\partial x_n\partial x_n}\end{bmatrix}=\mathbf A

此时：

  • 若A\mathbf A正定，则x∗\mathbf x^*为全局唯一最小值点（极小值点）
  • 若A\mathbf A负定，则x∗\mathbf x^*为全局唯一最大值点（极大值点）
  • 若A\mathbf A不定，则x∗\mathbf x^*为梯度为零的非极值点，即鞍点

最速下降法

此时，假设A\mathbf A对称正定，则我们需要求ϕ(x)\phi(\mathbf x)的最小值点

我们采取直线搜索的策略，即从x(k)\mathbf x^{(k)}出发，给定搜索方向p(k)\mathbf p^{(k)}，我们希望找到在直线

x(k)+αkp(k)\mathbf x^{(k)}+\alpha_k\mathbf p^{(k)}

上的ϕ(x)\phi(\mathbf x)的最小值点

考虑到标量场的梯度为

∇ϕ(x)=(Ax−b)T\nabla\phi(\mathbf x)=(\mathbf {Ax}-\mathbf b)^T

为了方便，我们忽略梯度中转置，将梯度按列向量处理。上述直线上每一点的梯度为

∇ϕ(x)T=Ax(k)+αkAp(k)−b\nabla\phi(\mathbf x)^T=\mathbf {Ax}^{(k)}+\alpha_k\mathbf{Ap}^{(k)}-\mathbf b

令残差向量

r(k)=b−Ax(k)\mathbf r^{(k)}=\mathbf b-\mathbf {Ax}^{(k)}

则

∇ϕ(x)T=αkAp(k)−r(k)\nabla\phi(\mathbf x)^T=\alpha_k\mathbf{Ap}^{(k)}-\mathbf r^{(k)}

直线上的方向导数为

∂ϕ∂p(k)=(∇ϕ(x))p(k)=αkp(k)TAp(k)−r(k)Tp(k)\frac{\partial\phi}{\partial\mathbf p^{(k)}}=(\nabla\phi(\mathbf x))\mathbf p^{(k)}=\alpha_k\mathbf p^{(k)T}\mathbf {Ap}^{(k)}-\mathbf{r}^{(k)T}\mathbf p^{(k)}

于是，直线上的最小值点满足

αk=r(k)Tp(k)p(k)TAp(k)\alpha_k=\frac{\mathbf{r}^{(k)T}\mathbf p^{(k)}}{\mathbf p^{(k)T}\mathbf {Ap}^{(k)}}

因此，每步迭代公式为

x(k+1)=x(k)+αkp(k)\mathbf x^{(k+1)}=\mathbf x^{(k)}+\alpha_k\mathbf p^{(k)}

其中

αk=r(k)Tp(k)p(k)TAp(k)\alpha_k=\frac{\mathbf{r}^{(k)T}\mathbf p^{(k)}}{\mathbf p^{(k)T}\mathbf {Ap}^{(k)}}

在上面的方法中，搜索方向p(k)\mathbf p^{(k)}的选取是任意的，实际上，为了满足最速下降，可以直接选择搜索方向为负梯度方向，即r(k)\mathbf r^{(k)}，于是迭代公式变为

x(k+1)=x(k)+αkr(k)\mathbf x^{(k+1)}=\mathbf x^{(k)}+\alpha_k\mathbf r^{(k)}

其中

r(k)=b−Ax(k)\mathbf r^{(k)}=\mathbf b-\mathbf{Ax}^{(k)}

αk=r(k)Tr(k)r(k)TAr(k)\alpha_k=\frac{\mathbf{r}^{(k)T}\mathbf r^{(k)}}{\mathbf r^{(k)T}\mathbf {Ar}^{(k)}}

这实际上是A\mathbf A对r(k)\mathbf r^{(k)}的Rayleigh商的倒数

此方法在实际使用时，常利用残差判据，因为算法本身需要计算残差，不需要额外计算

实际判断常使用相对残差

∣∣r(k)∣∣∣∣b∣∣\frac{||\mathbf r^{(k)}||}{||\mathbf b||}

根据ϕ(x)\phi(\mathbf x)的凸性，最速下降法一定收敛

与一阶定常迭代法相同，每步核心计算量为一步矩阵与向量乘法，这里需要计算的是Ar\mathbf {Ar}

共轭梯度法

事实上，沿最速下降方向搜索并不一定是最优的，实际最速下降法经常出现 迂回 的情况。为了减少这种迂回，使用共轭梯度法：

我们将实际搜索方向在 最速下降方向 的基础上添加少许 偏移 ，而此偏移基于 上一步的实际搜索方向 ：

x(k+1)=x(k)+ξr(k)+ηp(k−1)\mathbf x^{(k+1)}=\mathbf x^{(k)}+\xi\mathbf r^{(k)}+\eta\mathbf p^{(k-1)}

其中p(k)\mathbf p^{(k)}表示实际搜索方向，上式的后两项代表了p(k)\mathbf p^{(k)}的迭代：

p(k)=r(k)+βk−1p(k−1)\mathbf p^{(k)}=\mathbf r^{(k)}+\beta_{k-1}\mathbf p^{(k-1)}

x(k+1)=x(k)+αkp(k)\mathbf x^{(k+1)}=\mathbf x^{(k)}+\alpha_k\mathbf p^{(k)}

为此，我们需要找βk−1\beta_{k-1}和αk\alpha_k，使ϕ(x(k+1))\phi(\mathbf x^{(k+1)})取在经过x(k)\mathbf x^{(k)}，由r(k)\mathbf r^{(k)}和p(k−1)\mathbf p^{(k-1)}张成的平面上的最小值（实际计算n>>2n>>2，因此这个平面仍然是高维空间中的一个小空间）

我们将ϕ(x)\phi(\mathbf x)看成是以下二元函数：

ϕ^(ξ,η)=ϕ(x(k)+ξr(k)+ηp(k−1))\hat\phi(\xi,\eta)=\phi(\mathbf x^{(k)}+\xi\mathbf r^{(k)}+\eta\mathbf p^{(k-1)})

其中三个向量均为确定值

求出上面函数的偏导数：

∂ϕ^∂ξ=(∇ϕ(x))r(k)=(r(k))TAx(k)+ξ(r(k))TAr(k)+η(r(k))TAp(k−1)−(r(k))Tb\frac{\partial\hat\phi}{\partial\xi}=\left(\nabla\phi(\mathbf x)\right)\mathbf r^{(k)}=\left(\mathbf r^{(k)}\right)^T\mathbf {Ax}^{(k)}+\xi\left(\mathbf r^{(k)}\right)^T\mathbf {Ar}^{(k)}+\eta\left(\mathbf r^{(k)}\right)^T\mathbf {Ap}^{(k-1)}-\left(\mathbf r^{(k)}\right)^T\mathbf b

=ξ(r(k))TAr(k)+η(r(k))TAp(k−1)−(r(k))Tr(k)=\xi\left(\mathbf r^{(k)}\right)^T\mathbf {Ar}^{(k)}+\eta\left(\mathbf r^{(k)}\right)^T\mathbf {Ap}^{(k-1)}-\left(\mathbf r^{(k)}\right)^T\mathbf r^{(k)}

∂ϕ^∂η=(∇ϕ(x))p(k−1)=(p(k−1))TAx(k)+ξ(p(k−1))TAr(k)+η(p(k−1))TAp(k−1)−(p(k−1))Tb\frac{\partial\hat\phi}{\partial\eta}=\left(\nabla\phi(\mathbf x)\right)\mathbf p^{(k-1)}=\left(\mathbf p^{(k-1)}\right)^T\mathbf {Ax}^{(k)}+\xi\left(\mathbf p^{(k-1)}\right)^T\mathbf {Ar}^{(k)}+\eta\left(\mathbf p^{(k-1)}\right)^T\mathbf {Ap}^{(k-1)}-\left(\mathbf p^{(k-1)}\right)^T\mathbf b

=ξ(p(k−1))TAr(k)+η(p(k−1))TAp(k−1)−(p(k−1))Tr(k)=\xi\left(\mathbf p^{(k-1)}\right)^T\mathbf {Ar}^{(k)}+\eta\left(\mathbf p^{(k-1)}\right)^T\mathbf {Ap}^{(k-1)}-\left(\mathbf p^{(k-1)}\right)^T\mathbf {r}^{(k)}

由于上一步搜索沿方向p(k−1)\mathbf p^{(k-1)}，得到了所在平面内的最小值，因此

∂ϕ∂p(k−1)(x(k))=0\frac{\partial\phi}{\partial\mathbf p^{(k-1)}}\left(\mathbf x^{(k)}\right)=0

而x(k)\mathbf x^{(k)}处ϕ\phi的梯度正好是−r(k)-\mathbf r^{(k)}，这说明

(p(k−1))Tr(k)=0\left(\mathbf p^{(k-1)}\right)^T\mathbf {r}^{(k)}=0

即：每一步搜索的搜索方向总是与搜索后的残差正交，这也符合我们对搜索后的结果是此直线上的极值点的印象

于是，平面上的最小值满足的方程

∂ϕ^∂ξ=∂ϕ^∂η=0\frac{\partial\hat\phi}{\partial\xi}=\frac{\partial\hat\phi}{\partial\eta}=0

可以化简为

{ξ(r(k))TAr(k)+η(r(k))TAp(k−1)−(r(k))Tr(k)=0ξ(p(k−1))TAr(k)+η(p(k−1))TAp(k−1)=0\begin{cases} \xi\left(\mathbf r^{(k)}\right)^T\mathbf {Ar}^{(k)}+\eta\left(\mathbf r^{(k)}\right)^T\mathbf {Ap}^{(k-1)}-\left(\mathbf r^{(k)}\right)^T\mathbf r^{(k)}=0\\ \xi\left(\mathbf p^{(k-1)}\right)^T\mathbf {Ar}^{(k)}+\eta\left(\mathbf p^{(k-1)}\right)^T\mathbf {Ap}^{(k-1)}=0 \end{cases}

可知最优情况下的

(ηξ)opt=−(p(k−1))TAr(k)(p(k−1))TAp(k−1)\left(\frac\eta\xi\right)_{\rm opt}=-\frac{\left(\mathbf p^{(k-1)}\right)^T\mathbf {Ar}^{(k)}}{\left(\mathbf p^{(k-1)}\right)^T\mathbf {Ap}^{(k-1)}}

注意到η/ξ\eta/\xi就是新的搜索方向中的βk−1\beta_{k-1}，因此

βk−1=−(p(k−1))TAr(k)(p(k−1))TAp(k−1)\beta_{k-1}=-\frac{\left(\mathbf p^{(k-1)}\right)^T\mathbf {Ar}^{(k)}}{\left(\mathbf p^{(k-1)}\right)^T\mathbf {Ap}^{(k-1)}}

于是，计算出下一步的搜索方向

p(k)=r(k)+βk−1p(k−1)\mathbf p^{(k)}=\mathbf r^{(k)}+\beta_{k-1}\mathbf p^{(k-1)}

然后，即可计算步长

αk=r(k)Tp(k)p(k)TAp(k)\alpha_k=\frac{\mathbf{r}^{(k)T}\mathbf p^{(k)}}{\mathbf p^{(k)T}\mathbf {Ap}^{(k)}}

进行迭代

x(k+1)=x(k)+αkp(k)\mathbf x^{(k+1)}=\mathbf x^{(k)}+\alpha_k\mathbf p^{(k)}

上述公式可以采用以下手段来记忆：α,β\alpha,\beta的公式都满足以下形式：

prpp\frac{\mathbf{pr}}{\mathbf{pp}}

但α\alpha和最速下降法的α\alpha相同，是一个Rayleigh商倒数的格式，于是

α=pTrpTAp\alpha=\frac{\mathbf p^T\mathbf{r}}{\mathbf p^T\mathbf {Ap}}

而β\beta对A\mathbf A是齐次的，带负号，于是

β=−pTArpTAp\beta=-\frac{\mathbf p^T\mathbf{Ar}}{\mathbf p^T\mathbf {Ap}}

同时也要注意两个公式里使用的p\mathbf p不是同一阶段的p\mathbf p，算法中几个量的更新顺序为x(k)→r(k)→βk−1→p(k)→αk→x(k+1)\mathbf x^{(k)}\to\mathbf r^{(k)}\to\beta_{k-1}\to\mathbf p^{(k)}\to\alpha_k\to\mathbf x^{(k+1)}

取p(0)=r(0)\mathbf p^{(0)}=\mathbf r^{(0)}

共轭梯度法的几个结论

相邻两步的搜索方向在A\mathbf A定义的内积下正交（不是常规意义下正交），即

(p(k))TAp(k−1)=0\left(\mathbf p^{(k)}\right)^T\mathbf {Ap}^{(k-1)}=0

x(k+1)\mathbf x^{(k+1)}一定存在于一个称为Krylov空间的仿射空间（即平移的线性子空间）中：

x(0)+Kk+1(A,r(0))=x(0)+span(r(0),Ar(0),⋯ ,Akr(0))\mathbf x^{(0)}+\mathcal K_{k+1}(\mathbf A,\mathbf r^{(0)})=\mathbf x^{(0)}+{\rm span}(\mathbf r^{(0)},\mathbf{Ar}^{(0)},\cdots,\mathbf A^k\mathbf r^{(0)})

同时可以确定，x(k+1)\mathbf x^{(k+1)}是上述仿射子空间内的最小值点，因此x(n)\mathbf x^{(n)}将得到精确解（不考虑数值误差）

共轭梯度法的特征

共轭梯度法只能应用于对称正定矩阵，且由于其经过nn步必然计算出精确值，实际上也是一种直接解法

然而，当计算次数较多时，舍入误差很严重；对于病态矩阵也可能因误差传递而无法收敛

为了改善系数矩阵的病态性，可以使用预条件技术：

选择矩阵M\mathbf M，将原问题等价转化为

M−1Ax=M−1b\mathbf M^{-1}\mathbf{Ax}=\mathbf M^{-1}\mathbf b

其中M−1A\mathbf M^{-1}\mathbf{A}的性质更好

然而，此时无法保持矩阵的对称正定性，要保证这一点，需要选择对称正定矩阵M=LLT\mathbf M=\mathbf {LL}^T，进行等价转化

L−1AL−Ty=L−1b\mathbf{L}^{-1}\mathbf{AL}^{-T}\mathbf y=\mathbf L^{-1}\mathbf b

x=L−Ty\mathbf x=\mathbf L^{-T}\mathbf y