数值分析 笔记 3

线性方程组

线性代数相应内容复习

对于一个线性方程组，当方程数和变元数不等时，其解的数量通常是无穷或零，因此在本课程中基本不考虑

相应的，我们考虑实方阵A∈Rn×n\boldsymbol{A} \in \mathbb{R}^{n \times n}作为方程系数的情况，且为了方程有解，我们考虑非奇异矩阵（即可逆矩阵）

几种特殊矩阵：

  1. 对称正定矩阵：A=AT\boldsymbol{A} = \boldsymbol{A}^{T}，∀x≠0,xTAx>0\forall \boldsymbol{x} \neq \boldsymbol{0}, \quad \boldsymbol{x}^{T}\boldsymbol{Ax} > 0
  2. 正交矩阵：AAT=I\boldsymbol{AA}^{T} = \boldsymbol{I}

范数

向量的范数

一般的，我们定义满足如下条件的算子∣∣⋅∣∣||\cdot||为范数：

  1. 正定性：∀x∈Rn,∣∣x∣∣≥0\forall \boldsymbol{x} \in \mathbb{R}^{n}, ||\boldsymbol{x}|| \geq 0，等号成立当且仅当x=0\boldsymbol{x} = \boldsymbol{0}
  2. 正齐次性：∀α∈R,∣∣αx∣∣=∣α∣⋅∣∣x∣∣\forall \alpha \in \mathbb{R}, ||\alpha \boldsymbol{x}|| = |\alpha|\cdot||\boldsymbol{x}||
  3. 三角不等式 ∀x,y∈Rn,∣∣x+y∣∣≤∣∣x∣∣+∣∣y∣∣\forall \boldsymbol{x}, \boldsymbol{y} \in \mathbb{R}^{n}, ||\boldsymbol{x} + \boldsymbol{y}|| \leq ||\boldsymbol{x}|| + ||\boldsymbol{y}||

范数一定是以向量元素为变元的nn元连续函数

p-范数和内积范数是两种常见的范数：

  1. p-范数:

    ∣∣x∣∣p=(∑i=1n∣xi∣p)1p||\boldsymbol{x}||_{p} = \Big(\sum\limits_{i=1}^{n}|x_{i}|^{p} \Big)^{\frac{1}{p}}

  2. 内积范数（也就是2-范数）：

    ∣∣x∣∣=⟨x,x⟩=xTx||\boldsymbol{x}|| = \sqrt{\langle \boldsymbol{x}, \boldsymbol{x} \rangle} = \sqrt{\boldsymbol{x}^{T}\boldsymbol{x}}

对于p-范数，通常使用几种特例，包括1-范数，2-范数与无穷范数（也可以分别成为曼哈顿范数、欧式范数、最大范数），公式分别为：

∣∣x∣∣1=∑i=1n∣xi∣∣∣x∣∣2=∑i=1nxi2∣∣x∣∣∞=max⁡1≤i≤n∣xi∣\begin{align*} ||\boldsymbol{x}||_{1} &= \sum\limits_{i=1}^{n} |x_{i}| \\ ||\boldsymbol{x}||_{2} &= \sqrt{\sum\limits_{i=1}^{n}x_{i}^{2}} \\ ||\boldsymbol{x}||_{\infty} &= \max\limits_{1\leq i \leq n} |x_{i}| \end{align*}

针对范数有如下定理：

  1. 等价性：对于任意两种范数s,ts, t，存在常数c1,c2c_{1}, c_{2}，使得：

    c1∣∣x∣∣s≤∣∣x∣∣t≤c2∣∣x∣∣sc_{1}||\boldsymbol{x}||_{s} \leq ||\boldsymbol{x}||_{t} \leq c_{2}||\boldsymbol{x}||_{s}

  2. 极限性：对于向量序列{x(k)}\{\boldsymbol{x}^{(k)}\}对于任意一种范数，有：

    lim⁡k→∞x(k)=x∗⇔lim⁡k→∞∣∣x−x∗∣∣=0\lim_{k\to \infty}\boldsymbol{x}^{(k)} = \boldsymbol{x}^{*} \Leftrightarrow \lim_{k\to\infty}||\boldsymbol{x} - \boldsymbol{x}^{*}|| = 0

    可以先由无穷范数证明，再使用范数的等价性扩展到其他范数

矩阵范数

在向量范数的定义上进行扩张，增加两个条件：

  1. ∣∣AB∣∣≤∣∣A∣∣⋅∣∣B∣∣||\boldsymbol{A}\boldsymbol{B}|| \leq ||\boldsymbol{A}||\cdot||\boldsymbol{B}||
  2. 相容性：∣∣Ax∣∣≤∣∣A∣∣⋅∣∣x∣∣||\boldsymbol{A}\boldsymbol{x}|| \leq ||\boldsymbol{A}||\cdot||\boldsymbol{x}||

我们还可以使用向量的范数∣∣⋅∣∣v||\cdot||_{v}来定义矩阵的范数，这种定义称为向量诱导范数

∣∣A∣∣v=max⁡x≠0∣∣Ax∣∣v∣∣x∣∣v||\boldsymbol{A}||_{v} = \max\limits_{\boldsymbol{x} \neq \boldsymbol{0}} \frac{||\boldsymbol{Ax}||_{v}}{||\boldsymbol{x}||_{v}}

可以证明，向量诱导范数符合前述的定义

利用向量的p-范数进行诱导，我们可以得到几种矩阵的p-范数：

∣∣A∣∣1=max⁡1≤j≤n∑i=1n∣aij∣∣∣A∣∣∞=max⁡1≤i≤n∑j=1n∣aij∣∣∣A∣∣2=λmax⁡(ATA)\begin{align*} ||\boldsymbol{A}||_{1} &= \max\limits_{1\leq j\leq n}\sum\limits_{i=1}^{n}|a_{ij}| \\ ||\boldsymbol{A}||_{\infty} &= \max\limits_{1\leq i\leq n}\sum\limits_{j=1}^{n}|a_{ij}| \\ ||\boldsymbol{A}||_{2} &= \sqrt{\lambda_{\max}(\boldsymbol{A}^{T}\boldsymbol{A})} \\ \end{align*}

此外还有常用的Frobenius范数：

∣∣A∣∣F=∑i=1n∑j=1naij2=tr(ATA)||\boldsymbol{A}||_{F} = \sqrt{\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}a_{ij}^{2}} = \sqrt{\mathrm{tr}(\boldsymbol{A}^{T}\boldsymbol{A})}

矩阵条件数

下面讨论解线性方程组问题的敏感性

考虑右端发生小扰动，形式为：

Ax=b⇒A(x+Δx)=b+Δb\boldsymbol{Ax}=\boldsymbol{b} \Rightarrow \boldsymbol{A}(\boldsymbol{x} + \Delta\boldsymbol{x}) = \boldsymbol{b} + \Delta\boldsymbol{b}

这种情况下解的扰动是Δx\Delta\boldsymbol{x}，因此根据条件数的定义：

cond=∣∣Δx∣∣/∣∣x∣∣∣∣Δb∣∣/∣∣b∣∣=∣∣Δx∣∣⋅∣∣b∣∣∣∣Δb∣∣⋅∣∣x∣∣\mathrm{cond} = \frac{||\Delta\boldsymbol{x}|| / ||\boldsymbol{x}||}{||\Delta\boldsymbol{b}||/||\boldsymbol{b}||} = \frac{||\Delta\boldsymbol{x}||\cdot||\boldsymbol{b}||}{||\Delta\boldsymbol{b}||\cdot||\boldsymbol{x}||}

根据范数的定义，可以有如下的推导：

∣∣Δx∣∣=∣∣A−1Δb∣∣≤∣∣A−1∣∣⋅∣∣Δb∣∣∣∣b∣∣=∣∣Ax∣∣≤∣∣A∣∣⋅∣∣x∣∣⇒∣∣b∣∣⋅∣∣Δx∣∣≤∣∣A−1∣∣⋅∣∣Δb∣∣⋅∣∣A∣∣⋅∣∣x∣∣\begin{align*} ||\Delta\boldsymbol{x}|| &= ||\boldsymbol{A}^{-1}\Delta\boldsymbol{b}|| \leq ||\boldsymbol{A}^{-1}||\cdot||\Delta\boldsymbol{b}|| \\ ||\boldsymbol{b}|| &= ||\boldsymbol{Ax}|| \leq ||\boldsymbol{A}||\cdot||\boldsymbol{x}|| \\ \Rightarrow\qquad\qquad ||\boldsymbol{b}||\cdot||\Delta\boldsymbol{x}|| &\leq ||\boldsymbol{A}^{-1}||\cdot||\Delta\boldsymbol{b}||\cdot||\boldsymbol{A}||\cdot||\boldsymbol{x}|| \\ \end{align*}

代入有：

cond≤∣∣A∣∣⋅∣∣A−1∣∣\mathrm{cond} \leq ||\boldsymbol{A}||\cdot||\boldsymbol{A}^{-1}||

于是我们定义这个上界为矩阵的条件数，即其代表了解线性方程问题敏感性的上界：

cond(A)=∣∣A∣∣⋅∣∣A−1∣∣\mathrm{cond}(\boldsymbol{A}) = ||\boldsymbol{A}||\cdot||\boldsymbol{A}^{-1}||

可以证明的是，如果发生扰动的并非b\boldsymbol{b}而是A\boldsymbol{A}，也可以近似到上述上界

条件数越大，矩阵越病态，而对于奇异矩阵，我们定义期条件数为正无穷大

例如对于Hilbert矩阵，其可以计算出其条件数随着nn的增大会急剧增大，因此会变得越来越病态

Hn=[112⋯1n1213⋯1n+1⋮⋮⋱⋮1n1n+1⋯12n−1]\boldsymbol{H}_{n} = \begin{bmatrix} 1 & \dfrac{1}{2} & \cdots & \dfrac{1}{n} \\[1em] \dfrac{1}{2} & \dfrac{1}{3} & \cdots & \dfrac{1}{n+1} \\[1em] \vdots & \vdots & \ddots & \vdots \\[1em] \dfrac{1}{n} & \dfrac{1}{n + 1} & \cdots & \dfrac{1}{2n - 1} \end{bmatrix}

一些定理

对于任意向量范数意义下的矩阵范数，有如下定理：

cond(A)=(max⁡x≠0∣∣Ax∣∣∣∣x∣∣) / (min⁡x≠0∣∣Ax∣∣∣∣x∣∣)≥1\mathrm{cond}(\boldsymbol{A}) = \Big(\max\limits_{\boldsymbol{x} \neq \boldsymbol{0}}\frac{||\boldsymbol{Ax}||}{||\boldsymbol{x}||}\Big)\,\Big/\,\Big(\min\limits_{\boldsymbol{x} \neq \boldsymbol{0}}\frac{||\boldsymbol{Ax}||}{||\boldsymbol{x}||}\Big) \geq 1

同时，根据定义很容易推导出矩阵条件数满足的性质：

  1. cond(A−1)=cond(cA)=cond(A)\mathrm{cond}(\boldsymbol{A}^{-1}) = \mathrm{cond}(c\boldsymbol{A}) = \mathrm{cond}(\boldsymbol{A})
  2. cond(I)\mathrm{cond}(\boldsymbol{I}) = 1
  3. 对于一个对角阵D\boldsymbol{D}，有 cond(D)≥max⁡∣dii∣min⁡∣dii∣\mathrm{cond}(\boldsymbol{D})\geq \dfrac{\max|d_{ii}|}{\min|d_{ii}|}，其中等号在p-范数下取到
  4. 正交变换不改变2-范数下的条件数

LU分解

LU分解即为将一个矩阵拆分成一个下三角与一个上三角矩阵的乘积的形式，分为两种：

  • Doolittle分解：LL为单位下三角阵，UU为上三角阵
  • Crout分解：LL为下三角阵，UU为单位上三角阵

本章讨论Doolittle分解

LU分解的存在性有定理：

对于方程Ax=b\boldsymbol{Ax} = \boldsymbol{b}，其中A∈Rn×n\boldsymbol{A} \in \mathbb{R}^{n\times n}，则执行高斯消元过程中的主元 akk(k)≠0a_{kk}^{(k)} \neq 0 等价于 A\boldsymbol{A}有唯一LU分解

其中 akk(k)a_{kk}^{(k)} 代表消除kk步之后的元素

LU分解的计算很简单，不再赘述，其主要用途在解方程

选主元

可以看出，LU分解的计算与主元的选取有关，而这是可以通过对换矩阵来进行的，并且对换操作不会改变方程的解（可以想象若干方程以任意的顺序排列，它的解当然是不会变的）

有定理：

对于矩阵 A∈Rn×n\boldsymbol{A} \in \mathbb{R}^{n\times n}，执行高斯消元过程中的主元 akk(k)a_{kk}^{(k)}均不为零 等价于 A\boldsymbol{A} 的前 n−1n - 1 个顺序主子式均不为零

因此对称正定或负定矩阵一定可以LU分解

如果在高斯消去的过程中，遇到了零主元的情况，可以从该行下方的行中取一行，满足 ajk(k)≠0(j>k)a_{jk}^{(k)} \neq 0\quad (j > k) 与这一行交换，之后就可以正常进行下去了，可以证明的是，如果 A\boldsymbol{A} 非奇异，那么一定可以完成这个交换的操作

同时，可以选择绝对值尽量大的主元，这样可以有效避免误差的扩大传播

部分主元

这种方法是为了满足上述 可以选择绝对值尽量大的主元，这样可以有效避免误差的扩大传播 的操作，具体做法是在每一步消元的时候，将第 iki_k 行与第 kk 行交换，其中，ik=argmaxk≤i≤n∣aik(k)∣i_{k} = \mathop{\mathrm{argmax}}\limits_{k\leq i\leq n} |a_{ik}^{(k)}|

下面推导部分主元的LU分解：

假设 Pk\boldsymbol{P}_{k} 代表第 kk 步消元过程中的交换过程，即交换第 iki_{k} 行与第 kk 行的对换矩阵，则消去过程可以表示为：

Mn−1Pn−1⋯M2P2M1P1A=U\boldsymbol{M}_{n-1}\boldsymbol{P}_{n-1}\cdots\boldsymbol{M}_{2}\boldsymbol{P}_{2}\boldsymbol{M}_{1}\boldsymbol{P}_{1}\boldsymbol{A} = \boldsymbol{U}

由于对换矩阵的逆是自身，因此我们令：

M‾k=(∏i=n−1k+1Pi)Mk(∏i=k+1n−1Pi)\overline{\boldsymbol{M}}_{k} = \Big(\prod_{i=n-1}^{k+1}\boldsymbol{P}_{i}\Big)\boldsymbol{M}_{k}\Big(\prod_{i=k+1}^{n-1}\boldsymbol{P}_{i}\Big)

这样原始消去可以变换为：

(∏i=n−11M‾k)(∏i=n−11Pk)A=U\Big(\prod_{i=n-1}^{1}\overline{\boldsymbol{M}}_{k}\Big)\Big(\prod_{i=n-1}^{1}\boldsymbol{P}_{k}\Big)\boldsymbol{A} = \boldsymbol{U}

记：

L=(∏i=n−11M‾k)−1P=(∏i=n−11Pk)\begin{align*} \boldsymbol{L} &= \Big(\prod_{i=n-1}^{1}\overline{\boldsymbol{M}}_{k}\Big)^{-1} \\ \boldsymbol{P} &= \Big(\prod_{i=n-1}^{1}\boldsymbol{P}_{k}\Big)\end{align*}

则部分主元的LU分解为：

PA=LU\boldsymbol{PA} = \boldsymbol{LU}

实际计算过程中，我们并不需要按照推导式中的方法，求出如此多的 M‾k\overline{\boldsymbol{M}}_{k} ，相反我们只需要正常的进行 LU 分解，但是在需要交换的时候，记录对换矩阵 Pk\boldsymbol{P}_{k} ，这样最终求出的分解结果就是所需的 LU\boldsymbol{LU} ，而 P\boldsymbol{P} 是容易求出的

全选主元

在上述过程中，由于我们只是选取了列中最大的元素，因此这个方法称为部分选主元，相应的，可以在第 kk 步选取右下角(n−k+1)×(n−k+1)(n - k + 1) \times (n - k + 1) 维矩阵中的绝对值最大元素，通过一次行变换与一次列变换将其挪到主元位置，这种方法称为全选主元，得到的分解结果为：

PAQ=LU\boldsymbol{PAQ} = \boldsymbol{LU}

应用

可以方便的通过 LU 分解来进行解方程、矩阵乘的计算，假设我们知道了矩阵的 LU 分解，尽管这个计算的复杂度是O(n3)O(n^{3})的，但是我们可以通过这个操作来将矩阵乘的计算降到 O(n2)O(n^{2})

例如我们需要计算 (A+B−1)x(\boldsymbol{A} + \boldsymbol{B}^{-1})\boldsymbol{x}，我们可以按如下方式计算：

PAx=LAUAx⇒Ax=P−1LAUAxB−1x=UB−1LB−1Px\begin{align*} \boldsymbol{PAx} = \boldsymbol{L}_{A}\boldsymbol{U}_{A}\boldsymbol{x} &\Rightarrow \boldsymbol{Ax} = \boldsymbol{P}^{-1}\boldsymbol{L}_{A}\boldsymbol{U}_{A}\boldsymbol{x} \\ \boldsymbol{B}^{-1}\boldsymbol{x} &= \boldsymbol{U}_{B}^{-1}\boldsymbol{L}_{B}^{-1}\boldsymbol{P}\boldsymbol{x} \end{align*}

其中，计算逆可以转换成解方程，例如：

LB−1Px=y⇒LB y=Px\boldsymbol{L}_{B}^{-1}\boldsymbol{Px} = \boldsymbol{y} \Rightarrow \boldsymbol{L}_{B}\,\boldsymbol{y} = \boldsymbol{Px}

由于LU\boldsymbol{LU}都是三角阵，因此解这个方程是 O(n2)O(n^{2}) 的，如果需要进行大量计算，我们可以对 LU 分解的结果做记忆化，这样可以显著减少后续计算的计算量

稳定性

对于 LU 分解算法，有结论：

∣∣ΔA∣∣∞∣∣A∣∣∞≲nρεmach\frac{||\Delta\boldsymbol{A}||_{\infty}}{||\boldsymbol{A}||_{\infty}} \lesssim n\rho\varepsilon_{\text{mach}}

增长因子 ρ\rho 定义为：

ρ=max⁡∣uij∣max⁡∣aij∣\rho = \frac{\max |u_{ij}|}{\max |a_{ij}|}

Cholesky 分解

对于一个可无需交换主元的对称矩阵，首先将其 LU 分解扩展为 LDU 分解，满足其中的 LU\boldsymbol{LU} 都是单位三角阵，D\boldsymbol{D} 是对角阵，这样的话：

A=LDU=AT=UTDLT\boldsymbol{A} = \boldsymbol{LDU} = \boldsymbol{A}^{T} = \boldsymbol{U}^{T}\boldsymbol{DL}^{T}

因此 L=UT\boldsymbol{L} = \boldsymbol{U^{T}}，即 A=LDLT\boldsymbol{A} = \boldsymbol{LDL}^{T}

就此引入 Cholesky 分解：

对于一个对称正定矩阵 A\boldsymbol{A}，分解 A=LDLT\boldsymbol{A} = \boldsymbol{LDL}^{T} 满足 D\boldsymbol{D} 也是对称正定的，因此我们可以令：

D=diag(d11,d22,…,dnn)\sqrt{\boldsymbol{D}} = \mathrm{diag}(\sqrt{d_{11}}, \sqrt{d_{22}}, \dots, \sqrt{d_{nn}})

这样的话：

A=LDDTLT=L1L1T\boldsymbol{A} = \boldsymbol{L}\sqrt{\boldsymbol{D}}\sqrt{\boldsymbol{D}}^{T}\boldsymbol{L}^{T} = \boldsymbol{L}_{1}\boldsymbol{L}_{1}^{T}

这就是 Cholesky 分解

其具体算法两层遍历，但是只用遍历矩阵的半角（因为是对称的）

1
2
3
4
for j in range(1, n + 1):
   (1)
   for i in range(j + 1, n + 1):
      (2)

其中：

ajj=ajj−∑k=1j−1ajk2aij=1ajj(aij−∑k=1j−1aikajk)\begin{align} a_{jj} &= \sqrt{a_{jj} - \sum\limits_{k=1}^{j-1}a_{jk}^{2}} \\ a_{ij} &= \frac{1}{a_{jj}}\Big(a_{ij} - \sum\limits_{k=1}^{j-1}a_{ik}a_{jk}\Big) \end{align}

可以证明其增长因子 ρ≤1\rho \leq 1

带状矩阵

定义：

若矩阵 A\boldsymbol{A} 满足：

  • ∀i,j,∣i−j∣>β⇒aij=0\forall i, j, |i-j|>\beta \Rightarrow a_{ij} = 0
  • ∃k,ak(k−β)2+ak(k+β)2≠0\exists k, a_{k(k-\beta)}^{^2} + a_{k(k+\beta)}^{2} \neq 0

则称 A\boldsymbol{A} 是半带宽为 β\beta 的带状矩阵

可以证明带状矩阵的 LU 分解结果满足，其非零元仍然在原始带宽内

考虑三对角矩阵，即β=1\beta = 1的带状矩阵，例如：

[12345678910]\begin{bmatrix} 1 & 2 & & \\ 3 & 4 & 5 & \\ & 6 & 7 & 8 \\ & & 9 & 10 \\ \end{bmatrix}

三对角矩阵的 LU 分解算法是 O(n)O(n) 的，令：

A=[b1c1a2b2c2⋱⋱⋱an−1bn−1cn−1anbn]\boldsymbol{A} = \begin{bmatrix} b_{1} & c_{1} & & & \\[1em] a_{2} & b_{2} & c_{2} & & \\[1em] & \ddots & \ddots & \ddots & \\[1em] & & a_{n-1} & b_{n-1} & c_{n-1} \\[1em] & & & a_{n} & b_{n} \end{bmatrix}

则其不选主元的 LU 分解一定是：

L=[1m21⋱⋱mn−11mn1]U=[d1c1d2c2⋱⋱dn−1cn−1dn]\boldsymbol{L} = \begin{bmatrix} 1 & & & & \\[1em] m_{2} & 1 & & & \\[1em] & \ddots & \ddots & & \\[1em] & & m_{n-1} & 1 & \\[1em] & & & m_{n} & 1 \end{bmatrix} \quad \boldsymbol{U} = \begin{bmatrix} d_{1} & c_{1} & & & \\[1em] & d_{2} & c_{2} & & \\[1em] & & \ddots & \ddots & \\[1em] & & & d_{n-1} & c_{n-1} \\[1em] & & & & d_{n} \end{bmatrix}

其中：

d1=b1mi=ai/di−1di=bi−mici−12≤i≤n\begin{align*} d_{1} &= b_{1} \\ m_{i} &= a_{i}/d_{i - 1} \\ d_{i} &= b_{i} - m_{i}c_{i - 1} \\ 2 &\leq i \leq n \end{align*}

这可以很方便由直接 LU 分解得到

而解决三对角矩阵为系数的线性方程组也是很高效的，算法称为追赶法，具体来说，针对上述定义的A\boldsymbol{A}和方程Ax⃗=f⃗\boldsymbol{A}\vec{x}=\vec{f}：

2≤i≤nn−1≥j≥1mi=ai/bi−1bi=bi−mici−1fi=fi−mifi−1xn=fn/bnxj=(fj−cjxj+1)/bj\begin{align*} 2\leq i\leq n&\quad n - 1 \geq j \geq 1 \\ m_{i} &= a_{i} / b_{i - 1} \\ b_{i} &= b_{i} - m_{i}c_{i - 1} \\ f_{i} &= f_{i} - m_{i}f_{i - 1} \\ x_{n} &= f_{n} / b_{n} \\ x_{j} &= (f_{j} - c_{j}x_{j + 1}) / b_{j} \end{align*}

一般带状矩阵的 LU 分解时间复杂度为 O(β2n)O(\beta^{2}n)，空间复杂度为 O(βn)O(\beta n)，这说明 β<<n\beta << n 时甚至是线性的，这是由带状矩阵的稀疏特性决定的，而带状矩阵的逆是稠密的，因此应该尽量避免求逆

可以证明，对称正定矩阵与严格对角占有矩阵是可以不选主元进行稳定 LU 分解的，如果使用部分选主元，LU\boldsymbol{LU}矩阵的带宽不会超过 2β2\beta

许可协议

本文采用 署名-非商业性使用-相同方式共享 4.0 国际 许可协议，转载请注明出处。