数值分析：线性方程组的直接解法

本篇文章适合个人复习翻阅，不建议新手入门使用
本专栏：数值分析复习 的前置知识主要有：数学分析、高等代数、泛函分析

线性方程组的直接解法

Gauss（高斯）消元法

设有线性方程组
$$\{\begin{array}{l}{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\\ ⋮\\ a_{n1}{x}_1+a_{n2}{x}_2+\cdots +a_{nn}{x}_n=b_n\end{array}$$

或写为矩阵形式
$$\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& & & \\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right)=\left(\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right)$$

简记为 $Ax=b$

解线性方程组的经典算法是以下所称的“Guass消元法”

定义：Guass（高斯）消元法
从第一列开始，依次将系数矩阵 $A$ 每列对角线以下的元素变为零，得到一个上三角阵；方法是用该列对角线元素所在的行乘以相应倍数消去下面几行；
方程组化为
$$\left(\begin{array}{cccc}{a}_{11}^{(1)}& a_{12}^{(1)}& \cdots & a_{1n}^{(1)}\\ & a_{22}^{(2)}& \cdots & a_{2n}^{(2)}\\ & & \ddots & \\ & & & a_{nn}^{(n)}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right)=\left(\begin{array}{c}{b}_1^{(1)}\\ b_2^{(2)}\\ ⋮\\ b_n^{(n)}\end{array}\right)$$

这一步称为消元

然后依次计算 $x_n,x_{n-1},\dots ,x_1$，得到求解公式
$$\{\begin{array}{l}{x}_n=\frac{b_n^{(n)}}{a_{nn}^{(n)}}\\ x_k=(b_k^{(k)}-\sum _{j=k+1}^n{a}_{kj}^{(k)}{x}_j)/a_{kk}^{(k)},k=n-1,n-2,\dots ,1\end{array}$$

这一步称为回代

Guass消元法的改进：选（列）主元的高斯消元法
一般的Guass消元法的缺点是：无法处理对角线元素为零的矩阵，原因是零无法做分母；并且如果对角线上的元素绝对值较小，可能会产生数量级较大的数，计算结果可能不稳定，从而造成较大误差；因此采用如下称为选主元的高斯消元法，其在原先方法的基础上，在对第 $k$ 列进行消元时，考虑第 $k$ 列的从第 $k$ 行到最后一行的元素，选出绝对值最大的元素，将该元素所在的行与第 $k$ 行交换，再进行消元

高斯消元法与矩阵的LU分解

高斯消元法的矩阵表示

高斯消元法的消元过程相当于对 $A$ 依次左乘 $n-1$ 个下三角阵 $L_1,\dots ,L_{n-1}$，使其成为上三角阵U，记 $L^{-1}=L_{n-1}\cdots L_1$
$$L_1=\left(\begin{array}{ccccc}1& & & & \\ \ast & 1& & & \\ \ast & & 1& & \\ ⋮& & & \ddots & \\ \ast & & & & 1\end{array}\right),L_2=\left(\begin{array}{ccccc}1& & & & \\ & 1& & & \\ & \ast & 1& & \\ & ⋮& & \ddots & \\ & \ast & & & 1\end{array}\right),\dots ,L_{n-1}=\left(\begin{array}{ccccc}1& & & & \\ & 1& & & \\ & & \ddots & & \\ & & & 1& \\ & & & \ast & 1\end{array}\right)$$

命题
设 $L^{-1}=L_{n-1}\cdots L_1$，则有
$$L^{-1}=\left(\begin{array}{ccccc}1& & & & \\ \ast & 1& & & \\ \ast & \ast & 1& & \\ ⋮& ⋮& ⋮& \ddots & \\ \ast & \ast & \ast & \cdots & 1\end{array}\right)$$

其中该矩阵的第一列与 $L_1$ 的第一列相同，第二列与 $L_2$ 的第二列相同……

证明
记 $e_i$ 为单位向量，$l_i$ 的前 $i$ 个元素为0，后面的元素来自矩阵 $L_i$ 的第 $i$ 列的第 $i+1$ 行到最后一行的元素
注意到
$$L_k=I-l_k\cdot e_k^{\prime }$$

$$L_k^{-1}=I+l_k\cdot e_k^{\prime }$$

$$L_k^{-1}\cdot L_{k+1}^{-1}=(I+l_k\cdot e_k^{\prime })(I+l_{k+1}\cdot e_{k+1}^{\prime })=I+l_k\cdot e_k^{\prime }+l_{k+1}\cdot e_{k+1}^{\prime }$$

定义：矩阵的LU分解
矩阵$A\in {\mathbb{C}}^{m\times m}$的LU分解：A=LU，其中L是单位下三角矩阵，U是上三角阵

Gauss消元法相当于先求出矩阵的 LU 分解，再解如下方程组
$$\{\begin{array}{l}Ly=b\\ Ux=y\end{array}$$

高斯消元法
写成元素的形式，得到高斯消元法的具体公式为

第一步：写出矩阵 $U$ 的第一行，$L$ 的第一列
$$u_{1i}=a_{1i}(i=1,2,\dots ,n),l_{i1}=a_{i1}/u_{11},i=2,3,\dots ,n$$

第二步：对每个 $r$，$r=2,\dots ,n$
计算 $U$ 的第 $r$ 行
$$u_{ri}=a_{ri}-\sum _{k=1}^{r-1}{l}_{rk}{u}_{ki},i=r,r+1,\dots ,n$$

计算 $L$ 的第 $r$ 列
$$l_{ir}=(a_{ir}-\sum _{k=1}^{r-1}{l}_{ik}{u}_{kr})/u_{rr}$$

得到完整的矩阵 $L,U$

第三步：求解 $Ly=b$
$$\{\begin{array}{l}{y}_1=b_1\\ y_i=b_i-\sum _{k=1}^{i-1}{l}_{ik}{y}_k,i=2,3,\dots ,n\end{array}$$

第四步：求解 $Ux=y$
$$\{\begin{array}{l}{x}_n=y_n/u_{nn}\\ x_i=(y_i-\sum _{k=i+1}^n{u}_{ik}{x}_k)/u_{ii},i=n-1,n-2,\dots ,1\end{array}$$

其中第二步的公式称为杜利特尔（Doolittle）分解

选主元的高斯消元法
在高斯消元法的基础上，修改为

第零步：初始设置一个 $n$ 维空数组 $I_p$ ，选取矩阵 $A$ 的第一列的主元
$$a_{i_1}=\underset{1\le i\le n}{\max }|a_{i1}|$$

将第1行与第 $i_1$ 行交换，把新的矩阵的 $i-j$ 元重新记为 $a_{ij}$，把 $i_1$ 写入 $I_p$ 的第一个位置

第一步：写出矩阵 $U$ 的第一行，$L$ 的第一列
$$u_{1i}=a_{1i}(i=1,2,\dots ,n),l_{i1}=a_{i1}/u_{11},i=2,3,\dots ,n$$

第二步：对每个 $r$，$r=2,\dots ,n$
计算 $U$ 的第 $r$ 行第 $r$ 列
$$s_i=a_{ir}-\sum _{k=1}^{r-1}{l}_{ik}{u}_{kr},i=r,r+1,\dots ,n$$
选取第 $r$ 列的主元
$$|s_{i_r}|=\underset{r\le i\le n}{\max }|s_i|$$

交换第 $i$ 行与第 $i_r$ 行，把新的矩阵的 $i-j$ 元重新记为 $a_{ij}$，把 $i_r$ 写入 $I_p$ 的第 $r$ 个位置

计算 $U$ 的第 $r$ 行
$$u_{ri}=a_{ri}-\sum _{k=1}^{r-1}{l}_{rk}{u}_{ki},i=r,r+1,\dots ,n$$

计算 $L$ 的第 $r$ 列
$$l_{ir}=(a_{ir}-\sum _{k=1}^{r-1}{l}_{ik}{u}_{kr})/u_{rr}$$
得到完整的矩阵 $L,U$ ，和记录排列阵信息的数组 $I_p$

第三步：对 $i=1,2,3,\dots ,n$ ，
交换向量 $b$ 的第 $i$ 行和第 $i_r$ 行，记新的向量的第 $i$ 个分量仍为 $b_i$ ，计算如下公式
$$y_i=\{\begin{array}{ll}{b}_1,& i=1\\ b_i-\sum _{k=1}^{i-1}{l}_{ik}{y}_k,& i=2,\dots ,n\end{array}$$

第四步：求解 $Ux=y$
$$\{\begin{array}{l}{x}_n=y_n/u_{nn}\\ x_i=(y_i-\sum _{k=i+1}^n{u}_{ik}{x}_k)/u_{ii},\end{array}$$

注：在实际算法中，常将矩阵 $L,U$ 记在同一个矩阵里，矩阵的上三角部分源于 $U$ ，下三角部分（除对角线外）源于 $L$

对称正定矩阵的楚列斯基分解

对称阵的三角分解定理
设 $A$ 是 $n$ 阶对称矩阵，且 $A$ 的所有顺序主子式均不为零，则 $A$ 可唯一分解为
$$A=LD{L}^T$$

其中 $L$ 为单位下三角阵，$D$ 为对角阵

证明
设矩阵 $A$ 的LU分解
$$A=L_1{U}_1$$

其中 $L_1$ 为单位下三角阵，$U_1$ 为对角阵；$U_1$ 可分解为
$$U_1=\left(\begin{array}{cccc}{u}_{11}& & & \\ & u_{22}& & \\ & & \ddots & \\ & & & u_{nn}\end{array}\right)\left(\begin{array}{cccc}1& \frac{u_{12}}{u_{11}}& \cdots & \frac{u_{1n}}{u_{11}}\\ & 1& \cdots & \frac{u_{2n}}{u_{22}}\\ & & \ddots & ⋮\\ & & & 1\end{array}\right)=D{U}_0$$

其中 $D$ 为对角阵，$U_0$ 为单位上三角阵，从而
$$A=L_{1}D{U}_0$$

其转置为
$$A=U_0^{T}D{L}_1^T$$

由 LU 分解的唯一性可得 $L_1=U_0^T$ ，即
$$A=L_{1}D{L}_1^T$$
命题：对称正定矩阵的cholesky（楚列斯基）分解
设 $A$ 为 $n$ 阶对称正定矩阵，则存在一个实的非奇异下三角阵 $L$，使得
$$A=L{L}^T$$

称为 $A$ 的cholesky分解，若限定 $L$ 的对角元为正，则分解唯一

证明思路
注意到以下事实是显然的，在上述 $A=LD{L}^T$ 中，有
$$D=\left(\begin{array}{cccc}{u}_{11}& & & \\ & u_{22}& & \\ & & \ddots & \\ & & & u_{nn}\end{array}\right)=\left(\begin{array}{cccc}\sqrt{u_{11}}& & & \\ & \sqrt{u_{22}}& & \\ & & \ddots & \\ & & & \sqrt{u_{nn}}\end{array}\right)\left(\begin{array}{cccc}\sqrt{u_{11}}& & & \\ & \sqrt{u_{22}}& & \\ & & \ddots & \\ & & & \sqrt{u_{nn}}\end{array}\right)\triangleq D_0^2$$

解对称正定方程组的平方根法

平方根法
解对称正定方程组 $Ax=b$ ，使用cholesky分解
$$A=\left(\begin{array}{cccc}{l}_{11}& & & \\ l_{21}& l_{22}& & \\ ⋮& ⋮& \ddots & \\ l_{n1}& l_{n2}& \cdots & l_{nn}\end{array}\right)\left(\begin{array}{cccc}{l}_{11}& l_{21}& \cdots & l_{n1}\\ & l_{22}& \cdots & l_{n2}\\ & & \ddots & ⋮\\ & & & l_{nn}\end{array}\right)$$
算法如下：

第一步：对于 $j=1,2,\dots ,n$，
$$l_{jj}=(a_{jj}-\sum _{k=1}^{j-1}{l}_{jk}^2{)}^{\frac{1}{2}}$$

$$l_{ij}=(a_{ij}-\sum _{k=1}^{j-1}{l}_{ik}{l}_{jk})/l_{jj},i=j+1,\dots ,n$$

第二步：求解 $Ly=b$
$$y_i=(b_i-\sum _{k=1}^{i-1}{l}_{ik}{y}_k)/l_{ii},i=1,2,\dots ,n$$

第三步：求解 $L^{T}x=y$
$$x_i=(b_i-\sum _{k=i+1}^n{l}_{ki}{x}_k)/l_{ii},i=n,n-1,\dots ,1$$

注：观察第一步中的公式，容易发现 $L$ 是对角占优的，故算法稳定

平方根法的公式需要用到开方运算，为了避免这一点，有如下的改进算法

改进的平方根法
解对称正定方程组 $Ax=b$ ，使用分解式
$$A=\left(\begin{array}{cccc}1& & & \\ l_{21}& 1& & \\ ⋮& ⋮& \ddots & \\ l_{n1}& l_{n2}& \cdots & 1\end{array}\right)\left(\begin{array}{cccc}{d}_1& & & \\ & d_2& & \\ & & \ddots & \\ & & & d_n\end{array}\right)\left(\begin{array}{cccc}1& l_{21}& \cdots & l_{n1}\\ & 1& \cdots & l_{n2}\\ & & \ddots & ⋮\\ & & & 1\end{array}\right)$$

算法如下：

第一步：对于 $i=1,2,\dots ,n$

记
$$t_{ij}=l_{ij}{d}_j$$
从而
$$l_{ij}=(a_{ij}-\sum _{k=1}^{j-1}{t}_{ik}{l}_{jk})d_j,j=1,2,\dots ,i-1$$

$$d_i=a_{ii}-\sum _{k=1}^{i-1}{t}_{ik}^2{l}_{ik}$$

第二步：求解 $Ly=b$
$$y_1=b_1$$

$$y_i=b_i-\sum _{k=1}^{i-1}{l}_{ik}{y}_k,i=2,\dots ,n$$

第三步：求解 $D{L}^{T}x=y$
$$x_n=y_n/d_n$$

$$x_i=y_i/d_i-\sum _{k=i+1}^n{l}_{ki}{x}_k,i=n-1,\dots ,1$$

注：实际算法中常将 $L,D$ 的元素储存在同一个矩阵中

追赶法

设对角占优的三对角方程组
$$\left(\begin{array}{ccccc}{b}_1& c_1& & & \\ a_2& b_2& c_2& & \\ & \ddots & \ddots & \ddots & \\ & & a_{n-1}& b_{n-1}& c_{n-1}\\ & & & a_n& b_n\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{n-1}\\ x_n\end{array}\right)=\left(\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ f_{n-1}\\ f_n\end{array}\right)$$

简记为 $Ax=f$ ，其中 $A$ 是行对角占优的，且 $b_1,c_1,a_n,b_n\ne 0$

追赶法
解对称正定的三对角方程组 $Ax=f$ ，使用分解式
$$A=\left(\begin{array}{ccccc}{b}_1& c_1& & & \\ a_2& b_2& c_2& & \\ & \ddots & \ddots & \ddots & \\ & & a_{n-1}& b_{n-1}& c_{n-1}\\ & & & a_n& b_n\end{array}\right)=\left(\begin{array}{cccc}{\alpha }_1& & & \\ r_2& {\alpha }_2& & \\ & \ddots & \ddots & \\ & & r_n& {\alpha }_n\end{array}\right)\left(\begin{array}{ccccc}1& {\beta }_1& & & \\ & 1& \ddots & & \\ & & \ddots & {\beta }_{n-1}& \\ & & & 1& \end{array}\right)$$

算法如下：

第一步：计算 ${\beta }_i$
$${\beta }_1=c_1/b_1$$

$${\beta }_i=c_i/(b_i-a_i{\beta }_{i-1}),i=2,3,\dots ,n-1$$

第二步：求解 $Ly=f$
$$y_1=f_1/b_1$$

$$y_i=(f_i-a_i{y}_{i-1})/(b_i-a{\beta }_{i-1}),i=2,3,\dots ,n$$

第三步：求解 $Ux=y$
$$x_n=y_n$$

$$x_i=y_i-{\beta }_i{x}_{i+1},i=n-1,n-2,\dots ,1$$

参考书籍：《数值分析》李庆扬 王能超 易大义 编