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

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

1.引言

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

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

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

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

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

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

- 数乘

- 倍加

- 互换

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

Def： 称Cn×nCn×n中如下形式的矩阵E(u,v;σ)E(u,v;σ)为初等矩阵:

E(u,v;σ)=I−σuvHE(u,v;σ)=I−σuvH

其中非零向量u,v∈Cn,σ≠0u,v∈Cn,σ≠0是实或者复数，即

E(u,v;σ)=⎛⎝⎜⎜⎜⎜⎜1−σu1v1−σu2v1⋮−σunv1−σu1v21−σu2v2⋮−σunv2⋯⋯⋱⋯−σu1vn−σu2vn⋮1−σunvn⎞⎠⎟⎟⎟⎟⎟E(u,v;σ)=(1−σu1v1−σu1v2⋯−σu1vn−σu2v11−σu2v2⋯−σu2vn⋮⋮⋱⋮−σunv1−σunv2⋯1−σunvn)

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

数乘（E1=E(ei,ei;1−α)E1=E(ei,ei;1−α)）

倍加（E2=E(ei,ej;−μ)E2=E(ei,ej;−μ)）

互换（E3=E(ei−ej,ei−ej;1)E3=E(ei−ej,ei−ej;1)）

下面引出初等变换矩阵的一些重要的数学性质

1.两相同向量u，v组成的初等变换矩阵可交换，其积仍然为一个初等矩阵

E(u,v;σ)E(u,v;τ)=E(u,v;σ+τ−στvHu)E(u,v;σ)E(u,v;τ)=E(u,v;σ+τ−στvHu)

证明：

E(u,v;σ)E(u,v;τ)=(I−σuvH))(I−τuvH)=I−σuvH−τuvH+(6)(7)(6)E(u,v;σ)E(u,v;τ)=(I−σuvH))(I−τuvH)(7)=I−σuvH−τuvH+

2.若1−σvHu≠01−σvHu≠0 ,则初等矩阵E(u,v;\sigma)可逆，其逆矩阵也是初等矩阵

E−1(u,v;σ)=E(u,v;τ),τ=σσvHu−1E−1(u,v;σ)=E(u,v;τ),τ=σσvHu−1

3.设v⊥v⊥表示和vv正交的（n-1）维子空间

a.若u∉v⊥u∉v⊥,则E(u,v;σ)E(u,v;σ)有n个线性无关的特征向量，该组特征向量由u和v⊥v⊥中任取一组基向量组成

b.若u∈v⊥u∈v⊥,则E(u,v;σ)E(u,v;σ)仅有n-1个线性无关的特征向量，该组特征向量由v⊥v⊥中任取一组基向量组成

4.det(E(u,v;σ))=1−σvHudet(E(u,v;σ))=1−σvHu

5.对任意非零向量a,b∈Cna,b∈Cn,必可适当选取u,v,σu,v,σ使得E(u,v;σ)a=bE(u,v;σ)a=b

事实上只需要取u,v,σu,v,σ满足vHa≠0,σu=a−bvHavHa≠0,σu=a−bvHa

由初等变换矩阵引出Guass变换矩阵，我们选取σ=−1,u=lk=(0,⋯,0,lk+1,⋯,lnk)T,v=ek=(0,⋯,0,1,0,⋯,0)T,k=1,2,⋯,n−1σ=−1,u=lk=(0,⋯,0,lk+1,⋯,lnk)T,v=ek=(0,⋯,0,1,0,⋯,0)T,k=1,2,⋯,n−1

得到n-1个Guass变换矩阵

det(Lk)=1+eTklk=1det(Lk)=1+ekTlk=1[/b]

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

3.L1(l1)L2(l2)⋯Ln−1(ln−1)=⎛⎝⎜⎜⎜⎜⎜⎜⎜1l21l31⋮ln11l32⋮ln21⋮ln3⋱⋯1⎞⎠⎟⎟⎟⎟⎟⎟⎟L1(l1)L2(l2)⋯Ln−1(ln−1)=(1l211l31l321⋮⋮⋮⋱ln1ln2ln3⋯1)

注意左乘的顺序

3.Gauss消元法

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

- 消元过程

- 回代过程

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

需要注意的是，在一列列的消元过程中，我们需保证aii≠0(i=1,2,…,n)aii≠0(i=1,2,…,n),所以需要利用行互换来保证此条件

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

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

L−1n−1⋯L−1k⋯L−12L−11A(1)=A(n)x=L−1n−1⋯L−1k⋯L−12L−11b(1)=b(n)Ln−1−1⋯Lk−1⋯L2−1L1−1A(1)=A(n)x=Ln−1−1⋯Lk−1⋯L2−1L1−1b(1)=b(n)

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

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

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

4.三角分解法

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

L−1n−1⋯L−1k⋯L−12L−11A(1)x=A(n)x=UxLn−1−1⋯Lk−1⋯L2−1L1−1A(1)x=A(n)x=Ux

故A=L1⋯Lk⋯Ln−2Ln−1U=(I+l1eT1)⋯(I+ln−2eTn−2)(I+ln−1eTn−1)U=(I+l1eT1+⋯+ln−1eTn−1)U(8)(9)(10)(8)A=L1⋯Lk⋯Ln−2Ln−1U(9)=(I+l1e1T)⋯(I+ln−2en−2T)(I+ln−1en−1T)U(10)=(I+l1e1T+⋯+ln−1en−1T)U

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

A=⎛⎝⎜⎜⎜⎜⎜⎜⎜1l21l31⋮ln11l32⋮ln21⋮ln3⋱⋯1⎞⎠⎟⎟⎟⎟⎟⎟⎟⎛⎝⎜⎜⎜⎜⎜⎜⎜u11u12u22u13u23u33⋯⋯⋯⋱u1nu2nu3n⋮unn⎞⎠⎟⎟⎟⎟⎟⎟⎟A=(1l211l31l321⋮⋮⋮⋱ln1ln2ln3⋯1)(u11u12u13⋯u1nu22u23⋯u2nu33⋯u3n⋱⋮unn)

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

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

A=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜a11⋮ar+1,1⋯⋱⋱a1,s+1ar+1,s+1an,n−r⋱⋱an−s,nann⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟A=(a11⋯a1,s+1⋮⋱⋱ar+1,1ar+1,s+1an−s,n⋱⋱an,n−rann)

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

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

首先做Crout分解

A=LU=⎛⎝⎜⎜⎜⎜⎜⎜⎜p1r20p2⋱⋱rn−1pn−1rn0pn⎞⎠⎟⎟⎟⎟⎟⎟⎟⎛⎝⎜⎜⎜⎜⎜⎜⎜10q11q2⋱⋱10qn−11⎞⎠⎟⎟⎟⎟⎟⎟⎟A=LU=(p10r2p2⋱⋱rn−1pn−10rnpn)(1q101q2⋱⋱1qn−101)

然后分两步解决此类问题

追：解Ly=bLy=b

赶：解Ux=yUx=y

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

首先利用Doolittle分解，得A=LUA=LU,对U进一步提取对角矩阵diag(u11,…,unn)diag(u11,…,unn),从而有

U=DD−1U=D(D−1U)=DU0U=DD−1U=D(D−1U)=DU0

故，A=LDU0A=LDU0，由于A对称正定，AT=AAT=A,所以有

AT=(LDU0)T=UT0DLT=A=LDU0AT=(LDU0)T=U0TDLT=A=LDU0

由于分解的唯一性，可知LT=U0LT=U0，从而有

A=LDLTA=LDLT

我们可以记，D1/2=diag(u11−−−√,…,unn−−−√)D1/2=diag(u11,…,unn),从而

A=LD1/2D1/2LT=LD1/2(LD1/2)T=L1LT1A=LD1/2D1/2LT=LD1/2(LD1/2)T=L1L1T

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

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

易知，A=LDLTA=LDLT,则Ax=LDLTxAx=LDLTx

先解Ly=bLy=b,后解LTx=D−1yLTx=D−1y,其中D的逆只需要将对角元素取倒数即可

5.向量和矩阵的范数

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

关于向量范数

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

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

对于齐次性，有∥αx∥=|α|∥x∥,α∈K‖αx‖=|α|‖x‖,α∈K

对于三角不等式，有∥x+y∥≤∥x∥+∥y∥,∀x,y∈V‖x+y‖≤‖x‖+‖y‖,∀x,y∈V

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

l1−范数l1−范数 向量中的元素的绝对值之和

l2−范数l2−范数 向量中的元素的绝对值的平方加起来然后开方

l∞−范数l∞−范数 向量元素中的最大绝对值（使用Cauchy-Schwarz不等式证明三角不等式）

lp−范数lp−范数 向量中的元素的绝对值的p次方加起来然后开p次方根（利用赫尔德不等式即可证明三角不等式）

在最优化理论中可能会涉及加权范数，A为对称正定矩阵，(xTAx)1/2(xTAx)1/2是一种向量范数，记为∥x∥A‖x‖A

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

1范数

∥f(x)∥1=∫ba|f(x)|dx‖f(x)‖1=∫ab|f(x)|dx

p范数

∥f(x)∥p=(∫ba|f(x)|pdx)1p‖f(x)‖p=(∫ab|f(x)|pdx)1p

∞范数

∥f(x)∥∞=max|f(x)|,a≤x≤b‖f(x)‖∞=max|f(x)|,a≤x≤b

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

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

c1∥x∥α≤∥x∥β≤c2∥x∥αc1‖x‖α≤‖x‖β≤c2‖x‖α

称两个范数等价

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

关于矩阵范数

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

∥AB∥≤∥A∥∥B∥‖AB‖≤‖A‖‖B‖

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

∥A∥F=(∑j=1m∑i=0n|aij|2)12‖A‖F=(∑j=1m∑i=0n|aij|2)12

对于Cm×nCm×n上面的任意一种向量诱导范数，都有∥I∥=max∥x∥=1{∥Ix∥=1}‖I‖=max‖x‖=1{‖Ix‖=1}

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

定义：设在两个向量空间Cm,CnCm,Cn中存在向量范数∥∙∥V‖∙‖V, 定义在Cm×nCm×n空间上的矩阵A的由向量范数 ∥∙∥V‖∙‖V诱导所给出的矩阵范数为（其中x不为零向量）∥A∥V=max∥Ax∥V∥x∥V‖A‖V=max‖Ax‖V‖x‖V

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

∥A∥V=max∥Ax∥V∥x∥V=max∥Ax∥V‖A‖V=max‖Ax‖V‖x‖V=max‖Ax‖V

其中第一个max条件为x不为零向量，第二个max条件为∥x∥V=1‖x‖V=1

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

1范数：对矩阵的每一列中的元素取绝对值之后求和，然后选取其中的最大列作为1范数

2范数：矩阵的最大奇异值，也就是矩阵与矩阵的转置的乘积的最大特征值

无穷范数：对于矩阵的每一行的元素取绝对值之后求和，然后选取其中的最大行作为无穷范数

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

设矩阵A属于n*m的复矩阵空间，对于该空间上的某种矩阵范数∥∙∥V‖∙‖V,有∥A∥V<1‖A‖V<1,则矩阵（I±A）（I±A）非奇异，且有∥(I−A)−1∥V≤∥I∥1−∥A∥‖(I−A)−1‖V≤‖I‖1−‖A‖

给出矩阵谱半径的定义

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

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

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

Cond(A)=∥A∥V×∥A−1∥VCond(A)=‖A‖V×‖A−1‖V

对于矩阵的条件数来说，它显然大于等于1，当矩阵恰好是正交矩阵的时候，矩阵的条件数恰好等于1

当矩阵为对称阵，对应的矩阵范数为2范数的时候，此时的条件数称之为谱条件数，其值等于最大特征值除以最小特征值，然后取绝对值