线性方程组的直接解法

线性方程组的直接解法

§1线性方程组的条件数

我们曾对一个矩阵计算问题的病态性这个概念作了概略的定性说明。现在我们来考虑如何对线性方程组的病态程度做出定量估计。

设ARnn是非奇异的，bR。我们来考虑线性方程组

Axb， （1）

n

之解x对数据A和b的微小扰动的敏感程度。为此，考虑如下的含参数方程组

(AδA)x()bδb, x(0)x， （2）

其中δAR

nn

，δbR

1

n

，充分小。显然，在0的充分小邻域内

x()(AδA)(bδb)是的可微向量值函数，且易知

x(0)A1(δbδAx)。 （3）



将（3）代入x()的Taylor展开式



x()x(0)x(0)O()，

2

并取范数，可得

x()x

x

A

1

δb

δA

x2

O()， （4）



（5）

并利用不等式

bA

x。

由（4）可得

x()x

x

δAδb2

(A)O()， （6）

Ab

（6）表明，解x的相对误差大约是A与b的相对误差之和的(A)倍。从这个意义来讲，(A)的大小反映了方程组（1）的解对微小扰动的敏感程度。因此我们称(A)为线性

方程组（1）的求解问题的条件数。常用的是对应于p范数，通常记作p(A)；特别当p=2时称作谱条件数。关于条件数的一些最基本的性质，可总结为如下定理。

定理1 设ARnn非奇异。则

（1）(A)1，(A)(A1)，(A)(aA), a0；

（2）2(A)1(A)/n(A)，其中1(A)，n(A)分别表示A的最大和最小奇异值； （3）若A是正规矩阵，则

2(A)max/min；

(A)

(A)

（4）若A是酉矩阵，则2(A)1； （5）若U是酉矩阵，则

2(A)2(UA)2(AU)。

这一性质称作2(A)是酉不变的。

证明留给读者。

定理1的（3）表明，对于一个正规矩阵而言。有“大”的谱条件数的充分必要条件是特征值的模的最大者与最小者之比“大”。但这对一般矩阵来讲是不成立的。即使特征值相等，它们的条件数也可能很大。例如

1A



21





R100100， （7） 21

的特征值都是1，而2(A)2

100

。

在矩阵计算概论中我们曾考察过的方程组

1.0001

0.9999

0.9999



1.0001

x12， x22

其系数矩阵A是对称正定矩阵，它的最大特征值和最小特征值分别为12和20.0002。因此由定理1的（3）知，其谱条件数为2(A)=10000，这表明其系数矩阵

的条件数是很大的，因而才会出现其解对扰动十分敏感的现象。这进一步说明条件数的大小确实在一定程度上反映了线性方程组求解问题的病态程度。

因此，当(A)很大时，我们就说线性方程组（1）的求解问题是病态的，或者说矩阵A是病态的；否则就说其是良态的。

由范数等价定理可知，对于任意两个不同范数定义的条件数(A)和(A)，必存在

两个正数c1和c2，使得

c1(A)(A)c2(A), ARn

nn

。

因此，若一个矩阵A在范数下是病态的，即(A)很大，则(A)/c1以很大，其中c1是与A无关的正数；反过来，若A在范数下有(A)很大，则c2(A)亦很大，c2亦是与A无关的正数。从这个意义上来说，一个矩阵病态与否与具体的范数无关。

一个十分典型的病态矩阵是所谓的Hilber矩阵，其定义为

1121

n11n

12131n1n1

13141n11n2





n

1n1

1()。 ij11



2n212n11

Hn

其条件数2(A)e

3.5n

随着n的增加而非常迅速地增加。因此，其阶数越高，病态程度就越

为严重。

阶数分别为3，5，6，8，10的Hilber矩阵的条件数分别约为5×102，5×105，1.5×107，1.5×1010，1.6×1012。

定理2 设AR

nn

非奇异，bR非零，且δAR

nnn

满足||A

1

||||δA||1。若x

和xδx分别是方程组

Axb 和 (AδA)(xδx)bδb

的解，则

||δx||||x||

1



||δA||||A||

(

||δA||||A||



||δb||||b||

)， （8）

其中(A)||A||||A

证明 注意到

1

||，矩阵范数是由相应的向量范数诱导出的算子范数。

AδA(IδAA)A 和 δAA

11

δAA

1

1，

据第一章的定理3.7和定理3.9知，可逆，而且

(IδAA)

1

1



11δA

A

1

， （9）

因此， xδx(AδA)1(b+δb)(AδA)1b(AδA)1δb。 将xAb代入上式，并移项，可得

11

b(AδA)1δb。 δx(AδA)A

1

上式两边取范数，可得

11

b(AδA)1δx(AδA)A

δb。 （10）

利用恒等式

(AδA)

1

A

1

A(IδAA)δAA

1111

事实上，由A(AδA)δA可得，

(AδA)

1

AI(AδA)δA，

1

1

1

1

即得 (AδA)

1

A

1

A(IδAA)δAA

1

1

1

1

。

1

1

事实上，由 A(AδA)δA(IA)(IA)A

A(IδAA)

(AδA)

1

1

1

1

O得 O，

A(IδAA)δAA

1

1

1

1

1111

A

1

1

A(IδAA)δAA

1

1

1

1

A

1

O，

因而 (AδA)

1

AA(IδAA)δAA

并注意到不等式（9），可得

(AδA)1A1bA1(IδAA1)1δAA1bA1(IδAA1)1δAx 

A

1

(IδAA)

11

δAx

δA1δA

A

11

A

（11） x。

而利用不等式（9），又可得

(AδA)

1

A(IδAA)

111



A1δA

1

1

A

。 （12）

将（11）和（12）代人（10），可得

δx

A1δA

1

1

A

δAxδb。 （13）

（13）两边同除以x，并注意到

bAxA

x，

就有

δxx



A

1

AA

1

1δA

δAδb

。

bA

即不等式（8）。

注

||A

1

定理2的证明过程也证明了：如果ARnn可逆，且δAR

nn

满足

||δ||A||。则1AδA可逆，并且有

||δA||||A||

, ||δA||||A||

||(AδA)

||A

1

A||

1

||





1

1

这表明(A)||A||||A

1

||亦可作为矩阵求逆问题的条件数。

定理2 设ARnn，bR且A非奇异，b非零， x是Axb的精确解，x是近似解，rbAx（称为对应于x的剩余向量）。则有

1

r

xx

x

cond(A)

rb

n

cond(A)b

。

证明 （1）由Axb，得AxAxbAx因而 xxA又有

bA

x,

1

，所以A(xx)r，故xxA1r，

r。

1x



Ab

。

故

xxx

A

1

A

rb

cond(A)

rb

。

（2）由A(xx)r，有

rA

rA(xx)，即(xx)。

1

b，因此，又因xAb，所以xA

1

1x



1A

1

。故

b

xxx



A

1

1

rA

b



1

r

cond(A)b

。

§2基本解法回顾

设ARnn，bR且A非奇异，这一节我们来简要地复习一下求解线性方程组

Axb (2.1)

n

的最基本的解法。

2.1 Gauss消元法

大家知道，当A是稠密矩阵时，目前求解（2.1）的最有效是选主元的Gauss消元法。因此，我们首先复习这一方法的要点。

1．Gauss消元法的基本步骤

（1）利用Gauss变换求A的LU分解：A=LU，其中L是单位下三角矩阵，U是上三角矩阵；

（2）求解方程组Ly=b，得y；

（3）求解方程组Ux=y，得（2.1）的解x。

这一方法简单易行，然而它却有两个致命弱点：

（1）适用范围小。能够对A进行LU分解的矩阵的前提是A从1到n-1阶顺序主子式皆不为零。因此，像

0A

1

1 0

这样一类非常良态（条件数1(A)=1）的矩阵，也不存在LU分解。

（2）数值稳定性差。误差分析的结果表明，按Gauss消去法计算得到的解x，满足

(AE)xb。

U其中

En3A5L

O(

2

)，

分别是L和U计算的结果。是机器精度，由于其中||L||||U||的可能很大，和U这里L

因此数值稳定性差。

理论分析的结果表明，产生上述两个问题的主意原因是零主元素和小主元素的出现。因此，选主元的Gauss消元法就随之而产生。

2．全选主元素的Gauss消元法

（1）求排列方阵P，Q和分解：

PAQ=LU,

其中U是上三角矩阵，L(lij)是满足|lij|1的单位下三角矩阵。

（2）将方程组（2.1）分解为四个简单易解的方程组进行求解。

这样做的结果是弥补了不选主元的Gauss的不足，然而付出的代价也是极其昂贵的。因

n

为选主元必须进行k1次两个元素之间的比较和相应的逻辑判断，这在计算机上是相

k1

2

当费时间的。为了尽可能地减少所进行的比较，人们提出了所谓的部分选主元素的Gauss消元法。

3．部分选主元素的Gauss消元法

在全选主元素的Gauss消元法中，只需取Q=I即可，即在当前的列中选主元素，而不涉及其它元素。

实际计算的经验和理论分析的结果都表明，部分选主元素的Gauss消元法与全选元素的Gauss消元法在数值稳定性方面完全可以媲美，但它的工作量却大为减少（只需进行（n-1）n/2次两个元素之间比较即可），它受到了人们的青睐，成为求解中小型稠密线性方程组最受欢迎的方法之一。

2.2 Cholesky分解法

对于一般方阵，为了消除LU分解的局限性和误差的过分积累，而采用了选主元素的方法。但对于正定矩阵而言，选主元素却是完全不必要的。

设线性方程组（2.1）的系数矩阵A是对称正定的。此时，求解（2.1）的行之有效的方法是所谓的Cholesky分解法，其实质是不选主元素的Gauss消元法。Cholesky分解法的基本步骤如下：

（1）求A的Cholesky分解： AGG，

其中G是对角元素均为正数的下三角矩阵；

（2）求解方程组Gy=b，得y；

（3）求解方程组GTx=y，得（2.1）的解x。

在实际使用时，计算A的Cholesky分解是采用如下方式进行的。 先验地记

g11g21

G

gn1

， gnn

T

g22gn2



然后比较A两边对应的元素，得关系式

j

aij



k1

gikgjk，1jin。

由这一公式，可以逐列求出gij（从第一列开始）。

16



例 对A=4

8

454

8



4进行Cholesky分解和不带平方根的Cholesky分解。 22

1

练习对B=1

2

158

2



8进行Cholesky分解和不带平方根的Cholesky分解。 29

应用Cholesky分解法求对称正定方程组，在数值稳定性方面完全可以与全选主元素的

Gauss消元法媲美，而它的运算量却仅是Gauss消元法的一半，并且省去了元素之间的比较和相应的逻辑判断。因此，Cholesky分解法是目前求解中小型稠密对称正定线性方程组的最佳方法之一。

§5 Toeplitz方程组的解法

设A(aij)R

nn

，如果存在常数1n,2n,,1,0,1,,n1，使得

aij

ji

i,j1,2,,n。

则称A是Toeplitz矩阵；即如果A是Toeplitz矩阵，则它具有如下形状

01

A

1n

10







n1







1

。 10

由此可见，Toeplitz矩阵关于东北-西南对角线是对称的；即关于副对角线对称。具有这样对称性的矩阵通常称作广对角矩阵，即若B(ij)R

nn

是广对角矩阵，则它满足

ijnj1,ni1 i,j1,2,,n。

它等价于B满足

BEBE

T

其中E(en,en1,,e1)是n阶反序单位矩阵。由广对称矩阵的等价定义，易证，非奇异的广对称矩阵的逆矩阵也是广对称矩阵。

在这一节里，我们假定TnR设Tn具有如下形状

nn

是给定的对称正定的Toeplitz矩阵。不失一般性，可

1Tn



1

1

2





n1







2。 （5.1） 11

而且在下面的讨论中，我们将用Tk表示Tn的k阶顺序主子矩阵，即

1

1

Tk

k1

1

1





k1







1



nn

R k1,2,,n。 11

显然，Tk亦是对称正定的Toeplitz矩阵。

下面我们分三部分来讨论系数矩阵为Tn的线性方程组的求解问题和Tn的计算问题。

1

5.1 Yule-Walker方程组

我们先来考虑一类特殊的Toeplitz方程组

Tny(1,2,,n1,n)。 （5.2）

T

其中的1,2,,n1就是（5.1）中确定矩阵的n-1个常数，n是任意给定的实数，这类方程组称作Yule-Walker方程组。

由于这类方程组之右端项的特殊性，所以，若记yk为k阶Yule-Walker方程组

Tkyk(1,,k) k1,2,,n。 （5.3）

T

之解，则可导出由yk确定yk1的关系式。为此，记

zk

yk1

ak

T

kT

, r(,,)。 k1k1

则Tk1yk1(1,,k,k1)可写作

TkT

rkEk

Ekrk

1

zkrk。 akk1

则有

TkzkakEkrkrk， （5.4）

rkEkzkakk1， （5.5）

T

其中Ek表示k阶反序单位矩阵。

由于Tk是对称正定的Toeplitz方程组蕴含着TkEkEkTk，从（5.4）就可得

zkTk(rkakEkrk)ykakEkyk， （5.6）

1

11

将（5.6）代入（5.5），并移项整理，得

(1rkyk)akk1rkEkyk， （5.7）

T

T

注意到

Ik

0

EkykTk

T

1rkEk

T

T

EkrkIk

10EkykTk

10 T

1rkyk

T

和Tk1的正定性，即知1rkyk0。因此，在（5.7）两边同除以1rkyk，即得

ak(k1rkEkyk)/(1rkyk)。 （5.8）

T

T

这样，我们就找到了yk与yk1的之间的关系。从而，我们就可从一阶Yule-Walker方程组的解y1出发，利用公式（5.8）和（5.6）递推地求得方程组（5.2）之解，而且容易算出这一过程所需的运算量为3n/2

此外，令

δk1rkyk, k1,2,,n1， （5.9）

T

2

则有

ykakEkykTT

δk11rk1yk11(rk,k1)

ak

1rkykak(k1rkEkyk)(1ak)δk。 （5.10）

T

T

2

因此，若计算ak时，利用（5.10），便可将计算ak的运算量减少k-2。从而就可将整个求解过程的运算量降为n。

综合上面的讨论，可设计求解（5.2）得算法（参阅课本P100页算法5.1）

2

5.2 一般右端项的Toeplitz方程组

现在我们来考虑一般右端项的Toeplitz方程组

Tnxb。 （5.11）

其中b(1,2,,n)R是已知向量。

类似于Yule-Walker方程组的求解过程，假定xkR是k阶方程组

Tkxk(1,,k) k1,2,,n

T

n

k

之解，则有

ykkEkyk

xk1， （5.12）

k



其中yk是k阶Yule-Walker方程组（5.3）之解，Ek表示k阶反序单位矩阵，k由下式确定

k(k1rkEkxk)/(1rkyk)， （1.13）

T

T

这里， rk(1,,k)。

这样我们便可从x1出发，递推地求得（5.11）之解x，这一求解过程可总结为如下算法： （参阅课本P101页算法5.2）

这一算法的运算量是2n。

2

T

5.3 Toeplitz矩阵的逆

最后，我们来考虑Tn的计算问题。设

Tn1

T

rn1En1

En1rn1

1

1

1

Tn

1

X

T

vvn1



1

n1 1

其中 rn1(1,,n1)，En1表示n-1阶反序单位矩阵，由

Tn1T

rn1En1

En1rn1

1

X

Tv

vIn1

In0

0

， 1

T

可得

Tn1XEn1rn1v

T

In1， （5.14）

Tn1vEn1rn10， （5.15） rn1En1v1， （5.16）

T

由（5.15）可得

vEn1yn1。 （5.17）

其中yn1是n-1阶Yule-Walker方程组的解，将（5.17）代入（5.16），并整理，可得

1/(1rn1yn1)。 （5.18）

T

这样，我们只要求得n-1阶Yule-Walker方程组的解yn1，就可由（5.18）和（5.17）求出Tn的最后一列和最后一行。

下面再来看X(ij)所具有的特性。从（5.14）可得

XTn1Tn1En1rn1v

1

1

11T

Tn1vv/， （5．19）

1T

1其中的最后一个等式用到了Tn1En1rn1En1yn1和（5.17）。由于Tn1(tij)是广对称的，

故从（5.19）可得

ijtijvivj/tni,njvivj/ni,nj(vivjvnivnj)/， （5.20）

这里vi表示v的第i个分量。这也就是说，虽然X并非广对称的，但它的元素ij可由它的关于东北-西南对角线的对称元素ni,nj确定。这样一来，我们就可利用Tn的广对称性和（5.20），从Tn的边缘出发，逐层向内计算，求得Tn的全部元素。为了清楚了解这一过程，我们已n=5为例来说明其具体的计算过程。

t11

t21

T5t31

t41t51

t12t22t32t42t52

t13t23t33t43t53

t14t24t34t44t54

t15t25

t35。 t45t55

1

1

1

„„„„„

最后需指出的是，误差分析的结果表明，上述三种算法与应用Cholesky分解来求解上述三类问题所引起的误差是差不多的。因此，这三种算法是数值稳定性。

TkTkTkTkTkTkTkTkTkTk

(A)(A)(A)(A)(A)(A)(A)(A)(A)(A)(A)(A)上

线性方程组的直接解法

§1线性方程组的条件数

我们曾对一个矩阵计算问题的病态性这个概念作了概略的定性说明。现在我们来考虑如何对线性方程组的病态程度做出定量估计。

设ARnn是非奇异的，bR。我们来考虑线性方程组

Axb， （1）

n

之解x对数据A和b的微小扰动的敏感程度。为此，考虑如下的含参数方程组

(AδA)x()bδb, x(0)x， （2）

其中δAR

nn

，δbR

1

n

，充分小。显然，在0的充分小邻域内

x()(AδA)(bδb)是的可微向量值函数，且易知

x(0)A1(δbδAx)。 （3）



将（3）代入x()的Taylor展开式



x()x(0)x(0)O()，

2

并取范数，可得

x()x

x

A

1

δb

δA

x2

O()， （4）



（5）

并利用不等式

bA

x。

由（4）可得

x()x

x

δAδb2

(A)O()， （6）

Ab

（6）表明，解x的相对误差大约是A与b的相对误差之和的(A)倍。从这个意义来讲，(A)的大小反映了方程组（1）的解对微小扰动的敏感程度。因此我们称(A)为线性

方程组（1）的求解问题的条件数。常用的是对应于p范数，通常记作p(A)；特别当p=2时称作谱条件数。关于条件数的一些最基本的性质，可总结为如下定理。

定理1 设ARnn非奇异。则

（1）(A)1，(A)(A1)，(A)(aA), a0；

（2）2(A)1(A)/n(A)，其中1(A)，n(A)分别表示A的最大和最小奇异值； （3）若A是正规矩阵，则

2(A)max/min；

(A)

(A)

（4）若A是酉矩阵，则2(A)1； （5）若U是酉矩阵，则

2(A)2(UA)2(AU)。

这一性质称作2(A)是酉不变的。

证明留给读者。

定理1的（3）表明，对于一个正规矩阵而言。有“大”的谱条件数的充分必要条件是特征值的模的最大者与最小者之比“大”。但这对一般矩阵来讲是不成立的。即使特征值相等，它们的条件数也可能很大。例如

1A



21





R100100， （7） 21

的特征值都是1，而2(A)2

100

。

在矩阵计算概论中我们曾考察过的方程组

1.0001

0.9999

0.9999



1.0001

x12， x22

其系数矩阵A是对称正定矩阵，它的最大特征值和最小特征值分别为12和20.0002。因此由定理1的（3）知，其谱条件数为2(A)=10000，这表明其系数矩阵

的条件数是很大的，因而才会出现其解对扰动十分敏感的现象。这进一步说明条件数的大小确实在一定程度上反映了线性方程组求解问题的病态程度。

因此，当(A)很大时，我们就说线性方程组（1）的求解问题是病态的，或者说矩阵A是病态的；否则就说其是良态的。

由范数等价定理可知，对于任意两个不同范数定义的条件数(A)和(A)，必存在

两个正数c1和c2，使得

c1(A)(A)c2(A), ARn

nn

。

因此，若一个矩阵A在范数下是病态的，即(A)很大，则(A)/c1以很大，其中c1是与A无关的正数；反过来，若A在范数下有(A)很大，则c2(A)亦很大，c2亦是与A无关的正数。从这个意义上来说，一个矩阵病态与否与具体的范数无关。

一个十分典型的病态矩阵是所谓的Hilber矩阵，其定义为

1121

n11n

12131n1n1

13141n11n2





n

1n1

1()。 ij11



2n212n11

Hn

其条件数2(A)e

3.5n

随着n的增加而非常迅速地增加。因此，其阶数越高，病态程度就越

为严重。

阶数分别为3，5，6，8，10的Hilber矩阵的条件数分别约为5×102，5×105，1.5×107，1.5×1010，1.6×1012。

定理2 设AR

nn

非奇异，bR非零，且δAR

nnn

满足||A

1

||||δA||1。若x

和xδx分别是方程组

Axb 和 (AδA)(xδx)bδb

的解，则

||δx||||x||

1



||δA||||A||

(

||δA||||A||



||δb||||b||

)， （8）

其中(A)||A||||A

证明 注意到

1

||，矩阵范数是由相应的向量范数诱导出的算子范数。

AδA(IδAA)A 和 δAA

11

δAA

1

1，

据第一章的定理3.7和定理3.9知，可逆，而且

(IδAA)

1

1



11δA

A

1

， （9）

因此， xδx(AδA)1(b+δb)(AδA)1b(AδA)1δb。 将xAb代入上式，并移项，可得

11

b(AδA)1δb。 δx(AδA)A

1

上式两边取范数，可得

11

b(AδA)1δx(AδA)A

δb。 （10）

利用恒等式

(AδA)

1

A

1

A(IδAA)δAA

1111

事实上，由A(AδA)δA可得，

(AδA)

1

AI(AδA)δA，

1

1

1

1

即得 (AδA)

1

A

1

A(IδAA)δAA

1

1

1

1

。

1

1

事实上，由 A(AδA)δA(IA)(IA)A

A(IδAA)

(AδA)

1

1

1

1

O得 O，

A(IδAA)δAA

1

1

1

1

1111

A

1

1

A(IδAA)δAA

1

1

1

1

A

1

O，

因而 (AδA)

1

AA(IδAA)δAA

并注意到不等式（9），可得

(AδA)1A1bA1(IδAA1)1δAA1bA1(IδAA1)1δAx 

A

1

(IδAA)

11

δAx

δA1δA

A

11

A

（11） x。

而利用不等式（9），又可得

(AδA)

1

A(IδAA)

111



A1δA

1

1

A

。 （12）

将（11）和（12）代人（10），可得

δx

A1δA

1

1

A

δAxδb。 （13）

（13）两边同除以x，并注意到

bAxA

x，

就有

δxx



A

1

AA

1

1δA

δAδb

。

bA

即不等式（8）。

注

||A

1

定理2的证明过程也证明了：如果ARnn可逆，且δAR

nn

满足

||δ||A||。则1AδA可逆，并且有

||δA||||A||

, ||δA||||A||

||(AδA)

||A

1

A||

1

||





1

1

这表明(A)||A||||A

1

||亦可作为矩阵求逆问题的条件数。

定理2 设ARnn，bR且A非奇异，b非零， x是Axb的精确解，x是近似解，rbAx（称为对应于x的剩余向量）。则有

1

r

xx

x

cond(A)

rb

n

cond(A)b

。

证明 （1）由Axb，得AxAxbAx因而 xxA又有

bA

x,

1

，所以A(xx)r，故xxA1r，

r。

1x



Ab

。

故

xxx

A

1

A

rb

cond(A)

rb

。

（2）由A(xx)r，有

rA

rA(xx)，即(xx)。

1

b，因此，又因xAb，所以xA

1

1x



1A

1

。故

b

xxx



A

1

1

rA

b



1

r

cond(A)b

。

§2基本解法回顾

设ARnn，bR且A非奇异，这一节我们来简要地复习一下求解线性方程组

Axb (2.1)

n

的最基本的解法。

2.1 Gauss消元法

大家知道，当A是稠密矩阵时，目前求解（2.1）的最有效是选主元的Gauss消元法。因此，我们首先复习这一方法的要点。

1．Gauss消元法的基本步骤

（1）利用Gauss变换求A的LU分解：A=LU，其中L是单位下三角矩阵，U是上三角矩阵；

（2）求解方程组Ly=b，得y；

（3）求解方程组Ux=y，得（2.1）的解x。

这一方法简单易行，然而它却有两个致命弱点：

（1）适用范围小。能够对A进行LU分解的矩阵的前提是A从1到n-1阶顺序主子式皆不为零。因此，像

0A

1

1 0

这样一类非常良态（条件数1(A)=1）的矩阵，也不存在LU分解。

（2）数值稳定性差。误差分析的结果表明，按Gauss消去法计算得到的解x，满足

(AE)xb。

U其中

En3A5L

O(

2

)，

分别是L和U计算的结果。是机器精度，由于其中||L||||U||的可能很大，和U这里L

因此数值稳定性差。

理论分析的结果表明，产生上述两个问题的主意原因是零主元素和小主元素的出现。因此，选主元的Gauss消元法就随之而产生。

2．全选主元素的Gauss消元法

（1）求排列方阵P，Q和分解：

PAQ=LU,

其中U是上三角矩阵，L(lij)是满足|lij|1的单位下三角矩阵。

（2）将方程组（2.1）分解为四个简单易解的方程组进行求解。

这样做的结果是弥补了不选主元的Gauss的不足，然而付出的代价也是极其昂贵的。因

n

为选主元必须进行k1次两个元素之间的比较和相应的逻辑判断，这在计算机上是相

k1

2

当费时间的。为了尽可能地减少所进行的比较，人们提出了所谓的部分选主元素的Gauss消元法。

3．部分选主元素的Gauss消元法

在全选主元素的Gauss消元法中，只需取Q=I即可，即在当前的列中选主元素，而不涉及其它元素。

实际计算的经验和理论分析的结果都表明，部分选主元素的Gauss消元法与全选元素的Gauss消元法在数值稳定性方面完全可以媲美，但它的工作量却大为减少（只需进行（n-1）n/2次两个元素之间比较即可），它受到了人们的青睐，成为求解中小型稠密线性方程组最受欢迎的方法之一。

2.2 Cholesky分解法

对于一般方阵，为了消除LU分解的局限性和误差的过分积累，而采用了选主元素的方法。但对于正定矩阵而言，选主元素却是完全不必要的。

设线性方程组（2.1）的系数矩阵A是对称正定的。此时，求解（2.1）的行之有效的方法是所谓的Cholesky分解法，其实质是不选主元素的Gauss消元法。Cholesky分解法的基本步骤如下：

（1）求A的Cholesky分解： AGG，

其中G是对角元素均为正数的下三角矩阵；

（2）求解方程组Gy=b，得y；

（3）求解方程组GTx=y，得（2.1）的解x。

在实际使用时，计算A的Cholesky分解是采用如下方式进行的。 先验地记

g11g21

G

gn1

， gnn

T

g22gn2



然后比较A两边对应的元素，得关系式

j

aij



k1

gikgjk，1jin。

由这一公式，可以逐列求出gij（从第一列开始）。

16



例 对A=4

8

454

8



4进行Cholesky分解和不带平方根的Cholesky分解。 22

1

练习对B=1

2

158

2



8进行Cholesky分解和不带平方根的Cholesky分解。 29

应用Cholesky分解法求对称正定方程组，在数值稳定性方面完全可以与全选主元素的

Gauss消元法媲美，而它的运算量却仅是Gauss消元法的一半，并且省去了元素之间的比较和相应的逻辑判断。因此，Cholesky分解法是目前求解中小型稠密对称正定线性方程组的最佳方法之一。

§5 Toeplitz方程组的解法

设A(aij)R

nn

，如果存在常数1n,2n,,1,0,1,,n1，使得

aij

ji

i,j1,2,,n。

则称A是Toeplitz矩阵；即如果A是Toeplitz矩阵，则它具有如下形状

01

A

1n

10







n1







1

。 10

由此可见，Toeplitz矩阵关于东北-西南对角线是对称的；即关于副对角线对称。具有这样对称性的矩阵通常称作广对角矩阵，即若B(ij)R

nn

是广对角矩阵，则它满足

ijnj1,ni1 i,j1,2,,n。

它等价于B满足

BEBE

T

其中E(en,en1,,e1)是n阶反序单位矩阵。由广对称矩阵的等价定义，易证，非奇异的广对称矩阵的逆矩阵也是广对称矩阵。

在这一节里，我们假定TnR设Tn具有如下形状

nn

是给定的对称正定的Toeplitz矩阵。不失一般性，可

1Tn



1

1

2





n1







2。 （5.1） 11

而且在下面的讨论中，我们将用Tk表示Tn的k阶顺序主子矩阵，即

1

1

Tk

k1

1

1





k1







1



nn

R k1,2,,n。 11

显然，Tk亦是对称正定的Toeplitz矩阵。

下面我们分三部分来讨论系数矩阵为Tn的线性方程组的求解问题和Tn的计算问题。

1

5.1 Yule-Walker方程组

我们先来考虑一类特殊的Toeplitz方程组

Tny(1,2,,n1,n)。 （5.2）

T

其中的1,2,,n1就是（5.1）中确定矩阵的n-1个常数，n是任意给定的实数，这类方程组称作Yule-Walker方程组。

由于这类方程组之右端项的特殊性，所以，若记yk为k阶Yule-Walker方程组

Tkyk(1,,k) k1,2,,n。 （5.3）

T

之解，则可导出由yk确定yk1的关系式。为此，记

zk

yk1

ak

T

kT

, r(,,)。 k1k1

则Tk1yk1(1,,k,k1)可写作

TkT

rkEk

Ekrk

1

zkrk。 akk1

则有

TkzkakEkrkrk， （5.4）

rkEkzkakk1， （5.5）

T

其中Ek表示k阶反序单位矩阵。

由于Tk是对称正定的Toeplitz方程组蕴含着TkEkEkTk，从（5.4）就可得

zkTk(rkakEkrk)ykakEkyk， （5.6）

1

11

将（5.6）代入（5.5），并移项整理，得

(1rkyk)akk1rkEkyk， （5.7）

T

T

注意到

Ik

0

EkykTk

T

1rkEk

T

T

EkrkIk

10EkykTk

10 T

1rkyk

T

和Tk1的正定性，即知1rkyk0。因此，在（5.7）两边同除以1rkyk，即得

ak(k1rkEkyk)/(1rkyk)。 （5.8）

T

T

这样，我们就找到了yk与yk1的之间的关系。从而，我们就可从一阶Yule-Walker方程组的解y1出发，利用公式（5.8）和（5.6）递推地求得方程组（5.2）之解，而且容易算出这一过程所需的运算量为3n/2

此外，令

δk1rkyk, k1,2,,n1， （5.9）

T

2

则有

ykakEkykTT

δk11rk1yk11(rk,k1)

ak

1rkykak(k1rkEkyk)(1ak)δk。 （5.10）

T

T

2

因此，若计算ak时，利用（5.10），便可将计算ak的运算量减少k-2。从而就可将整个求解过程的运算量降为n。

综合上面的讨论，可设计求解（5.2）得算法（参阅课本P100页算法5.1）

2

5.2 一般右端项的Toeplitz方程组

现在我们来考虑一般右端项的Toeplitz方程组

Tnxb。 （5.11）

其中b(1,2,,n)R是已知向量。

类似于Yule-Walker方程组的求解过程，假定xkR是k阶方程组

Tkxk(1,,k) k1,2,,n

T

n

k

之解，则有

ykkEkyk

xk1， （5.12）

k



其中yk是k阶Yule-Walker方程组（5.3）之解，Ek表示k阶反序单位矩阵，k由下式确定

k(k1rkEkxk)/(1rkyk)， （1.13）

T

T

这里， rk(1,,k)。

这样我们便可从x1出发，递推地求得（5.11）之解x，这一求解过程可总结为如下算法： （参阅课本P101页算法5.2）

这一算法的运算量是2n。

2

T

5.3 Toeplitz矩阵的逆

最后，我们来考虑Tn的计算问题。设

Tn1

T

rn1En1

En1rn1

1

1

1

Tn

1

X

T

vvn1



1

n1 1

其中 rn1(1,,n1)，En1表示n-1阶反序单位矩阵，由

Tn1T

rn1En1

En1rn1

1

X

Tv

vIn1

In0

0

， 1

T

可得

Tn1XEn1rn1v

T

In1， （5.14）

Tn1vEn1rn10， （5.15） rn1En1v1， （5.16）

T

由（5.15）可得

vEn1yn1。 （5.17）

其中yn1是n-1阶Yule-Walker方程组的解，将（5.17）代入（5.16），并整理，可得

1/(1rn1yn1)。 （5.18）

T

这样，我们只要求得n-1阶Yule-Walker方程组的解yn1，就可由（5.18）和（5.17）求出Tn的最后一列和最后一行。

下面再来看X(ij)所具有的特性。从（5.14）可得

XTn1Tn1En1rn1v

1

1

11T

Tn1vv/， （5．19）

1T

1其中的最后一个等式用到了Tn1En1rn1En1yn1和（5.17）。由于Tn1(tij)是广对称的，

故从（5.19）可得

ijtijvivj/tni,njvivj/ni,nj(vivjvnivnj)/， （5.20）

这里vi表示v的第i个分量。这也就是说，虽然X并非广对称的，但它的元素ij可由它的关于东北-西南对角线的对称元素ni,nj确定。这样一来，我们就可利用Tn的广对称性和（5.20），从Tn的边缘出发，逐层向内计算，求得Tn的全部元素。为了清楚了解这一过程，我们已n=5为例来说明其具体的计算过程。

t11

t21

T5t31

t41t51

t12t22t32t42t52

t13t23t33t43t53

t14t24t34t44t54

t15t25

t35。 t45t55

1

1

1

„„„„„

最后需指出的是，误差分析的结果表明，上述三种算法与应用Cholesky分解来求解上述三类问题所引起的误差是差不多的。因此，这三种算法是数值稳定性。

TkTkTkTkTkTkTkTkTkTk

(A)(A)(A)(A)(A)(A)(A)(A)(A)(A)(A)(A)上