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

数值分析读书笔记（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_

下面给出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范数的时候，此时的条件数称之为谱条件数，其值等于最大特征值除以最小特征值，然后取绝对值