数值分析(3) 线性方程组的直接解法

线性方程组与矩阵的基本概念

线性方程组形如下面的样子：

$$\begin{cases}a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1\\ a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2\\ \cdots\\a_{m1}x_1+a_{m2}x_2+\cdots+a_{mn}x_n=b_m\\ \end{cases}$$

或矩阵形式：

$$\mathbf{A}\mathbf{x}=\mathbf{b}$$

• 若$m>n$，称为超定方程组，一般无解，可以考虑最小二乘解
• 若$m<n$，一般有无穷多解，可以考虑约束优化问题

---

现在讨论的情况均假定$m=n$，即要求

$$\mathbf{A}\in \mathbb{R}^{n\times n}$$

$$\mathbf{b}\in\mathbb{R}^n$$

基础矩阵概念复习

• 对称矩阵：$\mathbf{A}^T=\mathbf{A}$
• 对称正定矩阵：$\mathbf{A}^T=\mathbf{A}$且$\forall \mathbf{x}\neq \mathbf{0},\ \mathbf{x}^T\mathbf{A}\mathbf{x}>0$
• 正交矩阵：$\mathbf{A}^T=\mathbf{A}^{-1}$
• 顺序主子阵：左上角的方块；其行列式称为顺序主子式

向量的范数

定义（一般的范数）

对线性空间$\mathcal{V}$，若一个映射$||\cdot||:\ \mathcal{V}\to\mathbb{R}$满足：

• 正定性：$\forall\mathbf{x}\in\mathcal{V},\ ||\mathbf{x}||\geq0$，当且仅当$\mathbf{x}=\mathbf{0}$时取等
• 正齐次性：$||\alpha\mathbf{x}||=|\alpha|\cdot||\mathbf{x}||$
• 三角不等式：$\forall \mathbf{x},\mathbf{y}\in\mathcal{V},\ ||\mathbf{x}+\mathbf{y}||\leq||\mathbf{x}||+||\mathbf{y}||$

则称该映射是一个范数，此空间是一个赋范线性空间

---

常用的向量范数包括：

• $p$-范数：$||\mathbf{x}||_p=(\sum\limits_{i=1}^n|x_i|^p)^{1/p}$
• 内积范数：对一个内积$\left<\cdot,\cdot\right>$，$||\mathbf{x}||=\sqrt{\left<\mathbf{x},\mathbf{x}\right>}$
  • 对常用内积$\left<\mathbf{x},\mathbf{y}\right>=\mathbf{x}^T\mathbf{y}$而言，这等价于2-范数

---

定义（范数的等价性）

若对两种范数$||\cdot||_s$和$||\cdot||_t$，有

$$\exists c_1,c_2>0,\ \forall \mathbf{x}\in\mathcal{V},\ c_1||\mathbf{x}||_s\leq||\mathbf{x}||_t\leq c_2||\mathbf{x}||_s$$

则称两种范数等价

注 范数等价意味着用这两种范数定义的极限、连续等等概念都相同

定理 $\mathbb{R}^n$上的一切范数等价

注 这让我们证明一些极限问题时可以任意选择一种范数

矩阵的范数

定义（矩阵范数）

矩阵范数需要满足向量范数的正定、正齐次、三角不等式三条性质，同时增加乘法限制：

$$\forall \mathbf{A},\mathbf{B}\in\mathbb{R}^{n\times n},\ ||\mathbf{A}\mathbf{B}||\leq||\mathbf{A}||\ ||\mathbf{B}||$$

$$\forall \mathbf{A}\in\mathbb{R}^{n\times n},\ \mathbf{x}\in\mathbb{R}^n,\ ||\mathbf{A}\mathbf{x}||\leq||\mathbf{A}||\ ||\mathbf{x}||$$

注 上述第二条称为相容性条件，它建立了矩阵范数和向量范数的联系

---

矩阵的默认范数为算子范数，它根据某种向量范数$||\cdot||_v$来定义：

$$||\mathbf{A}||_v=\max\limits_{\mathbf{v}\neq\mathbf{0}}\frac{||\mathbf{Av}||_{v}}{||\mathbf{v}||_v}$$

算子范数的含义是矩阵作为一种作用于向量上的算子，对向量的最大拉伸倍数

向量范数使用2-范数时，$||\mathbf{A}||_v$表示椭圆$\mathbf{x}^T\mathbf{A}^{-T}\mathbf{A}^{-1}\mathbf{x}=1$的半长轴长度

根据算子范数的定义，很容易证明算子范数是一种合要求的矩阵范数

---

常用矩阵范数：

• 1-范数：$||\mathbf{A}||_1=\max\limits_j\sum\limits_{i=1}^n|a_{ij}|$
• ∞-范数：$||\mathbf{A}||_\infty=\max\limits_i\sum\limits_{j=1}^n|a_{ij}|$
• 2-范数：矩阵$\mathbf{A}^T\mathbf{A}$的最大特征值的算术平方根
• Frobenius范数：$||\mathbf{A}||_F=\sqrt{\sum\limits_i\sum\limits_ja_{ij}^2}=\sqrt{\mathrm{tr}(\mathbf{A}^T\mathbf{A})}$

---

上述除了Frobenius范数以外都是算子范数，矩阵的$p$-范数表示利用向量$p$-范数得到的算子范数。证明如下：

矩阵的1-范数

$$||\mathbf{A}||_1=\max\limits_{\mathbf{v}\neq\mathbf{0}}\frac{\left|\left|\sum\limits_{j}v_j\mathbf{a}_j\right|\right|_1}{||\mathbf{v}||_1}$$

由于

$$\frac{\left|\left|\sum\limits_{j}v_j\mathbf{a}_j\right|\right|_1}{||\mathbf{v}||_1}\leq\frac{|v_1|\ ||\mathbf{a}_1||_1+\cdots+|v_n|\ ||\mathbf{a}_n||_1}{|v_1|+\cdots+|v_n|}\leq\max_{j}||\mathbf{a}_j||_1$$

且此不等式在$\mathbf{v}=e_{j_m}$，其中$j_m=\mathrm{argmax}_j(||\mathbf{a}_j||_1)$时取等，于是

$$||\mathbf{A}||_1=\max\limits_{\mathbf{v}\neq\mathbf{0}}\frac{\left|\left|\sum\limits_{j}v_j\mathbf{a}_j\right|\right|_1}{||\mathbf{v}||_1}=\max\limits_j||\mathbf{a}_j||_1$$

矩阵的∞-范数

$$||\mathbf A||_\infty=\max_{\mathbf v\neq\mathbf 0}\frac{\max\limits_i\left|\sum\limits_{j}a_{ij}v_j\right|}{\max\limits_j|v_j|}$$

考虑当$i$固定，$\mathbf{v}$可变时，下式的最大值：

$$\frac{\left|\sum\limits_{j}a_{ij}v_j\right|}{\max\limits_j|v_j|}$$

显然，当所有$v_j$取绝对值相同，且符号使分子的各和项同号时取最大，最大值为

$$\sum\limits_j|a_{ij}|$$

于是

$$||\mathbf A||_\infty=\max\limits_i\sum\limits_j|a_{ij}|$$

矩阵的2-范数

取$\frac{||\mathbf{Av}||_2}{||\mathbf{v}||_2}$的最大值等价于取

$$\frac{|\mathbf v^T\mathbf A^T\mathbf{Av}|}{|\mathbf v^T\mathbf v|}$$

的最大值

由于$\mathbf A^T\mathbf A$实对称，可以用其特征向量构成一组标准正交基$\{\mathbf x_1,\cdots,\mathbf x_n\}$，在此基上分解$\mathbf v$，得到

$$\mathbf v=v_1\mathbf x_1+\cdots+v_n\mathbf x_n$$

于是

$$\mathbf A^T\mathbf{Av}=\lambda_1v_1\mathbf x_1+\cdots+\lambda_nv_n\mathbf x_n$$

由于

$$\mathbf x_i^T\mathbf x_j=\begin{cases}0,&i\neq j\\1,&i=j\end{cases}$$

因此

$$\frac{|\mathbf v^T\mathbf A^T\mathbf{Av}|}{|\mathbf v^T\mathbf v|}=\frac{|\lambda_1v_1^2+\cdots+\lambda_nv_n^2|}{|v_1^2+\cdots+v_n^2|}$$

由于$\mathbf A^T\mathbf A$半正定

$$\frac{|\mathbf v^T\mathbf A^T\mathbf{Av}|}{|\mathbf v^T\mathbf v|}=\frac{|\lambda_1v_1^2+\cdots+\lambda_nv_n^2|}{|v_1^2+\cdots+v_n^2|}=\frac{\lambda_1v_1^2+\cdots+\lambda_nv_n^2}{v_1^2+\cdots+v_n^2}\leq \lambda_\mathrm{max}$$

当$\mathbf v$是$\lambda_\mathrm{max}$对应的特征向量时取等

线性方程组求解问题的敏感性

考虑以下问题：

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

当右端发生扰动时：

$$\mathbf{A}(\mathbf{x}+\Delta\mathbf{x})=\mathbf{b}+\Delta\mathbf{b}$$

根据线性性：

$$\mathbf{A}\Delta\mathbf{x}=\Delta\mathbf{b}$$

于是问题的条件数

$${\rm cond}=\frac{||\Delta\mathbf{x}||/||\mathbf{x}||}{||\Delta\mathbf b||/||\mathbf b||}\leq||\mathbf{A}||\ \left|\left|\mathbf A^{-1}\right|\right|$$

定义矩阵的条件数为

$${\rm cond}(\mathbf{A})=||\mathbf{A}||\ \left|\left|\mathbf A^{-1}\right|\right|$$

条件数较大的矩阵对应着很病态的求解问题，称之为病态矩阵

---

定理

$$\mathrm{cond}(\mathbf A)=\frac{\max\limits_{\mathbf v\neq\mathbf 0}\frac{||\mathbf{Av}||}{||\mathbf v||}}{\min\limits_{\mathbf v\neq\mathbf 0}\frac{||\mathbf{Av}||}{||\mathbf v||}}$$

证明是容易的

注 定理说明条件数表示矩阵对单位圆的扭曲程度，它的值为比矩阵作用在单位圆上得到的椭圆的半长轴与半短轴之比

定义奇异矩阵的条件数为$+\infty$

---

定理（矩阵条件数的性质）

1. ${\rm cond}(\mathbf A)\geq1$
2. ${\rm cond}(\mathbf A^{-1})={\rm cond}(\mathbf A)$
3. ${\rm cond}(c\mathbf A)={\rm cond}(\mathbf A),\ c\neq0$
4. 单位阵条件数为$1$
5. 对角阵条件数为对角元的最大绝对值与最小绝对值之比
6. 2-范数下条件数是$\mathbf A^T\mathbf A$的最大特征值与最小特征值之比的算术平方根
7. 2-范数下乘正交矩阵不改变条件数

高斯消去法

消去过程：将增广矩阵的系数矩阵部分变为上三角

回代过程：解上三角方程组

在第$k$步消去前的主元$a_{kk}$不能为$0$，否则无法消去

消去过程可以使用原地工作，不消耗额外的存储空间

计算复杂度为$O(n^3)$

LU分解

以下过程为矩阵的LU分解：

$$\mathbf A =\mathbf {LU}$$

其中$\mathbf L$为单位下三角矩阵，$\mathbf U$为上三角矩阵

---

LU分解可以通过在高斯消去法的同时记录消去的过程而得到，这种方法在线性代数课程中已经有过充分介绍，这里不再赘述

定理（LU分解的条件） 对方阵$\mathbf A\in \mathbb R^{n\times n}$，其存在唯一的LU分解当且仅当执行高斯消去过程中第$k$步消去前的主元$a_{kk}^{(k)}\neq 0,\ \forall k=1,\cdots,n-1$

直接LU分解算法

如果我们不通过高斯消去的思路来进行LU分解，而是直接列出逐元素的表达式（这里以3阶为例）：

$$\begin{bmatrix}1\\l_{21}&1\\l_{31}&l_{32}&1\end{bmatrix}\begin{bmatrix}u_{11}&u_{12}&u_{13}\\&u_{22}&u_{23}\\&&u_{33}\end{bmatrix}=\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{bmatrix}$$

可以列出如下方程：

$$u_{1j}=a_{1j},\ j=1,2,3$$

$$l_{i1}u_{11}=a_{i1},\ i=2,3$$

$$l_{21}u_{1j}+u_{2j}=a_{2j},\ j=2,3$$

$$l_{31}u_{12}+l_{32}u_{22}=a_{32}$$

$$l_{31}u_{13}+l_{32}u_{23}+u_{33}=a_{33}$$

于是可以循着

$$u_{1j}\to l_{i1}\to u_{2j}\to l_{i2}\to u_{3j}\to l_{i3}\to\cdots\to l_{i(n-1)}\to u_{nn}$$

的顺序来解出所有$\mathbf L$和$\mathbf U$的元素，每个元素的求解过程类似于简单的回代法

---

特别地，注意到在上面的方程中，形如

$$a_{ij}-l_{ik}u_{kj}$$

的形状会频繁出现在方程右端，因此可以使用循环更新$\mathbf A$的方式来减少重复计算，具体地，对：

$$\mathbf A=\begin{bmatrix} 2&1&1&2\\ 2&2&2&3\\ 4&2&4&3\\ 0&0&6&-1 \end{bmatrix}$$

首先解出$\mathbf U$的第一行，它就是$\mathbf A$的第一行：

$$\mathbf L=\begin{bmatrix} 1\\ \cdot&1\\ \cdot&\cdot&1\\ \cdot&\cdot&\cdot&1\\ \end{bmatrix} \quad \mathbf U=\begin{bmatrix} 2&1&1&2\\ &\cdot&\cdot&\cdot\\ &&\cdot&\cdot\\ &&&\cdot\\ \end{bmatrix} \quad \mathbf A=\begin{bmatrix} 2&1&1&2\\ 2&2&2&3\\ 4&2&4&3\\ 0&0&6&-1 \end{bmatrix}$$

之后计算$\mathbf L$的第一列，它由$\mathbf A$的第一列除以$u_{11}=a_{11}$得到：

$$\mathbf L=\begin{bmatrix} 1\\ 1&1\\ 2&\cdot&1\\ 0&\cdot&\cdot&1\\ \end{bmatrix} \quad \mathbf U=\begin{bmatrix} 2&1&1&2\\ &\cdot&\cdot&\cdot\\ &&\cdot&\cdot\\ &&&\cdot\\ \end{bmatrix} \quad \mathbf A=\begin{bmatrix} 2&1&1&2\\ 2&2&2&3\\ 4&2&4&3\\ 0&0&6&-1 \end{bmatrix}$$

此时，我们希望迭代地用上述方法计算$\mathbf A$右下角块的LU分解，但实际上$\mathbf A$的右下角块中包含已经被解出的第一行和第一列的成分，因此需要把对应的成分减去，即：在$\mathbf A$的右下角块中减去$\mathbf L$第一列和$\mathbf U$第一行对应位置的乘积

$$\mathbf L=\begin{bmatrix} 1\\ 1&1\\ 2&\cdot&1\\ 0&\cdot&\cdot&1\\ \end{bmatrix} \quad \mathbf U=\begin{bmatrix} 2&1&1&2\\ &\cdot&\cdot&\cdot\\ &&\cdot&\cdot\\ &&&\cdot\\ \end{bmatrix} \quad \mathbf A^{(1)}=\begin{bmatrix} 2&1&1&2\\ 2&1&1&1\\ 4&0&2&-1\\ 0&0&6&-1 \end{bmatrix}$$

此后，可以直接进行递归求解，求出$\mathbf U$第二行：

$$\mathbf L=\begin{bmatrix} 1\\ 1&1\\ 2&\cdot&1\\ 0&\cdot&\cdot&1\\ \end{bmatrix} \quad \mathbf U=\begin{bmatrix} 2&1&1&2\\ &1&1&1\\ &&\cdot&\cdot\\ &&&\cdot\\ \end{bmatrix} \quad \mathbf A^{(1)}=\begin{bmatrix} 2&1&1&2\\ 0&1&1&1\\ 0&0&2&-1\\ 0&0&6&-1 \end{bmatrix}$$

求出$\mathbf L$第二列：

$$\mathbf L=\begin{bmatrix} 1\\ 1&1\\ 2&0&1\\ 0&0&\cdot&1\\ \end{bmatrix} \quad \mathbf U=\begin{bmatrix} 2&1&1&2\\ &1&1&1\\ &&\cdot&\cdot\\ &&&\cdot\\ \end{bmatrix} \quad \mathbf A^{(1)}=\begin{bmatrix} 2&1&1&2\\ 0&1&1&1\\ 0&0&2&-1\\ 0&0&6&-1 \end{bmatrix}$$

减掉$\mathbf A$中第二行和第二列带来的影响（实际上这里没有影响，因此不变）；然后求出$\mathbf U$第三行：

$$\mathbf L=\begin{bmatrix} 1\\ 1&1\\ 2&0&1\\ 0&0&\cdot&1\\ \end{bmatrix} \quad \mathbf U=\begin{bmatrix} 2&1&1&2\\ &1&1&1\\ &&2&-1\\ &&&\cdot\\ \end{bmatrix} \quad \mathbf A^{(2)}=\begin{bmatrix} 2&1&1&2\\ 0&1&1&1\\ 0&0&2&-1\\ 0&0&6&-1 \end{bmatrix}$$

求出$\mathbf L$第三列：

$$\mathbf L=\begin{bmatrix} 1\\ 1&1\\ 2&0&1\\ 0&0&3&1\\ \end{bmatrix} \quad \mathbf U=\begin{bmatrix} 2&1&1&2\\ &1&1&1\\ &&2&-1\\ &&&\cdot\\ \end{bmatrix} \quad \mathbf A^{(2)}=\begin{bmatrix} 2&1&1&2\\ 0&1&1&1\\ 0&0&2&-1\\ 0&0&6&-1 \end{bmatrix}$$

更新$\mathbf A$：

$$\mathbf L=\begin{bmatrix} 1\\ 1&1\\ 2&0&1\\ 0&0&3&1\\ \end{bmatrix} \quad \mathbf U=\begin{bmatrix} 2&1&1&2\\ &1&1&1\\ &&2&-1\\ &&&\cdot\\ \end{bmatrix} \quad \mathbf A^{(3)}=\begin{bmatrix} 2&1&1&2\\ 0&1&1&1\\ 0&0&2&-1\\ 0&0&0&2 \end{bmatrix}$$

于是

$$\mathbf L=\begin{bmatrix} 1\\ 1&1\\ 2&0&1\\ 0&0&3&1\\ \end{bmatrix} \quad \mathbf U=\begin{bmatrix} 2&1&1&2\\ &1&1&1\\ &&2&-1\\ &&&2\\ \end{bmatrix}$$

注意到循环更新$\mathbf A$的过程其实就是将$\mathbf A$化为$\mathbf U$的形式，与高斯消去法在理论上完全一致，但可以直接采取同时消去$\mathbf A$右下角整个块的方式，人工计算更加迅速。计算机计算时两种方法对稠密矩阵表现相近，稀疏矩阵会有差别

LU分解的用途

尽管LU分解的复杂度仍为$O(n^3)$，但如果已经完成LU分解，再进行解方程（即解$\mathbf {Ly}=\mathbf b$和$\mathbf {Ux}=\mathbf y$两个方程）的复杂度只有$O(n^2)$，因此，对同一个矩阵多个右端项的问题，只需进行一次LU分解，就可以对$p$个右端项达到$O(n^3+pn^2)$的复杂度

选主元技术

在LU分解过程中可能会出现遇到零主元使LU分解无法继续的情况，然而这种情况其实可以通过调换原方程组的方程顺序解决，下面讨论这种技术

---

定理（零主元出现的条件） 对$\mathbf A\in \mathbb R^{n\times n}$，高斯消去过程中不出现零主元当且仅当$\mathbf A$的前$n-1$个顺序主子式均不为零

---

事实上，在进行方程组求解的时候，交换两个方程的顺序并不会影响方程组求解，尽管这会改变方程组的矩阵。因此，对于以求解方程组为目标的LU分解，我们可以适当交换矩阵两行的顺序

注意到在LU分解进行了$k-1$步时，第$k,\cdots,n$行在消去状态上是完全等价的，因此如果此时发现$a_{kk}^{(k)}=0$，但$a_{jk}^{(k)}\neq 0, j>k$，则可以交换第$j$行与第$k$行而不影响方程求解

这被称为选主元技术

事实上，主元较大时，数值误差的传播会减小，因此我们可以强制选择绝对值较大的主元，来增强算法稳定性

但为了保证方程求解的顺利进行，我们需要记录下来发生交换的行。对于单一右端项的情况，我们可以同时交换右端项的元素顺序而不做这样的记录，但这对多右端项的问题不适用，因此较好的办法是用一个置换矩阵$\mathbf P$来记录，形成PLU分解

设PLU分解过程中第$k$步的交换矩阵为$\mathbf P_k$，消去矩阵为$\mathbf M_k$，则有

$$\mathbf M_{n-1}\mathbf P_{n-1}\cdots \mathbf M_2\mathbf P_2\mathbf M_1\mathbf P_1\mathbf A =\mathbf U$$

于是

$$\mathbf P_{n-1}\cdots \mathbf M_2\mathbf P_2\mathbf M_1\mathbf P_1\mathbf A =\mathbf M_{(n-1)}^{-1}\mathbf U$$

由于交换矩阵满足

$$\mathbf P_k^2=\mathbf I$$

由于$\mathbf M_k$起消去第$k$列的作用，只在第$k$列有非对角元素，而$\mathbf P_k$交换第$k$行和某一个第$k$行下方的行，于是

$$\mathbf{P}_j\mathbf{M}_k\mathbf{P}_j=\bar{\mathbf{M}}_k$$

仍然是消去第$k$列的消去阵

因此可以作如下变形

$$(\mathbf P_{n-1}\mathbf M_{n-2}\mathbf P_{n-1})(\mathbf P_{n-1}\mathbf P_{n-2}\mathbf M_{n-3}\mathbf P_{n-2}\mathbf P_{n-1})\cdots (\mathbf P_{n-1}\cdots\mathbf P_3\mathbf P_2\mathbf M_1\mathbf P_2\mathbf P_3\cdots\mathbf P_{n-1})(\mathbf P_{n-1}\cdots\mathbf P_2\mathbf P_1)\mathbf A =\mathbf M_{(n-1)}^{-1}\mathbf U$$

左端记为

$$\bar{\mathbf M}_{n-2}\bar{\mathbf M}_{n-3}\cdots\bar{\mathbf M}_{1}\mathbf P\mathbf A=\mathbf M_{n-1}^{-1}\mathbf U$$

于是，将左侧所有带横线的矩阵逆到右边，可以得到以下形状：

$$\mathbf {PA}=\mathbf {LU}$$

这就得到了一个PLU分解

---

实际计算时，我们只需要在计算过程中加入行交换，并使用向量$\mathbf p$来记录行置换矩阵$\mathbf P$（$\mathbf P$极其稀疏，只需要记录$n$个整数即可），具体地，我们仍然以上一节的例子进行计算（这次我们省略$\mathbf U$的中间表示）：

$$\mathbf p=[1,2,3,4]\quad \mathbf L=\begin{bmatrix} 1\\ \cdot&1\\ \cdot&\cdot&1\\ \cdot&\cdot&\cdot&1\\ \end{bmatrix} \quad \mathbf A=\begin{bmatrix} 2&1&1&2\\ 2&2&2&3\\ 4&2&4&3\\ 0&0&6&-1 \end{bmatrix}$$

交换第1、3行：

$$\mathbf p=[3,2,1,4]\quad \mathbf L=\begin{bmatrix} 1\\ \cdot&1\\ \cdot&\cdot&1\\ \cdot&\cdot&\cdot&1\\ \end{bmatrix} \quad \mathbf A=\begin{bmatrix} 4&2&4&3\\ 2&2&2&3\\ 2&1&1&2\\ 0&0&6&-1 \end{bmatrix}$$

消去第1行、第1列：

$$\mathbf p=[3,2,1,4]\quad \mathbf L=\begin{bmatrix} 1\\ 1/2&1\\ 1/2&\cdot&1\\ 0&\cdot&\cdot&1\\ \end{bmatrix} \quad \mathbf A^{(1)}=\begin{bmatrix} 4&2&4&3\\ 0&1&0&3/2\\ 0&0&-1&1/2\\ 0&0&6&-1 \end{bmatrix}$$

第2行、第2列已经消除完毕：

$$\mathbf p=[3,2,1,4]\quad \mathbf L=\begin{bmatrix} 1\\ 1/2&1\\ 1/2&0&1\\ 0&0&\cdot&1\\ \end{bmatrix} \quad \mathbf A^{(2)}=\begin{bmatrix} 4&2&4&3\\ 0&1&0&3/2\\ 0&0&-1&1/2\\ 0&0&6&-1 \end{bmatrix}$$

交换第3行、第4行（注意此时要同时在$\mathbf L$中重排）：

$$\mathbf p=[3,2,4,1]\quad \mathbf L=\begin{bmatrix} 1\\ 1/2&1\\ 0&0&1\\ 1/2&0&\cdot&1\\ \end{bmatrix} \quad \mathbf A^{(2)}=\begin{bmatrix} 4&2&4&3\\ 0&1&0&3/2\\ 0&0&6&-1\\ 0&0&-1&1/2 \end{bmatrix}$$

消去第3行、第3列：

$$\mathbf p=[3,2,4,1]\quad \mathbf L=\begin{bmatrix} 1\\ 1/2&1\\ 0&0&1\\ 1/2&0&-1/6&1\\ \end{bmatrix} \quad \mathbf A^{(3)}=\begin{bmatrix} 4&2&4&3\\ 0&1&0&3/2\\ 0&0&6&-1\\ 0&0&0&1/3 \end{bmatrix}$$

得到最终结果：

$$\mathbf P=\begin{bmatrix} &&1&\\ &1&&\\ &&&1\\ 1&&& \end{bmatrix} \quad\mathbf L=\begin{bmatrix} 1\\ 1/2&1\\ 0&0&1\\ 1/2&0&-1/6&1\\ \end{bmatrix} \quad \mathbf U=\begin{bmatrix} 4&2&4&3\\ &1&0&3/2\\ &&6&-1\\ &&&1/3 \end{bmatrix}$$

---

在实际计算中，由于$\mathbf L$的对角元素不需要存储，可以直接将其下三角元素存储在$\mathbf A$的已经消去的下三角元素的位置；同时，为了避免大量数据交换，不会真的对矩阵进行重排，而是通过向量$\mathbf p$维护下标映射，进行“懒重排”

上述算法称为部分选主元技术，因为其主元选取只选择待消去列中绝对值最大的元素，实际上可以有更加稳定的全选主元技术，同时进行行、列交换使得整个待消去的右下角子阵中绝对值最大的元素作为主元。这样的技术稳定性更好，但开销也响应更大，得到的分解形如$\mathbf{PAQ}=\mathbf{LU}$

LU分解的稳定性与其增长因子有关，定义为

$$\rho=\frac{\max\limits_{i,j}|u_{ij}|}{\max\limits_{i,j}|a_{ij}|}$$

对解方程的过程作向后误差分析，即使$(\mathbf A+\Delta\mathbf {A})\hat{\mathbf x}=\mathbf b$的$\Delta\mathbf A$满足

$$\frac{||\Delta\mathbf A||_\infty}{||\mathbf A||_\infty}\leq n\rho\varepsilon_{\rm mach}$$

对部分选主元技术，可以证明$\rho\leq 2^{n-1}$，而一般情况下$\rho<10$；如果不选主元则$\rho$不受控制

注 $\rho=2^{n-1}$的一个例子：

$$\mathbf A=\begin{bmatrix}1&-0.5&0&1\\-1&0.75&-0.5&1 \\-1&0.25&0.75&1\\-1&0.25&0.25&1 \end{bmatrix}$$

此矩阵无需置换即可达到部分选主元，但其分解结果为

$$\mathbf L=\begin{bmatrix}1\\-1&1\\-1&-1&1\\-1&-1&-1&1\end{bmatrix}$$

$$\mathbf U=\begin{bmatrix}1&-0.5&0&1\\&0.25&-0.5&2\\&&0.25&4\\&&&8\end{bmatrix}$$

每次消去时，最后一列的元素都恰好翻一倍。可以证明这也是它能够变大的最大倍数

对称正定矩阵的Cholesky分解

对称正定矩阵一定可以进行LU分解，且其LU分解可以写成更好的对称形式：

$$\mathbf A=\mathbf {LDL}^T$$

其中$\mathbf D$为对角阵

如果令$\mathbf D$的所有对角元素取算术平方根得到的矩阵为$\mathbf D^{1/2}$，则

$$\mathbf A=(\mathbf {LD}^{1/2})(\mathbf {LD}^{1/2})^T$$

令

$$\tilde{\mathbf L}=\mathbf{LD}^{1/2}$$

则

$$\mathbf A=\tilde{\mathbf L}\tilde{\mathbf L}^T$$

定理（Cholesky分解） 对称正定矩阵的$\mathbf A=\mathbf {LL}^T$形式分解存在且唯一

---

下面我们仿照之前的思路进行Cholesky分解，注意由于矩阵对陈，我们只需要记录其一半的元素即可

$$\begin{bmatrix} 4&2&4&2\\ \times&10&8&1\\ \times&\times&9&5\\ \times&\times&\times&19 \end{bmatrix}$$

处理第一行和第一列（先对对角元求平方根，然后对其它元素等比例缩小）

$$\begin{bmatrix} 2&1&2&1\\ &10&8&1\\ &\times&9&5\\ &\times&\times&19 \end{bmatrix}$$

更新右下角（用右下角的原始矩阵减去第一列乘第一行的对应位置）

$$\begin{bmatrix} 2&1&2&1\\ &9&6&0\\ &\times&5&3\\ &\times&\times&18 \end{bmatrix}$$

处理第二行和第二列

$$\begin{bmatrix} 2&1&2&1\\ &3&2&0\\ &&5&3\\ &&\times&18 \end{bmatrix}$$

更新右下角

$$\begin{bmatrix} 2&1&2&1\\ &3&2&0\\ &&1&3\\ &&\times&18 \end{bmatrix}$$

处理第三行和第三列

$$\begin{bmatrix} 2&1&2&1\\ &3&2&0\\ &&1&3\\ &&&18 \end{bmatrix}$$

更新右下角

$$\begin{bmatrix} 2&1&2&1\\ &3&2&0\\ &&1&3\\ &&&9 \end{bmatrix}$$

处理右下角元素（开平方即可）

$$\begin{bmatrix} 2&1&2&1\\ &3&2&0\\ &&1&3\\ &&&3 \end{bmatrix}$$

完成计算，得到的分解为

$$\mathbf L=\begin{bmatrix} 2&&&\\ 1&3&&\\ 2&2&1&\\ 1&0&3&3 \end{bmatrix} \quad \mathbf L^T=\begin{bmatrix} 2&1&2&1\\ &3&2&0\\ &&1&3\\ &&&3 \end{bmatrix}$$

Cholesky分解算法非常稳定，其增长因子$\rho\leq1$，存储量与计算量与一般LU分解在渐进阶上一致，但常数约为一半

带状矩阵

对角阵→三对角阵→带状阵

定义（带状矩阵） 对$\mathbf A\in\mathbb R^{n\times n}$，若$\forall i,j:|i-j|>\beta$都有$a_{ij}=0$且$\exists k$，$a_{k,k-\beta},a_{k,k+\beta}$至少一个非零，则称$\mathbf A$为半带宽为$\beta$的带状矩阵

带状矩阵的LU分解结果仍然是不超过原始带宽的带状矩阵

特别地，对三对角矩阵，可以只考虑三对角的元素，用线性向量存储，同时计算LU分解也只计算三对角元素，复杂度降低为线性

对一般的带状矩阵，在不选主元的情况下，时间复杂度$O(\beta^2n)$，空间复杂度$O(\beta n)$

然而带状矩阵的逆是稠密的，因此计算逆矩阵永远是不好的选择

---

不选主元也有稳定LU分解的矩阵：

• 对称正定矩阵
• 对角占优矩阵

定义（按行对角占优矩阵） $\forall i,\ |a_{ii}|\geq\sum\limits_{j=1,j\neq i}^{n}|a_{ij}|$，且至少一个大于号成立

类似可以定义按列对角占优和按行/列严格对角占优

对于带状矩阵，即使选主元，L和U的带宽也不会超过原始矩阵带宽的两倍

稀疏矩阵

稀疏矩阵要求使用非常规的数据结构存储，以达到在计算和存储时优化掉零元素带来的大量浪费的目的

常见的数据结构有：

• COO格式：包含$N$（非零元个数）个三元组(row,column,value)
• 压缩稀疏行（CSR）格式：包含$N$个按行顺序存储的二元组(column,value)和一个$n$长度数组prow，后者记载了每一行在前者中的起始点；相比COO格式节省了$N-n$个整数
• 压缩稀疏列（CSC）
• 若干个一维数组（例如带状矩阵）
• 分块CSR

非线性方程组求解

这实际上是非线性方程求解的多维拓展，下面介绍之前的几种方法的拓展

不动点迭代法

对$\mathbf f(\mathbf x)=\mathbf 0$，其中$\mathbf x\in \mathbb R^n$，$\mathbf f:\mathbb R^n\to\mathbb R^n$，构造一个$\mathbf g$使$\mathbf g(\mathbf x)=\mathbf x\Leftrightarrow \mathbf f(\mathbf x)=\mathbf 0$，则仍然可以进行不动点迭代

回忆：一元不动点迭代的局部收敛准则为$|g'(x^*)|<1$

一元函数的导数对应多元函数的Jacobi矩阵：

$$\mathbf J_g(\mathbf x):\quad J_{g,ij}(\mathbf x)=\frac{\partial g_i(\mathbf x)}{\partial x_j}$$

则局部收敛准则变为$\mathbf J_g(\mathbf x^*)$的所有特征值都满足$|\lambda|<1$

牛顿法

回忆一元牛顿法的迭代公式：

$$x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}$$

类似得到多元牛顿法：

$$\mathbf x_{k+1}=\mathbf x_k-[\mathbf J_f(\mathbf x_k)]^{-1}\mathbf f(\mathbf x_k)$$