数值分析读书笔记（2）求解线性代数方程组的直接方法

1.引言

矩阵的数值计算一般可以分为直接法和间接法

本章主要介绍$Ax=b$这类线性方程组求解的直接法，数值求解该方程组的基础思想是Gauss消元法

实质是通过一组满秩的初等行变换，将A保秩变换成一个三角矩阵U，此变换过程称为矩阵A的非奇异上三角化

我们的目的就是寻求一个矩阵P，使得PA=U，其中U是一个三角矩阵，其中$Ax=b$和$Ux=\overline b$同解（$\overline b=Pb$）,有效的生成一个P是我们主要研究的问题

2.初等下三角矩阵—Guass变换矩阵

回顾一下线性代数中的三个初等线性变换

• 数乘
• 倍加
• 互换

我们引入一个一般意义上的初等变换矩阵，它把许多常用的线性变换统一在一个框架里面，在数值线性代数中起着重要的意义

Def： 称$C^{n \times n}$中如下形式的矩阵$E(u,v;\sigma)$为初等矩阵:

$$E(u,v;\sigma)=I-\sigma u v^{H}$$

其中非零向量$u,v\in C^{n},\sigma \neq 0$是实或者复数，即

$$E(u,v;\sigma)=\begin{pmatrix} 1-\sigma u_1v_1 & -\sigma u_1v_2& \cdots &-\sigma u_1v_n\\ -\sigma u_2v_1& 1-\sigma u_2v_2 & \cdots & -\sigma u_2v_n\\ \vdots & \vdots &\ddots &\vdots\\ -\sigma u_nv_1&-\sigma u_nv_2 &\cdots &1-\sigma u_nv_n \end{pmatrix}$$

选取不同的$u，v，\sigma$，可以得到许多常用的线性变换矩阵

• 数乘（$E_1=E(e_i,e_i;1-\alpha)$）
• 倍加（$E_2=E(e_i,e_j;-\mu)$）
• 互换（$E_3=E(e_i-e_j,e_i-e_j;1)$）

下面引出初等变换矩阵的一些重要的数学性质
1.两相同向量u，v组成的初等变换矩阵可交换，其积仍然为一个初等矩阵

$$E(u,v;\sigma)E(u,v;\tau)=E(u,v;\sigma+\tau-\sigma \tau v^{H}u)$$

证明：

$$\begin{align} E(u,v;\sigma)E(u,v;\tau) &=(I-\sigma uv^{H}))(I-\tau uv^{H})\\ &=I-\sigma uv^{H}-\tau uv^{H}+ \end{align}$$

2.若$1-\sigma v^{H} u \neq 0$ ,则初等矩阵E(u,v;\sigma)可逆，其逆矩阵也是初等矩阵

$$E^{-1}(u,v;\sigma)=E(u,v;\tau),\tau=\frac{\sigma}{\sigma v^{H}u-1}$$

3.设$v^{\bot}$表示和$v$正交的（n-1）维子空间
a.若$u\notin v^{\bot}$,则$E(u,v;\sigma)$有n个线性无关的特征向量，该组特征向量由u和$v^{\bot}$中任取一组基向量组成
b.若$u\in v^{\bot}$,则$E(u,v;\sigma)$仅有n-1个线性无关的特征向量，该组特征向量由$v^{\bot}$中任取一组基向量组成

4.$det(E(u,v;\sigma))=1-\sigma v^{H}u$

5.对任意非零向量$a,b\in C^{n}$,必可适当选取$u,v,\sigma$使得$E(u,v;\sigma)a=b$

事实上只需要取$u,v,\sigma $满足$v^{H}a\neq 0,\sigma u=\frac{a-b}{v^{H}a}$

由初等变换矩阵引出Guass变换矩阵，我们选取$\sigma =-1,u=l_{k}=(0,\cdots,0,l_{k+1},\cdots,l_{nk})^{T},v=e_{k}=(0,\cdots,0,1,0,\cdots,0)^{T},k=1,2,\cdots,n-1$

得到n-1个Guass变换矩阵

$$L_{k}(l_{k})\equiv E(l_{k},e_{k};-1)=I+l_{k}e_{k}=\begin{pmatrix} 1\\& \ddots\\&&1\\&&l_{k+1,k}&1\\&&l_{k+2,k}&&1\\&&\vdots&&&\ddots\\&&l_"{nk}&&&&1 \end{pmatrix}$$

下面给出Guass变换矩阵的一些性质

1.$det(L_{k})=1+e_{k}^{T}l_{k}=1$

2.Guass变换矩阵的逆只需要将$\sigma$从-1变成+1

3.$L_{1}(l_{1})L_{2}(l_{2})\cdots L_{n-1}(l_{n-1})=\begin{pmatrix}1\\l_{21}&1\\l_{31}&l_{32}&1\\\vdots&\vdots&\vdots&\ddots\\l_{n1}&l_{n2}&l_{n3}&\cdots&1\end{pmatrix}$

注意左乘的顺序

3.Gauss消元法

先介绍一下顺序Gauss消元法，大概分两步

• 消元过程
• 回代过程

在消元过程中，我们不断去左乘Gauss变换矩阵，不断将原矩阵的下三角部分一列列变成0，从而最终变换成一个上三角矩阵

需要注意的是，在一列列的消元过程中，我们需保证$a_{ii}\neq 0 (i=1,2,\dots,n)$,所以需要利用行互换来保证此条件

当然这一切消元过程的前提是，矩阵A应该是非奇异的

经过n-1次的Gauss消元，我们可以得到一个上三角矩阵

$$L_{n-1}^{-1}\cdots L_{k}^{-1}\cdots L_{2}^{-1}L_{1}^{-1}A^{(1)}=A^{(n)}x=L_{n-1}^{-1}\cdots L_{k}^{-1}\cdots L_{2}^{-1}L_{1}^{-1}b^{(1)}=b^{(n)}$$

在回代过程中，由于我们得到了一个上三角矩阵，那么就可以从最底行开始逐步解出x

Gauss消元法的复杂度是$O(n^{3})$，高阶状态下比起克拉默法则运算量要小得多

Gauss消元法过程中，在对各列进行消元的时候，如果主元比较小的话，运算的结果会产生较大的误差，故引入Gauss列主元消元法，即在每一次利用主元消元的步骤之前，把该列中绝对值最大的数所在的行与主元所在的行进行交换

4.三角分解法

我们利用Gauss变换矩阵对Gauss消元法进行进一步的分析

$$L_{n-1}^{-1}\cdots L_{k}^{-1}\cdots L_{2}^{-1}L_{1}^{-1}A^{(1)}x=A^{(n)}x=Ux$$

故$\begin{align}A&=L_{1}\cdots L_{k}\cdots L_{n-2}L_{n-1}U\\&=(I+l_{1}e_{1}^{T})\cdots(I+l_{n-2}e_{n-2}^{T})(I+l_{n-1}e_{n-1}^{T})U\\&=(I+l_{1}e_{1}^{T}+\cdots+l_{n-1}e_{n-1}^{T})U\end{align}$

由此引出矩阵的LU分解，又称Doolittle分解

$$A=\begin{pmatrix}1\\l_{21}&1\\l_{31}&l_{32}&1\\\vdots&\vdots&\vdots&\ddots\\l_{n1}&l_{n2}&l_{n3}&\cdots&1\end{pmatrix}\begin{pmatrix}u_{11}&u_{12}&u_{13}&\cdots &u_{1n}\\&u_{22}&u_{23}&\cdots&u_{2n}\\&&u_{33}&\cdots &u_{3n}\\&&&\ddots &\vdots \\&&&&u_{nn}\end{pmatrix}$$

这里再介绍一下Crout分解，即A=LU中的L是一个下三角矩阵，U是单位上三角矩阵

注意到某些特殊矩阵的三角分解也是比较特殊的，这里引入一类带状对角形矩阵

$$A=\begin{pmatrix}a_{11}&\cdots &a_{1,s+1}&\\\vdots&\ddots&&\ddots\\a_{r+1,1}&&a_{r+1,s+1}&&a_{n-s,n}\\&\ddots&&\ddots\\&&a_{n,n-r}&&a_{nn}\end{pmatrix}$$

上半带宽为s，下半带宽为r，存在LU分解，其中L是下半带宽为r的单位下三角矩阵，U是上半带宽为s的上三角矩阵

对于r=s=1的这一类更加特殊的矩阵，称为三对角矩阵，对于此类矩阵的三角分解，介绍一种“追赶法”

首先做Crout分解

$$A=LU=\begin{pmatrix}p_{1} &&&&0\\r_{2}&p_{2}\\&\ddots&\ddots\\&&r_{n-1}&p_{n-1}\\0&&&r_{n}&p_{n}\end{pmatrix}\begin{pmatrix}1&q_{1}&&&0\\&1&q_{2}\\&&\ddots&\ddots\\&&&1&q_{n-1}\\0&&&&1\end{pmatrix}$$

然后分两步解决此类问题
追：解$Ly=b$
赶：解$Ux=y$

注意到正定对称矩阵的三角分解也是特殊的，这里引入Cholesky分解

首先利用Doolittle分解，得$A=LU$,对U进一步提取对角矩阵$diag(u_{11},\dots,u_{nn})$,从而有

$$U=DD^{-1}U=D(D^{-1}U)=DU_{0}$$

故，$A=LDU_{0}$，由于A对称正定，$A^{T}=A$,所以有

$$A^{T}=(LDU_{0})^{T}=U_{0}^{T}DL^{T}=A=LDU_{0}$$

由于分解的唯一性，可知$L^{T}=U_{0}$，从而有

$$A=LDL^{T}$$

我们可以记，$D^{1/2}=diag(\sqrt{u_{11}},\dots,\sqrt{u_{nn}})$,从而

$$A=LD^{1/2}D^{1/2}L^{T}=LD^{1/2}(LD^{1/2})^{T}=L_{1}L_{1}^{T}$$

此种分解手段称为Cholesky分解，限定对角元素为正，此类分解唯一

上述的Cholesky分解中涉及了开方的运算，下面介绍一种改进的平方根法

易知，$A=LDL^{T}$,则$Ax=LDL^{T}x$

先解$Ly=b$,后解$L^{T}x=D^{-1}y$,其中D的逆只需要将对角元素取倒数即可

5.向量和矩阵的范数

范数是比长度更为一般的概念，有了范数就可以更好的去测度误差的大小

关于向量范数

$$Def: V是数域R/C上的线性空间，对于V中任意的元素x，\\if 存在一个唯一的实函数N（x）与之对应，记为\begin{Vmatrix}x\end{Vmatrix},而且需满足三个条件\\1.非负正定，2.齐次性，3.三角不等式$$

对于非负正定，当仅当x=0，有N（x）=0，否则N（x）> 0;

对于齐次性，有$\begin{Vmatrix}\alpha x\end{Vmatrix}=\begin{vmatrix}\alpha\end{vmatrix} \begin{Vmatrix}x\end{Vmatrix} ,\alpha \in K$

对于三角不等式，有$\begin{Vmatrix}x+y\end{Vmatrix}\leq\begin{Vmatrix}x\end{Vmatrix}+\begin{Vmatrix}y\end{Vmatrix},\forall x,y \in V$

这里介绍几种常见的向量范数

• $l_{1}-范数$ 向量中的元素的绝对值之和
• $l_{2}-范数$ 向量中的元素的绝对值的平方加起来然后开方
• $l_{\infty}-范数$ 向量元素中的最大绝对值（使用Cauchy-Schwarz不等式证明三角不等式）
• $l_{p}-范数$ 向量中的元素的绝对值的p次方加起来然后开p次方根（利用赫尔德不等式即可证明三角不等式）

在最优化理论中可能会涉及加权范数，A为对称正定矩阵，$(x^{T}Ax)^{1/2}$是一种向量范数，记为$\begin{Vmatrix}x\end{Vmatrix}_{A}$

在无限维线性空间中，比如在[a,b]区间中，对于所有的实连续函数集合C[a,b],对于其中的一个元素f(x)也是有类似定义的范数

• 1范数 $$\begin{Vmatrix}f(x)\end{Vmatrix}_1=\int _{a}^{b}\begin{vmatrix}f(x)\end{vmatrix}dx$$
• p范数 $$\begin{Vmatrix}f(x)\end{Vmatrix}_p=\left(\int _{a}^{b}\begin{vmatrix}f(x)\end{vmatrix}^{p}dx \right)^{\frac{1}{p}}$$
• ∞范数 $$\begin{Vmatrix}f(x)\end{Vmatrix}_\infty=max\begin{vmatrix}f(x)\end{vmatrix} ,a\leq x\leq b$$

下面介绍一下范数的等价性

对于任意两个定义好的范数，存在两个与向量x无关的非零正常数c1，c2，有

$$c_1\begin{Vmatrix}x\end{Vmatrix}_\alpha \leq \begin{Vmatrix}x\end{Vmatrix}_\beta\leq c_2\begin{Vmatrix}x\end{Vmatrix}_\alpha$$

称两个范数等价

不难验证，此处的等价性满足数学定义中的等价性的三个条件，即自反，对称，传递

关于矩阵范数

矩阵范数不仅仅满足非负正定，齐次和三角不等式，而且须满足矩阵相乘的相容性，即

$$\begin{Vmatrix}AB\end{Vmatrix} \leq \begin{Vmatrix}A\end{Vmatrix}\begin{Vmatrix}B\end{Vmatrix}$$

这里给出一类特殊的范数， Frobenius范数

$$\begin{Vmatrix}A\end{Vmatrix}_F=\left(\sum_{j=1}^{m}\sum_{i=0}^{n}\begin{vmatrix}a_{ij}\end{vmatrix}^{2}\right)^{\frac{1}{2}}$$

对于$C^{m \times n}$上面的任意一种向量诱导范数，都有$\begin{Vmatrix} I\end{Vmatrix} = \max \limits_{\begin{Vmatrix} x\end{Vmatrix}=1} \{ \begin{Vmatrix} Ix\end{Vmatrix}=1 \}$

这里给出一种范数的定义，即诱导矩阵范数，诱导矩阵范数和向量范数密切相关

定义：设在两个向量空间$C^{m},C^{n}$中存在向量范数$\begin{Vmatrix} \bullet\end{Vmatrix}_{V}$, 定义在$C^{m \times n}$空间上的矩阵A的由向量范数 $\begin{Vmatrix}\bullet\end{Vmatrix}_{V}$诱导所给出的矩阵范数为（其中x不为零向量）$\begin{Vmatrix} A\end{Vmatrix}_{V}=max\frac{\begin{Vmatrix} Ax\end{Vmatrix}_{V}}{\begin{Vmatrix} x\end{Vmatrix}_{V}}$

我们为了解决这个最大值的问题，继续等价定义来优化这个问题

$\begin{Vmatrix} A\end{Vmatrix}_{V}=max\frac{\begin{Vmatrix} Ax\end{Vmatrix}_{V}}{\begin{Vmatrix} x\end{Vmatrix}_{V}}=max\begin{Vmatrix} Ax\end{Vmatrix}_{V}$
其中第一个max条件为x不为零向量，第二个max条件为$\begin{Vmatrix} x\end{Vmatrix}_{V}=1$

我们利用诱导范数的定义可以从原来的向量范数中诱导出三种范数，分别是

1范数：对矩阵的每一列中的元素取绝对值之后求和，然后选取其中的最大列作为1范数
2范数：矩阵的最大奇异值，也就是矩阵与矩阵的转置的乘积的最大特征值
无穷范数：对于矩阵的每一行的元素取绝对值之后求和，然后选取其中的最大行作为无穷范数

关于矩阵的应用，这里引入一个Banach引理

设矩阵A属于n*m的复矩阵空间，对于该空间上的某种矩阵范数$\begin{Vmatrix} \bullet \end{Vmatrix}_{V}$,有$\begin{Vmatrix} A\end{Vmatrix}_{V}<1$,则矩阵$（I±A）$非奇异，且有$\begin{Vmatrix} (I-A)^{-1}\end{Vmatrix}_{V}\leq \frac{\begin{Vmatrix} I\end{Vmatrix}}{1-\begin{Vmatrix} A\end{Vmatrix}}$

给出矩阵谱半径的定义

矩阵的谱半径为矩阵的最大特征值，关于矩阵的谱半径，它不超过其任意一种矩阵范数（当矩阵是Hermite矩阵时，矩阵的2范数恰好等于矩阵的谱半径）

继续给出线性方程组中条件数的定义

在某一矩阵空间中，对于某一矩阵范数，矩阵的条件数=矩阵的范数×矩阵的逆的范数，即

$$Cond(A)=\begin{Vmatrix} A\end{Vmatrix}_{V} \times \begin{Vmatrix} A^{-1}\end{Vmatrix}_{V}$$

对于矩阵的条件数来说，它显然大于等于1，当矩阵恰好是正交矩阵的时候，矩阵的条件数恰好等于1
当矩阵为对称阵，对应的矩阵范数为2范数的时候，此时的条件数称之为谱条件数，其值等于最大特征值除以最小特征值，然后取绝对值