miniyuan

知识点梳理

L01 误差分析

误差

  • 绝对误差： e(x∗)=x−x∗e(x^*) = x - x^* 有时也简称为误差。
  • 绝对误差限： ∣e(x∗)∣≤ε|e(x^*)| \leq \varepsilon 其中绝对误差绝对值的上界 ε\varepsilon 即为绝对误差限。
  • 相对误差： er(x∗)=e(x∗)x≈e(x∗)x∗e_r(x^*) = \frac{e(x^*)}{x} \approx \frac{e(x^*)}{x^*}
  • 相对误差限： ∣er(x∗)∣≤εr=ε∣x∗∣|e_r(x^*)| \leq \varepsilon_r = \frac{\varepsilon}{|x^*|} 其中相对误差绝对值的上界 εr\varepsilon_r 即为相对误差限。

误差传播

  • 绝对误差传播：

    e(y∗)≈∑i=1n∂f∂xi(x∗)e(xi∗)=∇f(x∗)⊤e(x∗)\begin{aligned} e(y^*) &\approx \sum_{i=1}^n \frac{\partial f}{\partial x_i}(x^*) e(x_i^*) \\ &= \nabla f(x^*)^\top \mathbf{e}(x^*) \end{aligned}

  • 相对误差传播：

    er(y∗)≈∑i=1n∂f∂xi(x∗)xi∗y∗er(xi∗)=∇ln⁡f(x∗)⊤er(x∗)\begin{aligned} e_r(y^*) &\approx \sum_{i=1}^n \frac{\partial f}{\partial x_i}(x^*) \frac{x_i^*}{y^*} e_r(x_i^*) \\ &= \nabla_{\ln} f(x^*)^\top \mathbf{e}_r(x^*) \end{aligned}

    其中

    ∇ln⁡f(x)∈Rn,(∇ln⁡f(x))i=∂ln⁡f∂ln⁡xi(x)=xif(x)∂f∂xi(x)\nabla_{\ln} f(x) \in \mathbb{R}^n, \quad \left( \nabla_{\ln} f(x) \right)_i = \frac{\partial \ln f}{\partial \ln x_i}(x) = \frac{x_i}{f(x)} \frac{\partial f}{\partial x_i}(x)

    本质上是因为 er(y∗)=y−y∗y∗≈Δ(ln⁡y)e_r(y^*) = \frac{y - y^*}{y^*} \approx \Delta (\ln y)。

有效数字

定义： x∗x^* 若可规格化为以下浮点数

±0.α1α2⋯αn×10m,α1≠0\pm 0.\alpha_1\alpha_2\cdots\alpha_n \times 10^m, \quad \alpha_1 \neq 0

其中写出的数字都是准确的，则称 x∗x^* 有 nn 位有效数字。

等价定义： x∗x^* 的绝对误差限若满足

∣e(x∗)∣≤12×10m−n|e(x^*)| \le \frac{1}{2} \times 10^{m-n}

其中 mm 为 x∗x^* 规格化表示中的指数，则称 x∗x^* 有 nn 位有效数字。

性质：

  • 绝对误差限：

    ∣e(x∗)∣≤12×10m−n|e(x^*)| \le \frac{1}{2} \times 10^{m-n}

  • 相对误差限：

    ∣er(x∗)∣≤1/2×10m−nα1×10m−1=12α1×10−(n−1)|e_r(x^*)| \le \frac{1/2 \times 10^{m-n}}{\alpha_1 \times 10^{m-1}} = \frac{1}{2\alpha_1} \times 10^{-(n-1)}

  • 相对误差限反推有效数字：

    若 x∗≠0x^* \ne 0 的相对误差限满足：

    ∣er(x∗)∣≤12(α1+1)×10−(n−1)|e_r(x^*)| \leq \frac{1}{2(\alpha_1+1)} \times 10^{-(n-1)}

    则 x∗x^* 至少有 nn 位有效数字。

数值问题和数值方法的性质

  1. 数值问题的适定性就是数值问题的解存在唯一，且关于参数连续。
  2. 数值方法的一致性就是近似函数接近准确函数。
  3. 数值方法的稳定性就是近似解关于近似参数连续。
  4. 数值方法的收敛性就是近似解接近准确解。

定义：数值计算问题是适定的，当且仅当满足：

  1. 解存在：对给定 dd，至少存在一个解 xx
  2. 解唯一：解在解空间中唯一
  3. 解关于参数连续： 参数 dd 的微小扰动 δd\delta d 仅引起解 xx 的微小变化 δx\delta x。 具体来说，设问题扰动为 F(x+δx,d+δd)=0F(x + \delta x, d + \delta d) = 0，则对任意 dd， 存在 η0>0\eta_0 > 0 和 K0>0K_0 > 0，使得当 ∥δd∥≤η0\|\delta d\| \leq \eta_0 时： ∥δx∥≤K0∥δd∥\|\delta x\| \leq K_{0} \|\delta d\|

定义：对近似问题序列 Fn(xn,dn)=0F_n(x_n, d_n) = 0，若

Fn(x,d)−F(x,d)→0(n→∞)F_n(x, d) - F(x, d) \to 0 \quad (n \to \infty)

其中 xx 是原问题的准确解，则称该序列一致。

定义：对近似问题序列 Fn(xn,dn)=0F_n(x_n, d_n) = 0，若对任意 dnd_n， 存在 η0>0\eta_0 > 0 和 K0>0K_0 > 0，使得当 ∥δdn∥≤η0\|\delta d_n\| \leq \eta_0 时：

∥δxn∥≤K0∥δdn∥\|\delta x_n\| \leq K_{0} \|\delta d_n\|

则称该序列稳定。

定义：对近似问题序列 Fn(xn,dn)=0F_n(x_n, d_n) = 0，若对任意 ε>0\varepsilon > 0， 存在 n0∈Nn_0 \in \mathbb N 和 δ>0\delta > 0，使得当 n>n0n > n_0 且 ∥δdn∥≤δ\|\delta d_n\| \leq \delta 时：

∥x(d)−xn(d+δdn)∥≤ε\|x(d) - x_n(d + \delta d_n)\| \leq \varepsilon

则称该序列收敛。

Lax-Richtmyer 定理：

对满足 Consistency\text{Consistency} 的数值方法，有：

Stability⇔Convergence\text{Stability} \Leftrightarrow \text{Convergence}

注：由于收敛性涉及到准确解，故一般不好判断。该定理可将收敛性转化为稳定性进行判断。

L02 解线性方程组的直接法（Gauss 消去法、列主元消去法）

GAUSS 消去法的矩阵视角

高斯消元的每一步可表示为左乘一个初等下三角矩阵 Mk\mathbf{M}_k，其逆的乘积构成 LU 分解中的下三角部分 L\mathbf{L}。

第 kk 步消元矩阵定义为：

Mk=I−mkek⊤\mathbf{M}_k = \mathbf{I} - \mathbf{m}_k \mathbf{e}_k^\top

其中 mk=[0,…,0,mk+1,k,…,mn,k]⊤\mathbf{m}_k = [0,\dots,0,m_{k+1,k},\dots,m_{n,k}]^\top，mj,k=aj,k(k)/ak,k(k)m_{j,k}=a_{j,k}^{(k)}/a_{k,k}^{(k)}，ek\mathbf{e}_k 为标准基向量。

易得：

Mk−1=I+mkek⊤\mathbf{M}_k^{-1} = \mathbf{I} + \mathbf{m}_k \mathbf{e}_k^\top

经过 n−1n-1 步消元：

U=Mn−1⋯M2M1A\mathbf{U} = \mathbf{M}_{n-1} \cdots \mathbf{M}_2 \mathbf{M}_1 \mathbf{A}

令：

L=M1−1M2−1⋯Mn−1−1\mathbf{L} = \mathbf{M}_1^{-1} \mathbf{M}_2^{-1} \cdots \mathbf{M}_{n-1}^{-1}

则得到 LU 分解：

A=LU\mathbf{A} = \mathbf{L} \mathbf{U}

L\mathbf{L} 为单位下三角矩阵，其第 jj 行第 kk 列（j>kj>k）的元素恰为第 kk 步的消元乘子 mj,km_{j,k}：

L=(1m211m31m321⋮⋮⋱⋱mn1mn2⋯mn,n−11)\mathbf{L} = \begin{pmatrix} 1 & & & & \\ m_{21} & 1 & & & \\ m_{31} & m_{32} & 1 & & \\ \vdots & \vdots & \ddots & \ddots & \\ m_{n1} & m_{n2} & \cdots & m_{n,n-1} & 1 \end{pmatrix}

列主元消去法

在每一步选择该列中绝对值最大的元素作主元并交换行，称为部分选主元（partial pivoting）。 由于 A\mathbf A 的非奇异性，这是存在的，否则该列会与前面的列线性相关。

L02 解线性方程组的直接法（LU 分解、Cholesky 分解）

LU 分解

一个方阵 AA （不一定可逆）分解为一个单位下三角矩阵 LL 和一个上三角矩阵 UU 的乘积：

A=LUA = LU

计算公式：

ukj=akj−∑m=1k−1lkmumj,j=k,k+1,…,nu_{kj} = a_{kj} - \sum_{m=1}^{k-1} l_{km}u_{mj}, \quad j=k,k+1,\dots,n lik=1ukk(aik−∑m=1k−1limumk),i=k+1,k+2,…,nl_{ik} = \frac{1}{u_{kk}}\left(a_{ik} - \sum_{m=1}^{k-1} l_{im}u_{mk}\right), \quad i=k+1,k+2,\dots,n

LU 分解求解线性方程组

给定线性方程组 Ax=bAx=b，若已知 A=LUA=LU，则可转化为两步求解：

  1. 前向代入（Forward Substitution）：求解 Ly=bLy = b

    yi=bi−∑j=1i−1lijyj,i=1,2,…,ny_i = b_i - \sum_{j=1}^{i-1} l_{ij}y_j, \quad i=1,2,\dots,n

  2. 后向代入（Backward Substitution）：求解 Ux=yUx = y

    xi=1uii(yi−∑j=i+1nuijxj),i=n,n−1,…,1x_i = \frac{1}{u_{ii}}\left(y_i - \sum_{j=i+1}^{n} u_{ij}x_j\right), \quad i=n,n-1,\dots,1

Cholesky 分解

若 A∈Rn×nA \in \mathbb{R}^{n \times n} 是对称正定矩阵，则存在唯一的下三角矩阵 LL（对角元为正），使得：

A=LL⊤A = LL^\top

计算公式：

lkk=akk−∑j=1k−1lkj2,k=1,2,3,…,nl_{kk} = \sqrt{a_{kk} - \sum_{j=1}^{k-1} l_{kj}^2}, \quad k=1,2,3,\dots,n lik=1lkk(aik−∑j=1k−1lijlkj),i=k+1,k+2,…,nl_{ik} = \frac{1}{l_{kk}}\left(a_{ik} - \sum_{j=1}^{k-1} l_{ij}l_{kj}\right), \quad i=k+1,k+2,\dots,n

LDL^T 分解

设 A∈Rn×nA \in \mathbb{R}^{n \times n} 是对称矩阵，且其所有顺序主子式 det⁡(Ak)≠0\det(A_k) \neq 0（其中 AkA_k 为 AA 的前 kk 阶顺序主子阵）， 则存在唯一的单位下三角矩阵 LL（对角元全为 11）和对角矩阵 DD（对角元 dk≠0d_k \neq 0），使得：

A=LDLTA = LDL^T

特别地，若 AA 对称正定（SPD），则所有 dk>0d_k > 0。

计算公式：

dk=akk−∑j=1k−1lkj2dj,k=1,2,…,nd_k = a_{kk} - \sum_{j=1}^{k-1} l_{kj}^2 d_j, \quad k=1,2,\dots,n lik=1dk(aik−∑j=1k−1lijlkjdj),i=k+1,…,nl_{ik} = \frac{1}{d_k} \left( a_{ik} - \sum_{j=1}^{k-1} l_{ij} l_{kj} d_j \right), \quad i=k+1,\dots,n

L04 解线性方程组的直接法（误差分析、超定方程组）

向量范数

映射 ∥⋅∥:Rn→R\|\cdot\|:\mathbb{R}^n\to\mathbb{R} 满足：

  1. 正定性：∥x∥≥0\|\mathbf{x}\|\geq 0，且 ∥x∥=0⇔x=0\|\mathbf{x}\|=0 \Leftrightarrow \mathbf{x}=\mathbf{0}
  2. 齐次性：∥αx∥=∣α∣∥x∥\|\alpha\mathbf{x}\| = |\alpha|\|\mathbf{x}\|，∀α∈R\forall \alpha\in\mathbb{R}
  3. 三角不等式：∥x+y∥≤∥x∥+∥y∥\|\mathbf{x}+\mathbf{y}\| \leq \|\mathbf{x}\|+\|\mathbf{y}\|

矩阵范数

映射 ∥⋅∥:Rn×n→R\|\cdot\|:\mathbb{R}^{n\times n}\to\mathbb{R} 满足：

  1. 正定性：∥A∥≥0\|\mathbf{A}\|\geq 0，且 ∥A∥=0⇔A=0\|\mathbf{A}\|=0 \Leftrightarrow \mathbf{A}=\mathbf{0}
  2. 齐次性：∥αA∥=∣α∣∥A∥\|\alpha\mathbf{A}\| = |\alpha|\|\mathbf{A}\|
  3. 三角不等式：∥A+B∥≤∥A∥+∥B∥\|\mathbf{A}+\mathbf{B}\| \leq \|\mathbf{A}\|+\|\mathbf{B}\|
  4. 相容性（次可乘性）：∥AB∥≤∥A∥∥B∥\|\mathbf{AB}\| \leq \|\mathbf{A}\|\|\mathbf{B}\|

向量范数与矩阵范数相容

若对 ∀A∈Rn×n\forall \mathbf{A} \in \mathbf{R}^{n \times n} 与 ∀x∈Rn\forall \mathbf{x} \in \mathbf{R}^n，都有：

∥Ax∥≤∥A∥∥x∥\| \mathbf{A} \mathbf{x} \| \le \| \mathbf{A} \| \| \mathbf{x} \|

则称式中的向量范数和矩阵范数相容。

不是任意向量范数与任意矩阵范数都相容的。但是我们可以做到：

  • 对任意向量范数，构造一个矩阵范数与之相容。即诱导范数。
  • 对任意矩阵范数，构造一个向量范数与之相容。给定矩阵范数 ∥⋅∥\|\cdot\|，定义向量范数如下： ∥x∥:=∥xuT∥\|\mathbf{x}\| := \|\mathbf{x}\mathbf{u}^T\| 其中 u\mathbf{u} 是任意一个固定的非零向量。

诱导范数

由向量范数 ∥⋅∥\|\cdot\| 诱导的：

∥A∥≡max⁡x≠0∥Ax∥∥x∥=max⁡∥x∥=1∥Ax∥\|\mathbf{A}\| \equiv \max_{\mathbf{x}\neq\mathbf{0}} \frac{\|\mathbf{Ax}\|}{\|\mathbf{x}\|} = \max_{\|\mathbf{x}\|=1} \|\mathbf{Ax}\|

是一个矩阵范数。

不是所有矩阵范数都是诱导范数。1,2,∞1, 2, \infty- 向量范数的诱导范数均为对应的矩阵范数。

谱半径

矩阵 A∈Cn×n\mathbf{A} \in \mathbb{C}^{n \times n} 的谱半径定义为：

ρ(A)=max⁡1≤i≤n∣λi∣\rho(\mathbf{A}) = \max_{1\leq i\leq n} |\lambda_i|

定理：对任意矩阵范数，ρ(A)≤∥A∥\rho(\mathbf{A}) \leq \|\mathbf{A}\|

几何意义：衡量特殊方向特征向量的最大拉伸率。

奇异值

矩阵 A∈Cm×n\mathbf{A} \in \mathbb{C}^{m \times n} 的奇异值定义为：

σi(A)=λi(AHA),i=1,2,…,min⁡(m,n)\sigma_i(\mathbf{A}) = \sqrt{\lambda_i(\mathbf{A}^H\mathbf{A})}, \quad i = 1, 2, \dots, \min(m, n)

其中 λi(AHA)\lambda_i(\mathbf{A}^H\mathbf{A}) 是 AHA\mathbf{A}^H\mathbf{A} 的特征值，通常按从大到小排列：

σ1(A)≥σ2(A)≥⋯≥σr(A)>0,σr+1=⋯=0\sigma_1(\mathbf{A}) \ge \sigma_2(\mathbf{A}) \ge \cdots \ge \sigma_r(\mathbf{A}) > 0, \quad \sigma_{r+1} = \cdots = 0

r=rank(A)r = \mathrm{rank}(\mathbf{A})。最大奇异值记为 σmax⁡(A)\sigma_{\max}(\mathbf{A})。

定理：对任意矩阵，最大奇异值等于矩阵的谱范数（即 22-范数）：

σmax⁡(A)=∥A∥2=max⁡∥x∥2=1∥Ax∥2\sigma_{\max}(\mathbf{A}) = \|\mathbf{A}\|_2 = \max_{\|\mathbf{x}\|_2 = 1} \|\mathbf{A}\mathbf{x}\|_2

几何意义：衡量任意方向（单位）向量的最大拉伸率。

我们知道对任意矩阵范数，有 ρ(A)≤∥A∥\rho(\mathbf{A}) \le \|\mathbf{A}\|。故：

ρ(A)≤σmax⁡(A)\rho(\mathbf{A}) \le\sigma_{\max}(\mathbf{A})

从几何意义上看显然奇异值寻找拉伸率的线性空间更大。

条件数

在矩阵范数 ∥⋅∥\|\cdot\| 下，非奇异方阵 A∈Rn×n\mathbf{A} \in \mathbf{R}^{n \times n} 的条件数为：

Cond(A)≡∥A∥∥A−1∥\text{Cond}(\mathbf{A}) \equiv \|\mathbf{A}\|\|\mathbf{A}^{-1}\|

特别地，Cond(A)p=∥A∥p∥A−1∥p\text{Cond}(\mathbf{A})_p = \|\mathbf{A}\|_p\|\mathbf{A}^{-1}\|_p（p=1,2,∞p=1,2,\infty）

  1. 下界： Cond(A)≥1\text{Cond}(\mathbf{A}) \geq 1
  2. 正交矩阵： 若 A\mathbf{A} 正交，则其谱范数下的条件数： Cond(A)2=1\text{Cond}(\mathbf{A})_2 = 1
  3. 齐次性： ∀α≠0\forall \alpha\neq 0 有： Cond(αA)=Cond(A)\text{Cond}(\alpha\mathbf{A}) = \text{Cond}(\mathbf{A})
  4. 谱条件数： Cond(A)2=λmax⁡(ATA)λmin⁡(ATA)=σmax⁡σmin⁡\text{Cond}(\mathbf{A})_2 = \sqrt{\frac{\lambda_{\max}(\mathbf{A}^T\mathbf{A})}{\lambda_{\min}(\mathbf{A}^T\mathbf{A})}} = \frac{\sigma_{\max}}{\sigma_{\min}}
  5. 等价性：Rn×n\mathbb{R}^{n \times n} 上的条件数都是等价的

条件数的几何意义

条件数的倒数 == 矩阵到最近奇异矩阵的相对距离

假设 A∈Rn×n\mathbf{A} \in \mathbb{R}^{n \times n} 可逆， 矩阵 δA∈Rn×n\delta\mathbf{A} \in \mathbb{R}^{n \times n} 使得 A+δA\mathbf{A} + \delta\mathbf{A} 奇异。那么：

∥δA∥∥A∥≥1Cond(A)\frac{\|\delta\mathbf{A}\|}{\|\mathbf{A}\|} \geq \frac{1}{\text{Cond}(\mathbf{A})}

并且如果矩阵范数 ∥⋅∥\|\cdot\| 由向量范数诱导定义，则存在矩阵 δA\delta\mathbf{A} 使得上述不等式取等号：

min⁡A+δA∈S∥δA∥∥A∥=1Cond(A)\min_{\mathbf{A}+\delta\mathbf{A}\in\mathcal{S}} \frac{\|\delta\mathbf{A}\|}{\|\mathbf{A}\|} = \frac{1}{\text{Cond}(\mathbf{A})}

其中 S={M∣det⁡(M)=0}\mathcal{S} = \{\mathbf{M} \mid \det(\mathbf{M}) = 0\} 为奇异矩阵集合。

误差分析

右端项扰动分析：

扰动方程：A(x+δx)=b+δb\mathbf{A}(\mathbf{x}+\delta\mathbf{x}) = \mathbf{b}+\delta\mathbf{b}

则有误差估计：

∥δx∥∥x∥≤∥A∥∥A−1∥∥δb∥∥b∥\frac{\|\delta\mathbf{x}\|}{\|\mathbf{x}\|} \leq \|\mathbf{A}\|\|\mathbf{A}^{-1}\| \frac{\|\delta\mathbf{b}\|}{\|\mathbf{b}\|}

也即解的相对误差不超过右端项相对误差的 ∥A∥∥A−1∥\|\mathbf{A}\|\|\mathbf{A}^{-1}\| 倍。

系数矩阵扰动分析：

扰动方程：(A+δA)(x+δx)=b(\mathbf{A}+\delta\mathbf{A})(\mathbf{x}+\delta\mathbf{x}) = \mathbf{b}

则有误差估计：

∥δx∥∥x∥≤∥A∥∥A−1∥∥δA∥∥A∥1−∥A∥∥A−1∥∥δA∥∥A∥\frac{\|\delta\mathbf{x}\|}{\|\mathbf{x}\|} \leq \frac{\|\mathbf{A}\|\|\mathbf{A}^{-1}\|\frac{\|\delta\mathbf{A}\|}{\|\mathbf{A}\|}}{1-\|\mathbf{A}\|\|\mathbf{A}^{-1}\|\frac{\|\delta\mathbf{A}\|}{\|\mathbf{A}\|}}

也即解的相对误差关于系数矩阵相对误差的函数，当扰动充分小时近似为 ∥A∥∥A−1∥\|\mathbf{A}\|\|\mathbf{A}^{-1}\| 倍。

综合误差分析：

对于非奇异矩阵 A∈Rn×n\mathbf{A} \in \mathbb{R}^{n \times n} 及其扰动 δA∈Rn×n\delta\mathbf{A} \in \mathbb{R}^{n \times n} 满足

∥A−1∥∥δA∥<1\|\mathbf{A}^{-1}\|\|\delta\mathbf{A}\| < 1

如果 x∈Rn\mathbf{x} \in \mathbb{R}^n 是 Ax=b\mathbf{Ax} = \mathbf{b} 的解，其中 b∈Rn,b≠0\mathbf{b} \in \mathbb{R}^n, \mathbf{b} \neq \mathbf{0}。 考虑扰动 δb∈Rn\delta\mathbf{b} \in \mathbb{R}^n，δx\delta\mathbf{x} 是

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

的解。此时有如下正向先验误差估计：

∥δx∥∥x∥≤Cond(A)1−Cond(A)∥δA∥∥A∥(∥δA∥∥A∥+∥δb∥∥b∥)\frac{\|\delta\mathbf{x}\|}{\|\mathbf{x}\|} \leq \frac{\text{Cond}(\mathbf{A})}{1-\text{Cond}(\mathbf{A})\frac{\|\delta\mathbf{A}\|}{\|\mathbf{A}\|}} \left(\frac{\|\delta\mathbf{A}\|}{\|\mathbf{A}\|} + \frac{\|\delta\mathbf{b}\|}{\|\mathbf{b}\|}\right)

L05 解线性方程组的迭代法（Jacobi，Gauss-Seidal，松弛法，收敛性）

统一迭代格式

三种方法可统一写成：

x(k+1)=Mx(k)+g\mathbf{x}^{(k+1)} = \mathbf{M}\mathbf{x}^{(k)} + \mathbf{g}
方法迭代矩阵 M\mathbf{M}向量 g\mathbf{g}
JacobiMJ=D−1(D−A)=I−D−1A\mathbf{M}_J = \mathbf{D}^{-1}(\mathbf{D}-\mathbf{A}) = \mathbf{I} - \mathbf{D}^{-1}\mathbf{A}gJ=D−1b\mathbf{g}_J = \mathbf{D}^{-1}\mathbf{b}
Gauss-SeidelMGS=(D−L)−1U\mathbf{M}_{GS} = (\mathbf{D}-\mathbf{L})^{-1}\mathbf{U}gGS=(D−L)−1b\mathbf{g}_{GS} = (\mathbf{D}-\mathbf{L})^{-1}\mathbf{b}
SORMSOR=(D−ωL)−1[(1−ω)D+ωU]\mathbf{M}_{SOR} = (\mathbf{D}-\omega\mathbf{L})^{-1}[(1-\omega)\mathbf{D} + \omega\mathbf{U}]gSOR=ω(D−ωL)−1b\mathbf{g}_{SOR} = \omega(\mathbf{D}-\omega\mathbf{L})^{-1}\mathbf{b}

Jacobi 迭代格式

将系数矩阵 A\mathbf{A} 分解为：

A=D−(D−A)=D−B\mathbf{A} = \mathbf{D} - (\mathbf{D} - \mathbf{A}) = \mathbf{D} - \mathbf{B}

其中 D=diag(a11,a22,…,ann)\mathbf{D} = \text{diag}(a_{11}, a_{22}, \dots, a_{nn}) 为对角矩阵，B=D−A\mathbf{B} = \mathbf{D} - \mathbf{A}。

原方程 Ax=b\mathbf{A}\mathbf{x} = \mathbf{b} 改写为：

Dx=Bx+b\mathbf{D}\mathbf{x} = \mathbf{B}\mathbf{x} + \mathbf{b}

若 aii≠0a_{ii} \neq 0（对所有 ii），则 D\mathbf{D} 可逆，得到 Jacobi 迭代格式：

x(k+1)=D−1(D−A)x(k)+D−1b=Bx(k)+g\mathbf{x}^{(k+1)} = \mathbf{D}^{-1}(\mathbf{D}-\mathbf{A})\mathbf{x}^{(k)} + \mathbf{D}^{-1}\mathbf{b} = \mathbf{B}\mathbf{x}^{(k)} + \mathbf{g}

Gauss-Seidal 迭代格式

将 A\mathbf{A} 分裂为：

A=D−L−U\mathbf{A} = \mathbf{D} - \mathbf{L} - \mathbf{U}

其中 L\mathbf{L} 为严格下三角矩阵（对角线为零），U\mathbf{U} 为严格上三角矩阵。

迭代格式变为：

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

即 Gauss-Seidal 迭代格式：

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}

松弛法迭代格式

Gauss-Seidel 的更新可看作在当前解 x(k)\mathbf{x}^{(k)} 上加上一个修正量：

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

其中：

Δx=D−1[Lx(k+1)+Ux(k)+b−Dx(k)]\quad \Delta\mathbf{x} = \mathbf{D}^{-1}[\mathbf{L}\mathbf{x}^{(k+1)} + \mathbf{U}\mathbf{x}^{(k)} + \mathbf{b} - \mathbf{D}\mathbf{x}^{(k)}]

引入松弛因子 ω\omega 对修正量进行加权：

x(k+1)=x(k)+ωΔx\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \omega \Delta\mathbf{x}

整理得 SOR 迭代格式：

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

收敛的充要条件

迭代格式对任意初始向量 x(0)\mathbf{x}^{(0)} 都收敛的充分必要条件是：

ρ(M)<1\rho(\mathbf{M}) < 1

迭代格式对任意初始向量 x(0)\mathbf{x}^{(0)} 都收敛的一个充分条件是： 存在某个矩阵范数 ∥⋅∥\|\cdot\| 使得

∥M∥<1\|\mathbf{M}\| < 1

对角占优与严格对角占优

矩阵 A\mathbf{A} 称为对角占优，若：

∣aii∣≥∑j≠i∣aij∣,∀i=1,…,n|a_{ii}| \geq \sum_{j\neq i}|a_{ij}|, \quad \forall i=1,\dots,n

且至少对一个 ii 严格不等式成立。若对所有 ii 都严格成立，则称为严格对角占优。

严格对角占优矩阵的收敛性

若迭代格式矩阵 A\mathbf{A} 严格对角占优，则

  • A\mathbf{A} 非奇异
  • Jacobi 迭代法对任意初始向量收敛
  • Gauss-Seidel 迭代法对任意初始向量收敛

L06 解线性方程组的迭代法（二次函数极值，最速下降法，共轭梯度法）

求解问题

本讲聚焦于 A∈Rn×n\mathbf{A} \in \mathbb{R}^{n \times n} 对称正定的情形。 此时求解线性方程组 Ax=b\mathbf{A}\mathbf{x} = \mathbf{b} 等价于求解如下二次函数的极小值问题：

ϕ(x)=12x⊤Ax−b⊤x\phi(\mathbf{x}) = \frac{1}{2} \mathbf{x}^\top \mathbf{A} \mathbf{x} - \mathbf{b}^\top \mathbf{x}

对于一般的方程 Ax=b\mathbf{A} \mathbf{x} = \mathbf{b}， 转化为 A⊤Ax=A⊤b\mathbf{A}^\top \mathbf{A} \mathbf{x} = \mathbf{A}^\top \mathbf{b} 即可同理求解。

最速下降法

每一步沿当前点处函数下降最快的方向（即负梯度方向）进行一维搜索。

初始化：给定初值 x(0)\mathbf{x}^{(0)}，计算梯度 g(0)=Ax(0)−b\mathbf{g}^{(0)} = \mathbf{A} \mathbf{x}^{(0)} - \mathbf{b}

循环：

  • 计算 t=Ag(k)\mathbf{t} = \mathbf{A} \mathbf{g}^{(k)}
  • 计算步长： α(k)=⟨g(k),g(k)⟩⟨g(k),t⟩\alpha^{(k)} = \dfrac{ \langle \mathbf{g}^{(k)}, \mathbf{g}^{(k)} \rangle } { \langle \mathbf{g}^{(k)}, \mathbf{t} \rangle }
  • 更新解：x(k+1)=x(k)−α(k)g(k)\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} - \alpha^{(k)} \mathbf{g}^{(k)}
  • 更新梯度：g(k+1)=g(k)−α(k)t\mathbf{g}^{(k+1)} = \mathbf{g}^{(k)} - \alpha^{(k)} \mathbf{t}
  • 若 ∥r(k+1)∥2=∥g(k+1)∥2<ε\|\mathbf{r}^{(k+1)}\|_2 = \|\mathbf{g}^{(k+1)}\|_2 < \varepsilon，终止；否则继续循环

共轭梯度法

不再使用负梯度方向，而是构造一组 A\mathbf{A}-共轭搜索方向 {d(k)}\{\mathbf{d}^{(k)}\}，使得在 nn 步内精确收敛，避免震荡。

初始化： 给定初值 x(0)\mathbf{x}^{(0)}， 计算梯度 g(0)=Ax(0)−b\mathbf{g}^{(0)} = \mathbf{A}\mathbf{x}^{(0)} - \mathbf{b}， 设搜索方向 d(0)=−g(0)\mathbf{d}^{(0)} = -\mathbf{g}^{(0)}

循环：

  • 计算 t=Ad(k)\mathbf{t} = \mathbf{A}\mathbf{d}^{(k)}
  • 计算 s=⟨d(k),t⟩s = \langle \mathbf{d}^{(k)}, \mathbf{t} \rangle（即 ⟨d(k),d(k)⟩A\langle \mathbf{d}^{(k)}, \mathbf{d}^{(k)} \rangle _\mathbf{A}）
  • 计算步长： α(k)=−⟨g(k),d(k)⟩s\alpha^{(k)} = -\dfrac{ \langle \mathbf{g}^{(k)}, \mathbf{d}^{(k)} \rangle }{ s }
  • 更新解：x(k+1)=x(k)+α(k)d(k)\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \alpha^{(k)} \mathbf{d}^{(k)}
  • 更新梯度：g(k+1)=g(k)+α(k)t\mathbf{g}^{(k+1)} = \mathbf{g}^{(k)} + \alpha^{(k)} \mathbf{t}
  • 若 ∥r(k+1)∥2=∥g(k+1)∥2<ε\|\mathbf{r}^{(k+1)}\|_2 = \|\mathbf{g}^{(k+1)}\|_2 < \varepsilon，终止
  • 计算 β(k)\beta^{(k)}：选择等价形式中的一种，一般采用 FR 形式。 β(k)=∥g(k+1)∥22∥g(k)∥22\beta^{(k)} = \dfrac{ \|\mathbf{g}^{(k+1)}\|_2^2 }{ \|\mathbf{g}^{(k)}\|_2^2 }
  • 更新搜索方向：d(k+1)=−g(k+1)+β(k)d(k)\mathbf{d}^{(k+1)} = -\mathbf{g}^{(k+1)} + \beta^{(k)} \mathbf{d}^{(k)}

L07 特征值和特征向量的计算（幂法，幂法加速，Jacobi 法）

方法对比与适用场景

方法适用矩阵目标特征值计算复杂度
幂法唯一主特征值最大模特征值O(n2)O(n^2)
位移幂法唯一主特征值离 μ\mu 最远的特征值O(n2)O(n^2)
降阶幂法对称矩阵前 kk 个特征值O(kn3)O(kn^3)
反幂法非奇异矩阵最小模特征值O(n3)O(n^3)（LU分解）+ O(n2)O(n^2)
位移反幂法非奇异矩阵离 μ\mu 最近的特征值O(n3)O(n^3)（LU分解）+ O(n2)O(n^2)
Jacobi 法实对称矩阵所有特征值未优化 O(n4)O(n^4)

幂法

设 A∈Rn×nA \in \mathbb{R}^{n \times n} 有唯一主特征值，即：

∣λ1∣>∣λ2∣≥∣λ3∣≥⋯≥∣λn∣|\lambda_1| > |\lambda_2| \ge |\lambda_3| \ge \cdots \ge |\lambda_n|

其中 λ1\lambda_1 为实数且几何重数为 11。初始向量 x(0)∈Rnx^{(0)} \in \mathbb{R}^n 在 λ1\lambda_1 的特征方向上的投影非零。

为避免溢出问题，每步用 ∞\infty-范数进行归一化：

{y(k)=x(k)∥x(k)∥∞,x(k+1)=Ay(k),k=0,1,2,…\begin{cases} y^{(k)} = \dfrac{x^{(k)}}{\|x^{(k)}\|_\infty}, \\[1em] x^{(k+1)} = A y^{(k)} \end{cases}, \quad k = 0,1,2,\dots

收敛性质：

  • y(k)y^{(k)} 收敛到 λ1\lambda_1 的一个特征向量。
  • 主特征值估计为： λ1(k)=(x(k+1))j,其中 j=arg⁡max⁡i∣(x(k))i∣\lambda_1^{(k)} = (x^{(k+1)})_j, \quad \text{其中 } j = \arg\max_i |(x^{(k)})_i|

位移幂法

注意到幂法收敛速度为：

(x(k+1))j(x(k))j−λ1=O ⁣(∣λ2λ1∣k)\frac{(x^{(k+1)})_j}{(x^{(k)})_j} - \lambda_1 = O\!\left( \left|\dfrac{\lambda_2}{\lambda_1}\right|^k \right)

利用特征值平移性质： 若 Aui=λiuiAu_i = \lambda_i u_i，则 (A−μI)ui=(λi−μ)ui(A - \mu I)u_i = (\lambda_i - \mu)u_i， 也即

  • λi−μ\lambda_i - \mu 为 A−μIA - \mu I 的特征值
  • uiu_i 仍为 A−μIA - \mu I 的特征向量

从而可以选择位移 μ\mu 使收敛更快：

∣λ2−μλ1−μ∣<∣λ2λ1∣\left| \frac{\lambda_2 - \mu}{\lambda_1 - \mu} \right| < \left| \frac{\lambda_2}{\lambda_1} \right|

特殊地，对称正定矩阵特征值均为实数，故最优位移为：

μ=12(λ2+λn)\mu = \frac{1}{2}(\lambda_2 + \lambda_n)

降阶幂法

假设对称矩阵 AA，则其特征向量相互正交。

已知主特征值 λ1\lambda_1 和单位主特征向量 u1u_1，构造：

A(2)=A−λ1u1u1Tu1Tu1A^{(2)} = A - \lambda_1 \frac{u_1 u_1^T}{u_1^T u_1}

则：

  • A(2)u1=0A^{(2)} u_1 = 0
  • A(2)ui=λiuiA^{(2)} u_i = \lambda_i u_i（i=2,3,…,ni=2,3,\dots,n）

即 A(2)A^{(2)} 的特征值为 0,λ2,…,λn0, \lambda_2, \dots, \lambda_n，可对 A(2)A^{(2)} 再次应用幂法求 λ2,u2\lambda_2, u_2。 同理可求更多特征对。

实际上 A(2)A^{(2)} 消除了主特征方向的所有贡献，从而让次特征方向变为 A(2)A^{(2)} 的主特征方向。

反幂法

适用于非奇异矩阵。由特征值倒数性质：若 Aui=λiuiAu_i = \lambda_i u_i，则 A−1ui=λi−1uiA^{-1}u_i = \lambda_i^{-1} u_i。

因此，对 A−1A^{-1} 应用幂法即得 ∣λn∣|\lambda_n| 最小的特征值。

为避免显式求逆，可以将 x(k)=A−1x(k−1)x^{(k)} = A^{-1} x^{(k-1)} 改为解线性方程组：

Ax(k)=x(k−1)A x^{(k)} = x^{(k-1)}

位移反幂法

结合位移技术，我们可以求已知近似特征值 λ∗\lambda^* 附近的精确特征值。

首先构造 B=A−λ∗IB = A - \lambda^* I，然后对 BB 应用反幂法求模最小特征值 μ\mu，则待求特征值即为：

λ=μ+λ∗\lambda = \mu + \lambda^*

由于 ∣λ−λ∗∣≪∣λi−λ∗∣|\lambda - \lambda^*| \ll |\lambda_i - \lambda^*|，收敛极快。

Jacobi 算法

对实对称矩阵 AA，存在正交矩阵 QQ 使 QTAQ=DQ^T A Q = D 为对角阵。Jacobi 法通过一系列 2D 正交变换（Givens 旋转）逐步消除非对角元。

算法流程：

  1. 初始化 V=InV = I_n，用于累积正交变换矩阵
  2. 重复以下步骤，直至所有非对角元绝对值小于给定阈值 ε\varepsilon：
    • 选择当前矩阵中绝对值最大的非对角元 ∣apq∣=max⁡i≠j∣aij∣|a_{pq}| = \max_{i \neq j} |a_{ij}|
    • 计算旋转参数：
    τ=app−aqq2apq,t=−sign(τ)∣τ∣+1+τ2,c=11+t2,s=c⋅t\tau = \frac{a_{pp} - a_{qq}}{2a_{pq}}, \quad t = -\frac{\text{sign}(\tau)}{|\tau| + \sqrt{1+\tau^2}}, \quad c = \frac{1}{\sqrt{1+t^2}}, \quad s = c \cdot t
    • 更新矩阵 A←QpqTAQpqA \leftarrow Q_{pq}^T A Q_{pq}：
      • 更新四个受影响的对角/非对角元： app←c2app+s2aqq−2sc apqaqq←s2app+c2aqq+2sc apqapq←0,aqp←0\begin{aligned} a_{pp} &\leftarrow c^2 a_{pp} + s^2 a_{qq} - 2sc\,a_{pq} \\ a_{qq} &\leftarrow s^2 a_{pp} + c^2 a_{qq} + 2sc\,a_{pq} \\ a_{pq} &\leftarrow 0, \quad a_{qp} \leftarrow 0 \end{aligned}
      • 对 i≠p,qi \neq p, q，更新第 ii 行和第 ii 列： aip←c aip−s aiqaiq←s aip+c aiq\begin{aligned} a_{ip} &\leftarrow c\,a_{ip} - s\,a_{iq} \\ a_{iq} &\leftarrow s\,a_{ip} + c\,a_{iq} \end{aligned} 实际操作中需用临时变量保存旧值。
    • 累积正交矩阵 V←VQpqV \leftarrow V Q_{pq}：
    vip←c vip−s viqviq←s vip+c viq\begin{aligned} v_{ip} &\leftarrow c\,v_{ip} - s\,v_{iq} \\ v_{iq} &\leftarrow s\,v_{ip} + c\,v_{iq} \end{aligned}
  3. 循环结束后，对角元 aiia_{ii} 即为特征值 λi\lambda_i，VV 的第 ii 列即为对应的特征向量

L08 特征值和特征向量的计算（QR 分解，SVD 分解）

Householder 变换

对单位向量 w∈Rnw \in \mathbb{R}^n，Householder 矩阵定义为：

H=I−2wwTH = I - 2ww^T

性质：

  • 对称：HT=HH^T = H
  • 正交：HTH=I⇒H2=IH^T H = I \Rightarrow H^2 = I

几何意义：将向量 xx 关于与 ww 正交的超平面 {y∣wTy=0}\{y \mid w^T y = 0\} 镜像反射。

定理：

对非零向量 x,y∈Rnx, y \in \mathbb{R}^n 且 ∥x∥2=∥y∥2\|x\|_2 = \|y\|_2，存在 Householder 矩阵 HH，使得：

Hx=yHx = y

Householder QR 分解

输入：A∈Rn×nA \in \mathbb{R}^{n \times n}

输出：正交矩阵 QQ 和上三角矩阵 RR，使得 A=QRA = QR

  1. 初始化 A(1)=AA^{(1)} = A，Q(1)=InQ^{(1)} = I_n
  2. 对 k=1,2,…,n−1k = 1, 2, \dots, n-1：
    1. 取 x(k)=A(k)[k:n,k]∈R(n−k+1)x^{(k)} = A^{(k)}[k:n, k] \in \mathbb{R}^{(n-k+1)}
    2. 构造 Householder 矩阵 H(k)∈R(n−k+1)×(n−k+1)H^{(k)} \in \mathbb{R}^{(n-k+1) \times (n-k+1)} 使 H(k)x(k)=−sign(x1(k))∥x(k)∥2e1H^{(k)} x^{(k)} = -\text{sign}(x_1^{(k)})\|x^{(k)}\|_2 e_1
    3. 令 H~(k)=[Ik−100H(k)]\tilde{H}^{(k)} = \begin{bmatrix} I_{k-1} & 0 \\ 0 & H^{(k)} \end{bmatrix}，注意其仍为 Householder 矩阵
    4. A(k+1)=H~(k)A(k)A^{(k+1)} = \tilde{H}^{(k)} A^{(k)}
    5. Q(k+1)=Q(k)(H~(k))TQ^{(k+1)} = Q^{(k)} (\tilde{H}^{(k)})^T
  3. 令 Q=Q(n)Q = Q^{(n)}，R=A(n)R = A^{(n)}，输出 Q,RQ, R 即为所求

Givens 旋转变换

Givens 旋转用于有选择地消去矩阵中的特定元素，每次只影响两行（列）。

具体同 Jacobi 算法。

Givens QR 分解

输入：A∈Rn×nA \in \mathbb{R}^{n \times n}

输出：正交矩阵 QQ 和上三角矩阵 RR，使得 A=QRA = QR

  1. 初始化 A(1)=AA^{(1)} = A，Q(1)=InQ^{(1)} = I_n
  2. 对 k=1,2,…,n−1k = 1, 2, \dots, n-1， 对 i=k+1,…,ni = k+1, \dots, n：
    1. 取 a=A(k)[k,k]a = A^{(k)}[k, k]，b=A(k)[i,k]b = A^{(k)}[i, k]
    2. 计算 c=aa2+b2c = \frac{a}{\sqrt{a^2+b^2}}，s=ba2+b2s = \frac{b}{\sqrt{a^2+b^2}}
    3. 构造 Givens 矩阵 Gik∈Rn×nG_{ik} \in \mathbb{R}^{n \times n}： Gik=[1⋱c⋯s⋮⋱⋮−s⋯c1]G_{ik} = \begin{bmatrix} 1 & & & & & \\ & \ddots & & & & \\ & & c & \cdots & s & \\ & & \vdots & \ddots & \vdots & \\ & & -s & \cdots & c & \\ & & & & & 1 \end{bmatrix} 其中非零元位于 (k,k),(k,i),(i,k),(i,i)(k,k), (k,i), (i,k), (i,i) 位置
    4. A(k+1)=GikTA(k)A^{(k+1)} = G_{ik}^T A^{(k)}
    5. Q(k+1)=Q(k)GikQ^{(k+1)} = Q^{(k)} G_{ik}
  3. 令 Q=Q(n)Q = Q^{(n)}，R=A(n)R = A^{(n)}，输出 Q,RQ, R 即为所求

两种 QR 分解对比

方法每次操作适用场景
Householder消去一整列稠密矩阵
Givens消去单个元素稀疏矩阵、并行计算

QR 算法

定义拟上三角矩阵为对角块为 1 阶或 2 阶的分块上三角矩阵。

实 Schur 分解：

对任意 A∈Rn×nA \in \mathbb{R}^{n \times n}，存在正交矩阵 Q∈Rn×nQ \in \mathbb{R}^{n \times n}，使：

QTAQ=SQ^T A Q = S

其中 SS 为拟上三角矩阵，且：

  • 1 阶对角块对应一个实特征值
  • 2 阶对角块对应一对共轭复特征值

迭代步骤：

  1. QR 分解：A(k)=Q(k)R(k)A^{(k)} = Q^{(k)} R^{(k)}
  2. 更新：A(k+1)=R(k)Q(k)=(Q(k))TA(k)Q(k)A^{(k+1)} = R^{(k)} Q^{(k)} = (Q^{(k)})^T A^{(k)} Q^{(k)}

显然 A(k+1)A^{(k+1)} 与 A(k)A^{(k)} 正交相似，故特征值不变。

收敛性：

设 AA 是 n×nn \times n 实矩阵，特征值满足：

∣λ1∣≥∣λ2∣≥⋯≥∣λn∣|\lambda_1| \ge |\lambda_2| \ge \cdots \ge |\lambda_n|

且等号仅出现在共轭复特征值对（即 λ=a±bi,b≠0\lambda = a \pm bi, b \ne 0）的情形。

则 QR 迭代产生的 A(k)A^{(k)} 收敛到拟上三角矩阵（实 Schur 标准型）。