基本概念

范数

其用于衡量向量和矩阵的大小。

常见的向量范数定义如下：

∥x∥1=∣x1∣+∣x2∣+⋯+∣xn∣∥x∥2=∣x1∣2+∣x2∣2+⋯+∣xn∣2∥x∥∞=max⁡1≤i≤n∣xi∣∥x∥p=(∑i=1n∣xi∣p)1p\begin{aligned} \Vert{\mathbf{x}}\Vert_{1} &= |x_{1}| + |x_{2}| + \dots + |x_{n}| \\ \Vert{\mathbf{x}}\Vert_{2} &= \sqrt{|x_{1}|^{2} + |x_{2}|^{2} + \dots + |x_{n}|^{2}} \\ \Vert{\mathbf{x}}\Vert_{\infty} &= \max_{1\leq i\leq n}|x_{i}| \\ \Vert{\mathbf{x}}\Vert_{p} &= \left(\sum_{i=1}^{n}|x_{i}|^{p} \right)^{\frac{1}{p}} \end{aligned}

矩阵范数的定义：对任意nn阶方阵A\mathbf{A}，若对应一个非负实数∥A∥\Vert\mathbf{A}\Vert，满足：

  1. ∥A∥≥0\Vert\mathbf{A}\Vert \geq 0，当且仅当A=0\mathbf{A} = 0时等号成立；
  2. 对任意实数α\alpha，∥αA∥=∣α∣⋅∥A∥\Vert\alpha\mathbf{A}\Vert = |\alpha| \cdot \Vert\mathbf{A}\Vert；
  3. 对任意两个nn阶方阵A\mathbf{A}和B\mathbf{B}，∥A+B∥≤∥A∥+∥B∥\Vert\mathbf{A}+\mathbf{B}\Vert \leq \Vert\mathbf{A}\Vert + \Vert\mathbf{B}\Vert；
  4. 对任意两个nn阶方阵A\mathbf{A}和B\mathbf{B}，∥AB∥≤∥A∥⋅∥B∥\Vert\mathbf{A}\mathbf{B}\Vert \leq \Vert\mathbf{A}\Vert \cdot \Vert\mathbf{B}\Vert。

与向量范数的定义比较，前三个性质只是向量范数的推广，第四个性质则是矩阵乘法的要求（矩阵范数的相容条件）。

为了介绍常见的矩阵范数，这里需要引入谱半径的定义：设A\mathbf{A}是nn阶矩阵，则称ρ(A)=max⁡1≤i≤n∣λi∣\rho(\mathbf{A})=\max_{1\leq i \leq n} |\lambda_{i}|为A\mathbf{A}的谱半径，这里的λi\lambda_{i}为A\mathbf{A}的特征值。

谱半径=矩阵按模的最大特征值。

根据矩阵范数的定义，可以定义出矩阵的算子范数：

∥A∥p=max⁡x∈Rn∥Ax∥p∥x∥p(1)\Vert\mathbf{A}\Vert_{p} = \max_{\mathbf{x}\in \mathbb{R}^{n}}\frac{\Vert\mathbf{Ax}\Vert_{p}}{\Vert\mathbf{x}\Vert_{p}} \tag{1}

该范数由向量范数导出，也称为导出范数或诱导范数，以下是一些常见的矩阵范数：

∥A∥1=max⁡1≤j≤n∑i=1n∣aij∣（沿着行方向累加后的最大值）∥A∥2=ρ(ATA)（谱范数）∥A∥∞=max⁡1≤i≤n∑j=1n∣aij∣（沿着列方向累加后的最大值）\begin{aligned} \Vert\mathbf{A}\Vert_{1} &= \max_{1\leq j \leq n}\sum_{i=1}^{n}|a_{ij}| \text{（沿着行方向累加后的最大值）}\\ \Vert\mathbf{A}\Vert_{2} &= \sqrt{\rho\left(\mathbf{A}^{\mathrm{T}}\mathbf{A}\right)} \text{（谱范数）}\\ \Vert\mathbf{A}\Vert_{\infty} &= \max_{1\leq i \leq n}\sum_{j=1}^{n}|a_{ij}| \text{（沿着列方向累加后的最大值）} \end{aligned}

单位阵的任何一种算子范数均为1。

矩阵的F\mathrm{F}范数是向量2范数的直接推广（不属于算子范数）：

∥A∥F=∑i=1n∑jn∣aij∣2（所有元素平方和的根）\Vert\mathbf{A}\Vert_{\mathrm{F}} = \sqrt{\sum_{i=1}^{n}\sum_{j}^{n}|a_{ij}|^{2}} \text{（所有元素平方和的根）}

以上这些范数都满足矩阵范数的定义要求。

由式(1)可引出相容的概念：若∥Ax∥α≤∥A∥β⋅∥x∥α\Vert\mathbf{Ax}\Vert_{\alpha} \leq \Vert\mathbf{A}\Vert_{\beta} \cdot \Vert\mathbf{x}\Vert_{\alpha}，则称矩阵范数和向量范数是相容的。

可以看出，矩阵的算子范数与相应的向量范数是相容的。

条件数

在实际问题中，对于线性方程组Ax=b\mathbf{Ax}=\mathbf{b}，数据A\mathbf{A}和b\mathbf{b}总会存在误差δA\delta \mathbf{A}和δb\delta \mathbf{b}，其会对所要求解的x\mathbf{x}产生一个扰动δx\delta \mathbf{x}。

一方面，若仅在方程组右端的b\mathbf{b}有扰动δb\delta \mathbf{b}，即

A(x+δx)=b+δb\mathbf{A}(\mathbf{x}+\delta\mathbf{x}) = \mathbf{b} + \delta\mathbf{b}

于是有

δx=A−1δb⟹∥δx∥=∥A−1δb∥⟹∥δx∥≤∥A−1∥⋅∥δb∥∵∥b∥≤∥A∥⋅∥x∥∴∥δx∥⋅∥b∥≤∥A−1∥⋅∥δb∥⋅∥A∥⋅∥x∥⟹∥δx∥∥x∥≤∥A∥∥A−1∥∥δb∥∥b∥(2)\begin{aligned} \delta\mathbf{x} &= \mathbf{A}^{-1}\delta\mathbf{b}\\ \Longrightarrow\Vert\delta\mathbf{x}\Vert &= \Vert\mathbf{A}^{-1} \delta\mathbf{b}\Vert \\ \Longrightarrow\Vert\delta\mathbf{x}\Vert &\leq \Vert\mathbf{A}^{-1}\Vert \cdot \Vert\delta\mathbf{b}\Vert \\ \because \Vert\mathbf{b}\Vert &\leq \Vert\mathbf{A}\Vert \cdot \Vert\mathbf{x}\Vert \\ \therefore \Vert\delta\mathbf{x}\Vert\cdot\Vert\mathbf{b}\Vert &\leq \Vert\mathbf{A}^{-1}\Vert \cdot \Vert\delta\mathbf{b}\Vert \cdot \Vert\mathbf{A}\Vert \cdot \Vert\mathbf{x}\Vert\\ \Longrightarrow \frac{\Vert\delta\mathbf{x}\Vert}{\Vert\mathbf{x}\Vert} &\leq \Vert\mathbf{A}\Vert \Vert\mathbf{A}^{-1}\Vert \frac{\Vert\delta\mathbf{b}\Vert}{\Vert\mathbf{b}\Vert} \quad (2) \end{aligned}

另一方面，若仅在方程组左端的A\mathbf{A}有扰动δA\delta \mathbf{A}，即

(A+δA)(x+δx)=b(\mathbf{A}+\delta\mathbf{A})(\mathbf{x}+\delta\mathbf{x}) = \mathbf{b}

于是有

Aδx+δAx+δAδx=0⟹Aδx=−δA(x+δx)⟹δx=−A−1δA(x+δx)⟹∥δx∥≤∥A−1∥⋅∥δA∥⋅∥x+δx∥⟹∥δx∥∥x+δx∥≤∥A∥∥A−1∥∥δA∥∥A∥(3)\begin{aligned} \mathbf{A} \delta\mathbf{x} &+ \delta\mathbf{A} \mathbf{x} + \delta\mathbf{A}\delta\mathbf{x} = 0\\ \Longrightarrow \mathbf{A} \delta\mathbf{x} &= -\delta\mathbf{A} (\mathbf{x} + \delta\mathbf{x})\\ \Longrightarrow \delta\mathbf{x} &= - \mathbf{A}^{-1} \delta\mathbf{A} (\mathbf{x} + \delta\mathbf{x})\\ \Longrightarrow \Vert\delta\mathbf{x}\Vert &\leq \Vert\mathbf{A}^{-1}\Vert \cdot \Vert\delta\mathbf{A}\Vert \cdot \Vert\mathbf{x} + \delta\mathbf{x}\Vert\\ \Longrightarrow \frac{\Vert\delta\mathbf{x}\Vert}{\Vert\mathbf{x} + \delta\mathbf{x}\Vert} &\leq \Vert\mathbf{A}\Vert \Vert\mathbf{A}^{-1}\Vert \frac{\Vert\delta\mathbf{A}\Vert}{\Vert\mathbf{A}\Vert} \quad (3) \end{aligned}

观察式(2)和(3)，可以看到，无论是左端A\mathbf{A}还是右端b\mathbf{b}有扰动，解x\mathbf{x}的相对误差除了受相应扰动的相对误差（∥δA∥∥A∥\dfrac{\Vert\delta\mathbf{A}\Vert}{\Vert\mathbf{A}\Vert}和∥δb∥∥b∥\dfrac{\Vert\delta\mathbf{b}\Vert}{\Vert\mathbf{b}\Vert}）影响外，还与∥A∥∥A−1∥\Vert\mathbf{A}\Vert \Vert\mathbf{A}^{-1}\Vert有关，其起到了放大倍数的作用。

于是，引入了条件数的概念：设A\mathbf{A}为nn阶非奇异矩阵，称数cond(A)=∥A∥∥A−1∥\text{cond}(\mathbf{A}) = \Vert\mathbf{A}\Vert \Vert\mathbf{A}^{-1}\Vert为线性方程组Ax=b\mathbf{A}\mathbf{x} = \mathbf{b}的条件数，或称为矩阵A\mathbf{A}的条件数。

条件数的大小与所选取的范数有关；

单位阵的条件数为1：cond(I)=1\text{cond}(\mathbf{I}) = 1；

若A\mathbf{A}对称正定，则cond2(A)=λ1(A)λn(A)\text{cond}_{2}(\mathbf{A}) = \dfrac{\lambda_{1}(\mathbf{A})}{\lambda_{n}(\mathbf{A})}，λ1(A)\lambda_{1}(\mathbf{A})和λn(A)\lambda_{n}(\mathbf{A})分别为最大、最小特征值。

对于一个确定的方程组，若系数矩阵的条件数相对的小，则称该方程组是良态的；若条件数相对的大，则称该方程组是病态的，其系数矩阵为病态矩阵。

基本迭代法

给定一个线性方程组

Ax=b(4)\mathbf{A} \mathbf{x} = \mathbf{b} \tag{4}

其中，A∈Rn×n\mathbf{A}\in\mathbb{R}^{n\times n}和b∈Rn\mathbf{b}\in\mathbb{R}^{n}已知。

假定A\mathbf{A}有分裂A=M−N\mathbf{A} = \mathbf{M} - \mathbf{N}，其中M\mathbf{M}是非奇异矩阵。那么上述线性方程组可以改写为：

Mx=Nx+b\mathbf{M}\mathbf{x} = \mathbf{N}\mathbf{x} + \mathbf{b}

或者

x=M−1Nx+M−1b⟹x=Bx+g\mathbf{x} = \mathbf{M}^{-1}\mathbf{N}\mathbf{x} + \mathbf{M}^{-1}\mathbf{b} \Longrightarrow\mathbf{x} = \mathbf{B}\mathbf{x} + \mathbf{g}

从而可以建立如下迭代公式：

Mx(k+1)=Nx(k)+b(5)\mathbf{M}\mathbf{x}^{(k+1)} = \mathbf{N}\mathbf{x}^{(k)} + \mathbf{b} \quad (5)

或者写成：

x(k+1)=Bx(k)+g(6)\mathbf{x}^{(k+1)} = \mathbf{B}\mathbf{x}^{(k)} + \mathbf{g} \quad (6)

其中，B=M−1N\mathbf{B} = \mathbf{M}^{-1}\mathbf{N}称为迭代矩阵。

基本迭代法：给定一个初值x(0)\mathbf{x}^{(0)}，按照式(6)进行迭代，可得到一个向量序列{x(k)}\{\mathbf{x}^{(k)}\}。若x(k)\mathbf{x}^{(k)}收敛于确定的向量x∗\mathbf{x}^{*}，则x∗=Bx∗+g⟹Ax∗=b\mathbf{x}^{*} = \mathbf{B}\mathbf{x}^{*} + \mathbf{g} \Longrightarrow \mathbf{A}\mathbf{x}^{*} = \mathbf{b}，即x∗\mathbf{x}^{*}为该线性方程组的解。

但这样每次迭代需要计算一个系数矩阵为M\mathbf{M}的线性方程组，因此，当M\mathbf{M}具有某些特殊性质时，系数矩阵为M\mathbf{M}的线性方程组将易于求解（即M\mathbf{M}要构造得好）。

例如，M\mathbf{M}是对角矩阵或上三角矩阵等矩阵。

令A=D−L−U\mathbf{A} = \mathbf{D} - \mathbf{L} - \mathbf{U}，其中D\mathbf{D}为对角矩阵，L\mathbf{L}为严格下三角矩阵，U\mathbf{U}为严格上三角矩阵：

D=diag(a11,a22,⋯ ,ann), L=−[00⋯0a210⋯0⋮⋮⋱⋮an1⋯an,n−10], U=−[0a12⋯a1n⋮⋮⋱⋮00⋯an−1,n0000]\mathbf{D} = \text{diag}(a_{11}, a_{22}, \cdots,a_{nn}) ,\ \mathbf{L} = -\begin{bmatrix} 0 & 0 & \cdots & 0 \\ a_{21} & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & \cdots & a_{n,n-1} & 0 \\ \end{bmatrix} ,\ \mathbf{U} = -\begin{bmatrix} 0 & a_{12} & \cdots & a_{1n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_{n-1,n} \\ 0 & 0 & 0 & 0 \\ \end{bmatrix}

根据不同的构造形式，对以下三种基本迭代方法（Jacobi迭代法、Gauss-Seidel迭代法、SOR迭代法）进行分析。

Jacobi迭代法

在式(5)中，取M=D\mathbf{M} = \mathbf{D}和N=L+U\mathbf{N} = \mathbf{L} + \mathbf{U}，则得到雅可比（Jacobi）迭代法，其迭代格式为：

Dx(k+1)=(L+U)x(k)+b（Jacobi迭代格式）\mathbf{D}\mathbf{x}^{(k+1)} = (\mathbf{L} + \mathbf{U})\mathbf{x}^{(k)} + \mathbf{b} \quad \text{（Jacobi迭代格式）}

此时，迭代矩阵为：

BJacobi=D−1(L+U)=D−1(D−A)=I−D−1A\mathbf{B}_{\text{Jacobi}} = \mathbf{D}^{-1}(\mathbf{L}+\mathbf{U}) = \mathbf{D}^{-1}(\mathbf{D}-\mathbf{A}) = \mathbf{I} - \mathbf{D}^{-1}\mathbf{A}

Jacobi迭代法算法：

  1. 给定初值x(0)\mathbf{x}^{(0)}；

  2. 对k=0,1,2,⋯k=0,1,2,\cdots，计算：

  3. xi(k+1)=1aii(bi−∑j=1,i≠jnaijxi(k))x_{i}^{(k+1)} = \dfrac{1}{a_{ii}}\left( b_{i} - \sum_{j=1,i\neq j}^{n}a_{ij}x_{i}^{(k)} \right)；

    这里写成了每个元素求解的形式，也可直接写成向量求解的形式x(k+1)=D−1(L+U)x(k)+D−1b\mathbf{x}^{(k+1)} = \mathbf{D}^{-1}(\mathbf{L} + \mathbf{U})\mathbf{x}^{(k)} + \mathbf{D}^{-1}\mathbf{b}。

  4. 若近似解达到收敛条件，则结束；否则，继续2~3步的计算。

Gauss-Seidel迭代法

若Jacobi迭代法是收敛的，那么在Jacobi迭代法算法的第3步中，将计算出来的x(k+1)\mathbf{x}^{(k+1)}的分量（xi(k+1)x_{i}^{(k+1)}）直接投入到下一个迭代方程中使用，可能会得到更好的效果。

即，将原来的∑j=1,i≠jnaijxi(k)\sum_{j=1,i\neq j}^{n}a_{ij}x_{i}^{(k)}分成了两部分叠加：∑j=1i−1aijxi(k+1)+∑j=i+1naijxi(k)=lr,ix(k+1)+ur,ix(k)\sum_{j=1}^{i-1}a_{ij}x_{i}^{(k+1)} + \sum_{j=i+1}^{n}a_{ij}x_{i}^{(k)} = \mathbf{l}_{\text{r},i} \mathbf{x}^{(k+1)} + \mathbf{u}_{\text{r},i} \mathbf{x}^{(k)}，其中，lr,i\mathbf{l}_{\text{r},i}和ur,i\mathbf{u}_{\text{r},i}分别表示L\mathbf{L}和U\mathbf{U}的第ii行向量。

这样建立起来的迭代格式即为高斯-赛德尔（Gauss-Seidel）迭代法（GS迭代），其迭代格式为：

Dx(k+1)=Lx(k+1)+Ux(k)+b(GS迭代格式1)\mathbf{D}\mathbf{x}^{(k+1)} = \mathbf{L}\mathbf{x}^{(k+1)} + \mathbf{U}\mathbf{x}^{(k)} + \mathbf{b} \quad \text{(GS迭代格式1)}

或者把x(k+1)\mathbf{x}^{(k+1)}全部移到左侧：

(D−L)x(k+1)=Ux(k)+b(GS迭代格式2)(\mathbf{D}-\mathbf{L})\mathbf{x}^{(k+1)} = \mathbf{U}\mathbf{x}^{(k)} + \mathbf{b} \quad \text{(GS迭代格式2)}

此时，迭代矩阵为：

BGS=(D−L)−1U=(D−L)−1(D−L−A)=I−(D−L)−1A\mathbf{B}_{\text{GS}} = (\mathbf{D}-\mathbf{L})^{-1}\mathbf{U} = (\mathbf{D}-\mathbf{L})^{-1}(\mathbf{D}-\mathbf{L}-\mathbf{A}) = \mathbf{I}-(\mathbf{D}-\mathbf{L})^{-1}\mathbf{A}

GS迭代算法：

  1. 给定初值x(0)\mathbf{x}^{(0)}；

  2. 对k=0,1,2,⋯k=0,1,2,\cdots，计算：

  3. xi(k+1)=1aii(bi−∑j=1i−1aijxi(k+1)−∑j=i+1naijxi(k))x_{i}^{(k+1)} = \dfrac{1}{a_{ii}}\left( b_{i} - \sum_{j=1}^{i-1}a_{ij}x_{i}^{(k+1)} - \sum_{j=i+1}^{n}a_{ij}x_{i}^{(k)} \right)；

    右侧括号中的第二项累加下标是j=i+1j=i+1；

    这里写成了每个元素求解的形式，也可直接写成向量求解的形式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}。

  4. 若近似解达到收敛条件，则结束；否则，继续2~3步的计算。

SOR迭代法

对GS迭代格式1进行改写得到：

x(k+1)=D−1(Lx(k+1)+Ux(k)+b)=x(k)+D−1(Lx(k+1)+Ux(k)−Dx(k)+b)\begin{aligned} \mathbf{x}^{(k+1)} &= \mathbf{D}^{-1}\left(\mathbf{L}\mathbf{x}^{(k+1)} + \mathbf{U}\mathbf{x}^{(k)} + \mathbf{b}\right) \\&= \mathbf{x}^{(k)} + \mathbf{D}^{-1}\left(\mathbf{L}\mathbf{x}^{(k+1)} + \mathbf{U}\mathbf{x}^{(k)} - \mathbf{D}\mathbf{x}^{(k)} + \mathbf{b}\right) \end{aligned}

这时，可以将等式右侧的第二项看作是修正量，为了获得更快的收敛效果，可以在该修正量前乘上一个参数ω\omega，即得到逐次超松弛（Successive over relaxation，SOR）迭代法。其中，ω\omega称为松弛因子，SOR迭代格式如下：

x(k+1)=x(k)+ωD−1(Lx(k+1)+Ux(k)−Dx(k)+b)(SOR迭代格式1)\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \omega\mathbf{D}^{-1}\left(\mathbf{L}\mathbf{x}^{(k+1)} + \mathbf{U}\mathbf{x}^{(k)} - \mathbf{D}\mathbf{x}^{(k)} + \mathbf{b}\right) \quad \text{(SOR迭代格式1)}

或者把x(k+1)\mathbf{x}^{(k+1)}全部移到左侧，并集中右侧的x(k)\mathbf{x}^{(k)}：

⟹(I−ωD−1L)x(k+1)=x(k)+ωD−1Ux(k)−ωx(k)+ωD−1b⟹(D−ωL)x(k+1)=Dx(k)+ωUx(k)−ωDx(k)+ωb⟹(D−ωL)x(k+1)=((1−ω)D+ωU)x(k)+ωb⟹x(k+1)=(D−ωL)−1((1−ω)D+ωU)x(k)+ω(D−ωL)−1b(SOR迭代格式2)\begin{aligned} &\Longrightarrow \left(\mathbf{I}-\omega\mathbf{D}^{-1}\mathbf{L}\right)\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \omega\mathbf{D}^{-1}\mathbf{U}\mathbf{x}^{(k)} - \omega\mathbf{x}^{(k)} + \omega\mathbf{D}^{-1}\mathbf{b} \\ &\Longrightarrow \left(\mathbf{D}-\omega\mathbf{L}\right)\mathbf{x}^{(k+1)} = \mathbf{D}\mathbf{x}^{(k)} + \omega\mathbf{U}\mathbf{x}^{(k)} - \omega\mathbf{D}\mathbf{x}^{(k)} + \omega\mathbf{b}\\ &\Longrightarrow \left(\mathbf{D}-\omega\mathbf{L}\right)\mathbf{x}^{(k+1)} = \Big((1-\omega)\mathbf{D} + \omega\mathbf{U}\Big)\mathbf{x}^{(k)} + \omega\mathbf{b}\\ &\Longrightarrow \mathbf{x}^{(k+1)} = \left(\mathbf{D}-\omega\mathbf{L}\right)^{-1}\Big((1-\omega)\mathbf{D} + \omega\mathbf{U}\Big)\mathbf{x}^{(k)} + \omega\left(\mathbf{D}-\omega\mathbf{L}\right)^{-1}\mathbf{b} \quad \text{(SOR迭代格式2)} \end{aligned}

此时，迭代矩阵为：

BSOR=(D−ωL)−1((1−ω)D+ωU)\mathbf{B}_{\text{SOR}} = \left(\mathbf{D}-\omega\mathbf{L}\right)^{-1}\Big((1-\omega)\mathbf{D} + \omega\mathbf{U}\Big)

ω=1\omega=1时，SOR迭代法即为GS迭代法。

SOR迭代算法：

  1. 给定初值x(0)\mathbf{x}^{(0)}；

  2. 对k=0,1,2,⋯k=0,1,2,\cdots，计算：

  3. xi(k+1)=xi(k)+ωaii(bi−∑j=1i−1aijxi(k+1)−∑j=inaijxi(k))x_{i}^{(k+1)} = x_{i}^{(k)}+\dfrac{\omega}{a_{ii}}\left( b_{i} - \sum_{j=1}^{i-1}a_{ij}x_{i}^{(k+1)} - \sum_{j=i}^{n}a_{ij}x_{i}^{(k)} \right)；

    右侧括号中的第二项累加下标是j=ij=i；

    这里写成了每个元素求解的形式，也可直接写成向量求解的形式，即(SOR迭代格式2)。

  4. 若近似解达到收敛条件，则结束；否则，继续2~3步的计算。

不定常迭代法

这里介绍两类最基本的不定常迭代方法：

  1. 求解对称正定线性方程组的最速下降法和共轭梯度法；

    这两种方法本质上是一种变分方法，对应于求一个二次函数的极值，也称之为极小化方法；

    预处理共轭梯度法已被广泛应用到各个领域，是求解大型稀疏对称正定方程组的最有效方法。

  2. 求解不对称线性方程组的广义极小残量法。

最速下降法

设A\mathbf{A}对称正定，线性方程组为

Ax=b(7)\mathbf{A} \mathbf{x} = \mathbf{b} \tag{7}

其中，A=(aij)∈Rn×n\mathbf{A}=(a_{ij})\in\mathbb{R}^{n\times n}，x=(x1,x2,⋯ ,xn)T\mathbf{x}=(x_{1}, x_{2}, \cdots, x_{n})^{\mathrm{T}}，b=(b1,b2,⋯ ,bn)T\mathbf{b}=(b_{1}, b_{2}, \cdots, b_{n})^{\mathrm{T}}。

首先，介绍上述线性方程组等价的变分问题。

任取x\mathbf{x}，对于给定的A\mathbf{A}和b\mathbf{b}，定义nn元二次函数φ:Rn→R\varphi:\mathbb{R}^{n}\rightarrow\mathbb{R}为：

φ(x)=12(x,Ax)−(x,b)=12∑i=1n∑j=1naijxixj−∑i=1nbixi(8)\begin{aligned} \varphi(\mathbf{x}) &= \frac{1}{2}(\mathbf{x}, \mathbf{Ax}) - (\mathbf{x}, \mathbf{b})\\ & = \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}x_{i}x_{j} - \sum_{i=1}^{n} b_{i}x_{i} \end{aligned}\tag{8}

其中，(x,y)=xTy=∑i=1nxiyi(\mathbf{x},\mathbf{y})=\mathbf{x}^{\mathrm{T}}\mathbf{y}=\sum_{i=1}^{n} x_{i}y_{i}（内积）。

该二次函数(8)具有如下性质：

  1. 对于x∈Rn\mathbf{x}\in\mathbb{R}^{n}，∂φ∂xi=12∑j=1n(aij+aji)xj−bi\dfrac{\partial\varphi}{\partial x_{i}} = \dfrac{1}{2} \sum_{j=1}^{n} \left(a_{ij} + a_{ji}\right)x_{j} - b_{i}，即φ(x)\varphi(\mathbf{x})的梯度为：

    ∇φ(x)=12(A+AT)x−b=Ax−b\nabla\varphi(\mathbf{x})=\frac{1}{2}(\mathbf{A} + \mathbf{A}^{\mathrm{T}})\mathbf{x} - \mathbf{b} = \mathbf{A}\mathbf{x} - \mathbf{b}

    从这个地方就可以大致看出，φ(x)\varphi(\mathbf{x})的极值点就是式(7)的解。

  2. 对于x,y∈Rn\mathbf{x},\mathbf{y}\in\mathbb{R}^{n}，a∈Ra\in\mathbb{R}，有：

    φ(x+ay)=12(x+ay,A(x+ay))−(x+ay,b)=12(x,Ax)+12(x,aAy)+12(ay,Ax)+12(ay,aAy)−(x,b)−(ay,b)=φ(x)+a(y,Ax)−a(y,b)+a22(y,Ay)=φ(x)+a(y,Ax−b)+a22(y,Ay)\begin{aligned} \varphi(\mathbf{x}+a\mathbf{y}) &= \frac{1}{2}\big(\mathbf{x}+a\mathbf{y}, \mathbf{A}(\mathbf{x}+a\mathbf{y})\big) - (\mathbf{x}+a\mathbf{y}, \mathbf{b})\\ &= \frac{1}{2}(\mathbf{x},\mathbf{A}\mathbf{x}) + \frac{1}{2}(\mathbf{x},a\mathbf{A}\mathbf{y}) + \frac{1}{2}(a\mathbf{y},\mathbf{A}\mathbf{x}) + \frac{1}{2}(a\mathbf{y},a\mathbf{A}\mathbf{y}) - (\mathbf{x}, \mathbf{b}) - (a\mathbf{y}, \mathbf{b})\\ &= \varphi(\mathbf{x}) + a(\mathbf{y}, \mathbf{A}\mathbf{x}) - a(\mathbf{y}, \mathbf{b}) + \frac{a^{2}}{2}(\mathbf{y},\mathbf{A}\mathbf{y})\\ &= \varphi(\mathbf{x}) + a(\mathbf{y},\mathbf{A}\mathbf{x}-\mathbf{b}) + \frac{a^{2}}{2}(\mathbf{y},\mathbf{A}\mathbf{y}) \end{aligned}

  3. 设x∗=A−1b\mathbf{x}^{*} = \mathbf{A}^{-1}\mathbf{b}为式(5)的解，那么φ(x∗)=12(x∗,b)−(x∗,b)=−12(x∗,b)\varphi(\mathbf{x}^{*}) = \dfrac{1}{2}(\mathbf{x}^{*}, \mathbf{b}) - (\mathbf{x}^{*}, \mathbf{b}) = -\dfrac{1}{2}(\mathbf{x}^{*}, \mathbf{b})，且对于x∈Rn\mathbf{x}\in\mathbb{R}^{n}，有：

    φ(x)−φ(x∗)=12(x,Ax)−(x,b)+12(x∗,b)=12(x,Ax)−(x,Ax∗)+12(x∗,Ax∗)=12(x,Ax)−12(x,Ax∗)−12(x,Ax∗)+12(x∗,Ax∗)=12(x,A(x−x∗))−12(x−x∗,Ax∗)=12(x,A(x−x∗))−12(x∗,A(x−x∗))=12(x−x∗,A(x−x∗))\begin{aligned} \varphi(\mathbf{x}) - \varphi(\mathbf{x}^{*}) &= \frac{1}{2}(\mathbf{x}, \mathbf{Ax}) - (\mathbf{x}, \mathbf{b}) + \frac{1}{2}(\mathbf{x}^{*}, \mathbf{b}) \\&= \frac{1}{2}(\mathbf{x}, \mathbf{Ax}) - (\mathbf{x}, \mathbf{A}\mathbf{x}^{*}) + \frac{1}{2}(\mathbf{x}^{*}, \mathbf{A}\mathbf{x}^{*}) \\&= \frac{1}{2}(\mathbf{x}, \mathbf{Ax}) - \frac{1}{2}(\mathbf{x}, \mathbf{A}\mathbf{x}^{*}) - \frac{1}{2}(\mathbf{x}, \mathbf{A}\mathbf{x}^{*}) + \frac{1}{2}(\mathbf{x}^{*}, \mathbf{A}\mathbf{x}^{*}) \\&= \frac{1}{2}\big(\mathbf{x}, \mathbf{A}(\mathbf{x}-\mathbf{x}^{*})\big) - \frac{1}{2}(\mathbf{x}-\mathbf{x}^{*}, \mathbf{A}\mathbf{x}^{*}) \\&= \frac{1}{2}\big(\mathbf{x}, \mathbf{A}(\mathbf{x}-\mathbf{x}^{*})\big) - \frac{1}{2}\big(\mathbf{x}^{*}, \mathbf{A}(\mathbf{x} - \mathbf{x}^{*})\big) \\&= \frac{1}{2}\big(\mathbf{x}-\mathbf{x}^{*}, \mathbf{A}(\mathbf{x}-\mathbf{x}^{*})\big) \end{aligned}

以上三个性质的证明都需要用到A\mathbf{A}的对称性。

这里引入如下定理：设A\mathbf{A}对称正定，x∗\mathbf{x}^{*}是方程组Ax=b\mathbf{Ax}=\mathbf{b}的解的充要条件是x∗\mathbf{x}^{*}为二次函数φ(x)\varphi(\mathbf{x})的极小值点，即φ(x∗)=min⁡x∈Rnφ(x)\varphi(\mathbf{x}^{*}) = \min_{\mathbf{x}\in\mathbb{R}^{n}}\varphi(\mathbf{x})。

那么，求解Ax=b\mathbf{Ax}=\mathbf{b}的问题可以转化为求函数φ(x)\varphi(\mathbf{x})的唯一极小值点的问题。

为了找到这个极小值点，可以从任一点x(k)\mathbf{x}^{(k)}出发，沿某一指定的方向y(k)∈Rn\mathbf{y}^{(k)}\in\mathbb{R}^{n}，搜索下一个近似点x(k+1)\mathbf{x}^{(k+1)}，使得φ(x(k+1))\varphi(\mathbf{x}^{(k+1)})在该方向上达到极小值。该搜索策略，即迭代格式如下：

x(k+1)=x(k)+αky(k)(9)\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \alpha_{k}\mathbf{y}^{(k)} \tag{9}

选择y(k)\mathbf{y}^{(k)}的方式不同，将会得到不同的算法。

令y(k)\mathbf{y}^{(k)}为某一搜索方向，r(k)=b−Ax(k)\mathbf{r}^{(k)} = \mathbf{b} - \mathbf{A}\mathbf{x}^{(k)}为x(k)\mathbf{x}^{(k)}对应的残量，有如下函数（套用式(8)的性质2）：

φ(x(k+1))=φ(x(k)+αky(k))=φ(x(k))−αk(y(k),r(k))+αk22(y(k),Ay(k))(10)\begin{aligned} \varphi\left(\mathbf{x}^{(k+1)}\right) &= \varphi\left(\mathbf{x}^{(k)} + \alpha_{k}\mathbf{y}^{(k)}\right) \\&= \varphi\left(\mathbf{x}^{(k)}\right) - \alpha_{k}\left(\mathbf{y}^{(k)},\mathbf{r}^{(k)}\right) + \frac{\alpha_{k}^{2}}{2}\left(\mathbf{y}^{(k)},\mathbf{A}\mathbf{y}^{(k)}\right) \end{aligned} \quad \tag{10}

在上式中，为了使φ(x(k+1))\varphi\left(\mathbf{x}^{(k+1)}\right)在已知x(k)\mathbf{x}^{(k)}和y(k)\mathbf{y}^{(k)}下能够达到最小，则说明步长αk\alpha_{k}需要满足如下要求（只剩下αk\alpha_{k}这个变量在影响φ(x(k+1))\varphi\left(\mathbf{x}^{(k+1)}\right)）：

ddαkφ(x(k)+αky(k))=0⟹αk(y(k),Ay(k))=(y(k),r(k))⟹αk=(y(k),r(k))(y(k),Ay(k))(11)\begin{aligned} &\frac{\mathrm{d}}{\mathrm{d}\alpha_{k}}\varphi\left(\mathbf{x}^{(k)} + \alpha_{k}\mathbf{y}^{(k)}\right) = 0 \\&\Longrightarrow \alpha_{k}\left(\mathbf{y}^{(k)},\mathbf{A}\mathbf{y}^{(k)}\right) = \left(\mathbf{y}^{(k)},\mathbf{r}^{(k)}\right) \\&\Longrightarrow \alpha_{k} = \frac{\left(\mathbf{y}^{(k)},\mathbf{r}^{(k)}\right)}{\left(\mathbf{y}^{(k)},\mathbf{A}\mathbf{y}^{(k)}\right)} \quad \text{(11)} \end{aligned}

由$$\mathbf{A}$$的正定性得到：

d2dαk2φ(x(k)+αky(k))=(y(k),Ay(k))>0\frac{\mathrm{d}^{2}}{\mathrm{d}\alpha_{k}^{2}}\varphi\left(\mathbf{x}^{(k)} + \alpha_{k}\mathbf{y}^{(k)}\right) = \left(\mathbf{y}^{(k)},\mathbf{A}\mathbf{y}^{(k)}\right) > 0

因此，可以说明αk=(y(k),r(k))(y(k),Ay(k))\alpha_{k} = \dfrac{\left(\mathbf{y}^{(k)},\mathbf{r}^{(k)}\right)}{\left(\mathbf{y}^{(k)},\mathbf{A}\mathbf{y}^{(k)}\right)}正是φ(x(k)+αky(k))\varphi\left(\mathbf{x}^{(k)} + \alpha_{k}\mathbf{y}^{(k)}\right)的极小值点。将式(11)代入式(10)，得到：

φ(x(k+1))−φ(x(k))=φ(x(k)+αky(k))−φ(x(k))=−(y(k),r(k))2(y(k),Ay(k))+12(y(k),r(k))2(y(k),Ay(k))=−12(y(k),r(k))2(y(k),Ay(k))(12)\begin{aligned} \varphi\left(\mathbf{x}^{(k+1)}\right) - \varphi\left(\mathbf{x}^{(k)}\right) &= \varphi\left(\mathbf{x}^{(k)}+ \alpha_{k}\mathbf{y}^{(k)}\right) - \varphi\left(\mathbf{x}^{(k)}\right) \\&= - \frac{\left(\mathbf{y}^{(k)},\mathbf{r}^{(k)}\right)^{2}}{\left(\mathbf{y}^{(k)},\mathbf{A}\mathbf{y}^{(k)}\right)} + \frac{1}{2}\frac{\left(\mathbf{y}^{(k)},\mathbf{r}^{(k)}\right)^{2}}{\left(\mathbf{y}^{(k)},\mathbf{A}\mathbf{y}^{(k)}\right)} \\&= - \frac{1}{2}\frac{\left(\mathbf{y}^{(k)},\mathbf{r}^{(k)}\right)^{2}}{\left(\mathbf{y}^{(k)},\mathbf{A}\mathbf{y}^{(k)}\right)} \quad\text{(12)} \end{aligned}

可以看到，当(y(k),r(k))≠0\left(\mathbf{y}^{(k)},\mathbf{r}^{(k)}\right)\neq 0，即y(k)\mathbf{y}^{(k)}和r(k)\mathbf{r}^{(k)}不正交时，φ(x(k+1))<φ(x(k))\varphi\left(\mathbf{x}^{(k+1)}\right)<\varphi\left(\mathbf{x}^{(k)}\right)成立。

即随着迭代进行，φ(x(k+1))\varphi\left(\mathbf{x}^{(k+1)}\right)是不断在减小的。

由式(12)可知， φ(x(k+1))\varphi\left(\mathbf{x}^{(k+1)}\right)相对于φ(x(k))\varphi\left(\mathbf{x}^{(k)}\right)的下降量是取决于y(k)\mathbf{y}^{(k)}的方向而不是y(k)\mathbf{y}^{(k)}的大小。

原因如下：

将式(12)中的y(k)\mathbf{y}^{(k)}进行单位化，可以得到：

φ(x(k+1))−φ(x(k))=−12∥y(k)∥2(ynorm(k),r(k))2∥y(k)∥2(ynorm(k),Aynorm(k))=−12(ynorm(k),r(k))2(ynorm(k),Aynorm(k))\varphi\left(\mathbf{x}^{(k+1)}\right) - \varphi\left(\mathbf{x}^{(k)}\right) = - \frac{1}{2} \frac{\Vert\mathbf{y}^{(k)}\Vert^{2}\left(\mathbf{y}_{\text{norm}}^{(k)},\mathbf{r}^{(k)}\right)^{2}}{\Vert\mathbf{y}^{(k)}\Vert^{2} \left(\mathbf{y}_{\text{norm}}^{(k)},\mathbf{A}\mathbf{y}_{\text{norm}}^{(k)}\right)} = - \frac{1}{2} \frac{\left(\mathbf{y}_{\text{norm}}^{(k)},\mathbf{r}^{(k)}\right)^{2}}{\left(\mathbf{y}_{\text{norm}}^{(k)},\mathbf{A}\mathbf{y}_{\text{norm}}^{(k)}\right)}

其中，ynorm(k)\mathbf{y}_{\text{norm}}^{(k)}表示单位化后的向量，∥⋅∥\Vert\cdot\Vert表示向量模长。容易看出，下降量只与y(k)\mathbf{y}^{(k)}的方向，即ynorm(k)\mathbf{y}_{\text{norm}}^{(k)}有关。

那么，当y(k)\mathbf{y}^{(k)}与r(k)\mathbf{r}^{(k)}同向时，下降量最多，也可称为，下降最快。

而下降最快的方向，正好是在该点x(k)\mathbf{x}^{(k)}的负梯度方向，换言之，在点x(k)\mathbf{x}^{(k)}处下降最快的方向为在该点的负梯度方向。

结合式(8)的性质1和残量r(k)\mathbf{r}^{(k)}的定义即可看出。

从而，取y(k)=r(k)\mathbf{y}^{(k)} = \mathbf{r}^{(k)}，并代入迭代格式(9)，可以得到求的极小值点的最速下降法。其中，αk=(r(k),r(k))(r(k),Ar(k))\alpha_{k} = \dfrac{\left(\mathbf{r}^{(k)},\mathbf{r}^{(k)}\right)}{\left(\mathbf{r}^{(k)},\mathbf{A}\mathbf{r}^{(k)}\right)}。

最速下降法算法：

  1. 选定初值x(0)\mathbf{x}^{(0)}；
  2. 对k=0,1,2,⋯k=0,1,2,\cdots，直至收敛，计算：
  3. r(k)=b−Ax(k)\mathbf{r}^{(k)} = \mathbf{b} - \mathbf{A}\mathbf{x}^{(k)}；
  4. αk=(r(k),r(k))(r(k),Ar(k))\alpha_{k} = \dfrac{\left(\mathbf{r}^{(k)},\mathbf{r}^{(k)}\right)}{\left(\mathbf{r}^{(k)},\mathbf{A}\mathbf{r}^{(k)}\right)}；
  5. x(k+1)=x(k)+αkr(k)\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \alpha_{k}\mathbf{r}^{(k)}。

但一次迭代出现了两次矩阵-向量乘法运算，即Ax(k)\mathbf{A}\mathbf{x}^{(k)}和Ar(k)\mathbf{A}\mathbf{r}^{(k)}。

由于x(k+1)=x(k)+αkr(k)⟹r(k+1)=r(k)−αkAr(k)\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \alpha_{k}\mathbf{r}^{(k)}\Longrightarrow \mathbf{r}^{(k+1)} = \mathbf{r}^{(k)} - \alpha_{k}\mathbf{A}\mathbf{r}^{(k)}，那么，可以不用算Ax(k)\mathbf{A}\mathbf{x}^{(k)}就可以完成第5步。

不过，这样子计算r(k+1)\mathbf{r}^{(k+1)}会导致舍入误差逐渐累积，使得最后计算出来的r(k+1)\mathbf{r}^{(k+1)}达不到精度。因此，一般可设置每几次迭代后再计算Ax(k+1)\mathbf{A}\mathbf{x}^{(k+1)}来修正r(k+1)\mathbf{r}^{(k+1)}。

共轭梯度法

从局部角度来看，最速下降法所选择的负梯度方向是最佳的搜索方向，但从整体角度来看并非最优。

为了寻找更好的搜索方向，且每次确定搜索方向时的计算也不会太过复杂，共轭梯度法诞生了。

假设A\mathbf{A}对称正定，但此时不再沿着负梯度方向r(0),r(1),⋯r(k)\mathbf{r}^{(0)},\mathbf{r}^{(1)},\cdots\mathbf{r}^{(k)}搜索，而是寻找另外一组方向，使得经过kk次搜索后，得到近似解x(k)\mathbf{x}^{(k)}。

这里引入共轭的定义：A\mathbf{A}对称正定，若Rn\mathbb{R}^{n}中向量组{p(0),p(1),⋯ ,p(l)}\{\mathbf{p}^{(0)}, \mathbf{p}^{(1)}, \cdots, \mathbf{p}^{(l)}\}满足：

(Ap(i),p(j))=(p(i),p(j))A=0,i≠j\left(\mathbf{A}\mathbf{p}^{(i)}, \mathbf{p}^{(j)}\right)=\left(\mathbf{p}^{(i)}, \mathbf{p}^{(j)}\right)_{\mathbf{A}}=0, \quad i\neq j

则称该向量组为Rn\mathbb{R}^{n}中的一个A\mathbf{A}-共轭向量组，或称为A\mathbf{A}-正交向量组，或称这些向量是A\mathbf{A}-共轭的（两两A\mathbf{A}-共轭）。

可理解为：从另一个角度上看，这些向量是”正交“的，而该角度由A\mathbf{A}确定。

可以证明，A\mathbf{A}-共轭向量组中的向量是线性无关的。

共轭梯度法所要寻找的方向便是这样的共轭方向{p(0),p(1),⋯ ,p(l)}\{\mathbf{p}^{(0)}, \mathbf{p}^{(1)}, \cdots, \mathbf{p}^{(l)}\}，并且要求p(k)\mathbf{p}^{(k)}是由当前负梯度r(k)\mathbf{r}^{(k)}和上一帧搜索方向p(k−1)\mathbf{p}^{(k-1)}的线性组合：

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

待有时间和更加深入的理解后，将进一步说明是如何引出该构造的。

Hestenes-Stiefel公式

现在问题则变成：如何确定出这样的βk−1\beta_{k-1}？根据{p(0),p(1),⋯ ,p(l)}\{\mathbf{p}^{(0)}, \mathbf{p}^{(1)}, \cdots, \mathbf{p}^{(l)}\}是A\mathbf{A}-共轭向量组，得到

⟹(p(k−1),p(k))A=(p(k−1),r(k))A+βk−1(p(k−1),p(k−1))A⟹0=(p(k−1),r(k))A+βk−1(p(k−1),p(k−1))A⟹βk−1=−(p(k−1),r(k))A(p(k−1),p(k−1))A(14)\begin{aligned} \Longrightarrow \left(\mathbf{p}^{(k-1)}, \mathbf{p}^{(k)}\right)_{\mathbf{A}} &= \left(\mathbf{p}^{(k-1)}, \mathbf{r}^{(k)}\right)_{\mathbf{A}} + \beta_{k-1} \left(\mathbf{p}^{(k-1)}, \mathbf{p}^{(k-1)}\right)_{\mathbf{A}} \\ \Longrightarrow 0 &= \left(\mathbf{p}^{(k-1)}, \mathbf{r}^{(k)}\right)_{\mathbf{A}} + \beta_{k-1} \left(\mathbf{p}^{(k-1)}, \mathbf{p}^{(k-1)}\right)_{\mathbf{A}} \\ \Longrightarrow \beta_{k-1} &= -\frac{\left(\mathbf{p}^{(k-1)}, \mathbf{r}^{(k)}\right)_{\mathbf{A}}}{\left(\mathbf{p}^{(k-1)}, \mathbf{p}^{(k-1)}\right)_{\mathbf{A}}} \quad\text{(14)} \end{aligned}

式(14)称为Hestenes-Stiefel公式。

虽然式(14)只用了p(k)\mathbf{p}^{(k)}和p(k−1)\mathbf{p}^{(k-1)}A\mathbf{A}-共轭的条件，（这里不加证明给出）但按照式(14)进行构造的{p(0),p(1),⋯ ,p(l)}\{\mathbf{p}^{(0)}, \mathbf{p}^{(1)}, \cdots, \mathbf{p}^{(l)}\}是能够满足A\mathbf{A}-共轭的。

Crowder-Wolfe公式

观察式(14)，里面还含有矩阵-向量运算，为了消去矩阵-向量运算，可以利用如下等式：

x(k)=x(k−1)+αkp(k−1)⟹Ax(k)=Ax(k−1)+αkAp(k−1)⟹b−Ax(k)=b−Ax(k−1)−αkAp(k−1)⟹r(k)−r(k−1)=−αkAp(k−1)(15)\begin{aligned} \mathbf{x}^{(k)} &= \mathbf{x}^{(k-1)} + \alpha_{k}\mathbf{p}^{(k-1)} \\ \Longrightarrow \mathbf{A}\mathbf{x}^{(k)} &= \mathbf{A}\mathbf{x}^{(k-1)} + \alpha_{k}\mathbf{A}\mathbf{p}^{(k-1)} \\ \Longrightarrow \mathbf{b}-\mathbf{A}\mathbf{x}^{(k)} &= \mathbf{b}-\mathbf{A}\mathbf{x}^{(k-1)} - \alpha_{k}\mathbf{A}\mathbf{p}^{(k-1)} \\ \Longrightarrow \mathbf{r}^{(k)} - \mathbf{r}^{(k-1)} &= -\alpha_{k}\mathbf{A}\mathbf{p}^{(k-1)} \quad\text{(15)} \end{aligned}

可以将Ap(k−1)\mathbf{A}\mathbf{p}^{(k-1)}运算替换为只有向量之差的计算，即Ap(k−1)=1−αk(r(k)−r(k−1))\mathbf{A}\mathbf{p}^{(k-1)}=\dfrac{1}{-\alpha_{k}}\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}\right)，那么，式(14)修正为

βk−1=−(p(k−1),r(k))A(p(k−1),p(k−1))A=−1−αk(r(k)−r(k−1),r(k))1−αk(r(k)−r(k−1),p(k−1))=−(r(k)−r(k−1),r(k))(r(k)−r(k−1),p(k−1))(16)\begin{aligned} \beta_{k-1} &= -\frac{\left(\mathbf{p}^{(k-1)}, \mathbf{r}^{(k)}\right)_{\mathbf{A}}}{\left(\mathbf{p}^{(k-1)}, \mathbf{p}^{(k-1)}\right)_{\mathbf{A}}} \\ &= -\frac{\frac{1}{-\alpha_{k}}\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}, \mathbf{r}^{(k)}\right)}{\frac{1}{-\alpha_{k}}\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}, \mathbf{p}^{(k-1)}\right)} \\ &= -\frac{\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}, \mathbf{r}^{(k)}\right)}{\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}, \mathbf{p}^{(k-1)}\right)} \quad\text{(16)} \end{aligned}

式(16)称为Crowder-Wolfe公式。

Dixon公式

这里不加证明给出两个性质：

  1. (r(i),r(j))=0, i≠j\left(\mathbf{r}^{(i)},\mathbf{r}^{(j)}\right)=0,\ i\neq j；
  2. (p(i),r(j))=0, i<j\left(\mathbf{p}^{(i)},\mathbf{r}^{(j)}\right)=0,\ i< j。

根据性质1，式(16)的分子则简化为：

(r(k)−r(k−1),r(k))=(r(k),r(k))(17)\left(\mathbf{r}^{(k)}-\mathbf{r}^{(k-1)}, \mathbf{r}^{(k)}\right) = \left(\mathbf{r}^{(k)}, \mathbf{r}^{(k)}\right) \tag{17}

根据性质2，式(16)的分母则简化为：

(r(k)−r(k−1),p(k−1))=−(r(k−1),p(k−1))(18)\begin{aligned} \left(\mathbf{r}^{(k)}-\mathbf{r}^{(k-1)}, \mathbf{p}^{(k-1)}\right) &= -\left(\mathbf{r}^{(k-1)}, \mathbf{p}^{(k-1)}\right) \end{aligned}\tag{18}

那么，根据式(17)和式(18)，式(16)可修正为：

βk−1=−(r(k)−r(k−1),r(k))(r(k)−r(k−1),p(k−1))=(r(k),r(k))(r(k−1),p(k−1))(19)\begin{aligned} \beta_{k-1} &= -\frac{\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}, \mathbf{r}^{(k)}\right)}{\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}, \mathbf{p}^{(k-1)}\right)} \\ &= \frac{\left(\mathbf{r}^{(k)}, \mathbf{r}^{(k)}\right)}{\left(\mathbf{r}^{(k-1)}, \mathbf{p}^{(k-1)}\right)} \end{aligned}\tag{19}

式(19)称为Dixon公式。

Fletcher-Reeves公式

观察式(18)，由于p(k−1)=r(k−1)+βk−2p(k−2)\mathbf{p}^{(k-1)}=\mathbf{r}^{(k-1)}+\beta_{k-2}\mathbf{p}^{(k-2)}，那么式(18)可进一步简化为：

(r(k)−r(k−1),p(k−1))=−(r(k−1),p(k−1))=−(r(k−1),r(k−1)+βk−2p(k−2))=−(r(k−1),r(k−1))(20)\begin{aligned} \left(\mathbf{r}^{(k)}-\mathbf{r}^{(k-1)}, \mathbf{p}^{(k-1)}\right) &= -\left(\mathbf{r}^{(k-1)}, \mathbf{p}^{(k-1)}\right) \\ &= -\left(\mathbf{r}^{(k-1)}, \mathbf{r}^{(k-1)}+\beta_{k-2}\mathbf{p}^{(k-2)}\right) \\ &= -\left(\mathbf{r}^{(k-1)}, \mathbf{r}^{(k-1)}\right) \end{aligned}\tag{20}

那么，根据式(17)和式(20)，式(16)可修正为：

βk−1=−(r(k)−r(k−1),r(k))(r(k)−r(k−1),p(k−1))=(r(k),r(k))(r(k−1),r(k−1))(21)\begin{aligned} \beta_{k-1} &= -\frac{\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}, \mathbf{r}^{(k)}\right)}{\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}, \mathbf{p}^{(k-1)}\right)} \\ &= \frac{\left(\mathbf{r}^{(k)}, \mathbf{r}^{(k)}\right)}{\left(\mathbf{r}^{(k-1)}, \mathbf{r}^{(k-1)}\right)} \end{aligned}\tag{21}

式(21)称为Fletcher-Reeves公式。

Polak-Ribiere-Polyak公式

若只对式(16)的分母使用式(20)进行简化，则式(16)修正为：

βk−1=−(r(k)−r(k−1),r(k))(r(k)−r(k−1),p(k−1))=(r(k)−r(k−1),r(k))(r(k−1),r(k−1))(22)\begin{aligned} \beta_{k-1} &= -\frac{\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}, \mathbf{r}^{(k)}\right)}{\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}, \mathbf{p}^{(k-1)}\right)} \\ &= \frac{\left(\mathbf{r}^{(k)} - \mathbf{r}^{(k-1)}, \mathbf{r}^{(k)}\right)}{\left(\mathbf{r}^{(k-1)}, \mathbf{r}^{(k-1)}\right)} \end{aligned}\tag{22}

式(22)称为Polak-Ribiere-Polyak公式。

在实际应用中，常使用FR和PRP公式，因为这两个公式只需要梯度（残量）信息便能直接计算βk−1\beta_{k-1}。

在第一步迭代时（k=1k=1），上述公式无法计算出p(0)\mathbf{p}^{(0)}，因此，一般就直接取为负梯度方向，即p(0)=r(0)\mathbf{p}^{(0)}=\mathbf{r}^{(0)}。

为了保证所寻找的方向能够使φ(x(k))\varphi\left(\mathbf{x}^{(k)}\right)减小，那么说明p(k)\mathbf{p}^{(k)}方向不会太过偏离梯度下降的方向，即(p(k),r(k))<0\left(\mathbf{p}^{(k)},\mathbf{r}^{(k)}\right)<0。

对于步长，同最速下降法的线性搜索一致：

ddαkφ(x(k)+αkp(k))=0⟹αk(p(k),Ap(k))=(p(k),r(k))\begin{aligned} &\frac{\mathrm{d}}{\mathrm{d}\alpha_{k}}\varphi\left(\mathbf{x}^{(k)} + \alpha_{k}\mathbf{p}^{(k)}\right) = 0 \Longrightarrow \alpha_{k}\left(\mathbf{p}^{(k)},\mathbf{A}\mathbf{p}^{(k)}\right) = \left(\mathbf{p}^{(k)},\mathbf{r}^{(k)}\right) \end{aligned}

可求得：

αk=(p(k),r(k))(p(k),Ap(k))=(r(k)+βk−1p(k−1),r(k))(p(k),p(k))A=(r(k),r(k))(p(k),p(k))A(23)\begin{aligned} \alpha_{k} &= \frac{\left(\mathbf{p}^{(k)},\mathbf{r}^{(k)}\right)}{\left(\mathbf{p}^{(k)},\mathbf{A}\mathbf{p}^{(k)}\right)} \\ &= \frac{\left(\mathbf{r}^{(k)}+\beta_{k-1}\mathbf{p}^{(k-1)},\mathbf{r}^{(k)}\right)}{\left(\mathbf{p}^{(k)},\mathbf{p}^{(k)}\right)_{\mathbf{A}}} \\ &= \frac{\left(\mathbf{r}^{(k)},\mathbf{r}^{(k)}\right)}{\left(\mathbf{p}^{(k)},\mathbf{p}^{(k)}\right)_{\mathbf{A}}} \end{aligned}\tag{23}

共轭梯度算法：

  1. 选定初值x(0)\mathbf{x}^{(0)}，设p(0)=r(0)=b−Ax(0)\mathbf{p}^{(0)}=\mathbf{r}^{(0)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(0)}；
  2. 对k=0,1,2,⋯k=0,1,2,\cdots，直至收敛，计算：
  3. αk=(r(k),r(k))(p(k),p(k))A\alpha_{k} = \dfrac{\left(\mathbf{r}^{(k)},\mathbf{r}^{(k)}\right)}{\left(\mathbf{p}^{(k)},\mathbf{p}^{(k)}\right)_{\mathbf{A}}}；
  4. x(k+1)=x(k)+αkp(k)\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \alpha_{k}\mathbf{p}^{(k)}
  5. r(k+1)=b−Ax(k+1)\mathbf{r}^{(k+1)} = \mathbf{b} - \mathbf{A}\mathbf{x}^{(k+1)}；
  6. βk=(r(k+1),r(k+1))(r(k),r(k))\beta_{k}=\dfrac{\left(\mathbf{r}^{(k+1)}, \mathbf{r}^{(k+1)}\right)}{\left(\mathbf{r}^{(k)}, \mathbf{r}^{(k)}\right)}（FR公式）或者βk=(r(k+1)−r(k),r(k+1))(r(k),r(k))\beta_{k}=\dfrac{\left(\mathbf{r}^{(k+1)} - \mathbf{r}^{(k)}, \mathbf{r}^{(k+1)}\right)}{\left(\mathbf{r}^{(k)}, \mathbf{r}^{(k)}\right)}（PRP公式）；
  7. p(k+1)=r(k+1)+βkp(k)\mathbf{p}^{(k+1)}=\mathbf{r}^{(k+1)}+\beta_{k}\mathbf{p}^{(k)}；

与最速下降法对比，共轭梯度法其实就多了一步搜索方向微调的步骤，即第7步。

预处理

为了保证迭代矩阵的条件数能够小点，一般会对原方程Ax=b\mathbf{Ax}=\mathbf{b}进行变换，即M−1Ax=M−1b\mathbf{M}^{-1}\mathbf{Ax}=\mathbf{M}^{-1}\mathbf{b}，例如可以取M=diag(A)\mathbf{M}=\mathrm{diag}(\mathbf{A})，可以使求解过程收敛更快。

参考资料

[1] 同济大学计算数学教研室. 现代数值计算(第2版)[M]// 现代数值计算（第2版）. 人民邮电出版社, 2014.

[2] 共轭梯度法通俗讲义（⭐详实的内容，强推）

[2] 最优化方法复习笔记（六）共轭梯度法

[4] 数值分析6(3共轭梯度法)