大型矩阵快速求逆算法的研究.doc
《大型矩阵快速求逆算法的研究.doc》由会员分享，可在线阅读，更多相关《大型矩阵快速求逆算法的研究.doc（29页珍藏版）》请在冰点文库上搜索。
说明书
矩阵是高等代数的重要组成部分，是许多数学分支研究的重要工具。
并且矩阵作为数学工具之一有其重要的使用价值，它常见于很多学科中，如：
线性代数、线性规划、统计分析，以及组合数学等。
在实际生活中，很多问题都可以借用矩阵抽象出来进行表述并进行运算，如在各循环赛中常用的赛况表格等，矩阵的概念和性质相对矩阵的运算较容易理解和掌握。
在矩阵的运算和应用中，仅逆矩阵的求解方法及算法问题就值得我们去好好研究，尤其是对大型矩阵的逆矩阵的研究。
近年来，随着互联网的高速发展，计算机内部的运算量也急剧增加，如何把浩瀚的数据准确地计算出结果，并且加快它的计算速度，己成为一个备受关注的研究课题。
随着计算机应用领域的发展，逆矩阵运算的需求越来越大。
现阶段大型逆矩阵运算都是由软件实现，如Matlab等数学软件。
逆矩阵运算软件的普及，将使计算机的逆矩阵运算性能得到几何级数的提高。
伴随矩阵作、初等变换等算法作为逆矩阵运算的一个重要组成部分，对其研究的意义也就很大。
我们所做的研究方向，仅是利用Matlab数学软件和我们所学到的求逆矩阵的知识编写几种求逆矩阵的算法。
为了使内容更完整充实，文中列举了几种基本的求逆矩阵的方法以及整理了一些特殊矩阵求逆的一些简便方法等。
最后简要说明了一般矩阵的广义逆矩阵。
由于考虑自身学识的不足，我们并未对广义逆矩阵做详细的展开。
希望有兴趣的读者能自己探索。
大
型
矩
阵
快
速
求
逆
算
法
的
研
究
信息09本
耀魂雨
大型矩阵快速求逆算法的研究
本文首先介绍了逆矩阵的定义以及逆矩阵的相关性质。
其次，根据逆矩阵的相关理论主要介绍了矩阵求逆的几种常用方法。
如定义法、伴随矩阵法、初等变换法、三角分解法、分块矩阵法等,并运用软件matlab7.0对一些方法实现了程序化。
且通过多次检验证明了所编程序的正确性。
文章最后简要阐述了对一些特殊矩阵的求逆算法（如对称正定矩阵、有理矩阵），还有针对一般矩阵的广义逆矩阵的作了浅层次的探究。
本文的相关研究不仅对提高矩阵求逆的速度和准确度起着较为重要的作用。
而且对逆矩阵应用的进一步发展有着深远的意义。
关键词：
逆矩阵求逆算法Matlab7.0广义矩阵
逆矩阵的定义：
对于n阶方阵A，如果有一个n阶方阵B，使AB=BA=E，则说方阵A是可逆的，并把方阵B称为A的逆矩阵（其中只有方阵才有逆矩阵的概念）。
矩阵可逆的条件：
定理1若方阵A可逆，则A的行列式不等于0。
定理2若A的行列式不等于0，则A可逆，且
故，求逆矩阵可以用逆矩阵的定义来求或者利用求伴随矩阵法的方法，还可以用行初等变换法（本文仅以行初等变换为例）、矩阵的三角分解、矩阵分块等方法求得逆矩阵。
一、五种基本求逆矩阵的方法
1.1定义法求逆矩阵
用一个例子来解释定义法求逆矩阵的方法：
例：
求矩阵A的逆矩阵，其中
解：
因为，所以存在.
设=，由定义知,
所以.
由矩阵乘法得
.
由矩阵相等可解得
故
1.2伴随矩阵求逆矩阵（附算法）
用求伴随矩阵的方法求逆矩阵的原理分析如下：
设A=（）为n阶矩阵，为|A|中元素的代数余子式,（i,j=1,2,…,n）,则称矩阵
÷
÷
÷
÷
÷
ø
ö
ç
ç
ç
ç
ç
è
æ
=
nn
n
n
n
n
A
A
A
A
A
A
A
A
A
A
L
L
L
L
L
L
L
L
2
1
2
22
12
1
21
11
*
为A的伴随矩阵。
根据上述定理2:
若A的行列式不等于0,则A可逆，且,
1.3初等变换法求逆矩阵（附算法）
用初等变换法的方法求逆矩阵的原理分析如下：
引理1.3.1对mxn阶矩阵A，施行一次初等行变换，相当于在A的左边乘以相应m阶初等矩阵；对A施行一次初等列变换，相当于在A的右边乘以相应的n阶初等矩阵。
公理1.3.1初等矩阵都是可逆矩阵，其逆矩阵还是初等矩阵。
推论1.3.1A可逆的充要条件为A可表为若干初等矩阵之积。
即
推论1.3.2A可逆，则A可由初等行变换化为单位矩阵。
（1）
由矩阵初等变换的这些性质可知，若A可逆，构造分块矩阵（A︱E），其中E为与A同阶的单位矩阵，那么
（2）
由
（1）式代入
（2）式左边，
有
上式说明分块矩阵（A︱E）经过初等行变换，原来A的位置变换为单位阵E，原来E的位置变换为我们所要求的,即
1.4三角分解求逆矩阵（附算法）
直接可以得到
从此方法把A矩阵进行LU分解，并得到A的逆矩阵。
1.5分块矩阵法求逆
此方法可将阶数较高的矩阵化为阶数较低的矩阵再求其逆，使计算简化
1、利用分块矩阵的乘法，求可逆矩阵的逆阵.
设分块矩阵其中A，C均可逆，求
解：
令，于是有
可得
所以，从而
2、利用分块矩阵的初等变换，降阶求逆
定理1.5.1设A和B分别是可逆矩阵，C是k×r矩阵，D是r×k矩阵,“可逆”或“可逆”,则矩阵可逆。
上例说明定理3可以应用到任何阶（奇数阶或偶数阶）可逆矩阵的求逆问题中。
通常它使2n阶可逆阵的求逆问题转化为一些n阶可逆阵的求逆问题，而使2m+1阶可逆阵的求逆问题转化为一些m阶和一些m+1阶可逆阵的求逆问题，从而使计算的难度降低。
因而这种求逆矩阵的方法也称之为“降阶法”。
应当指出，当A,B可逆时，相平行地，对于形如的矩阵（这里A,B分别是k阶和r阶可逆矩阵,C是k×r矩阵，D是r×k矩阵）有着类似的结论：
二、算法要求及算法分析：
这里仅运用Matlab7.0软件对求伴随矩阵、初等行变换及直接三角分解法求逆矩阵这三种方法编写了算法，并进行了适当的算法分析。
2.1算法要求
首先，不管是用求伴随矩阵A*还是利用初等变换求A的逆矩阵，都要先判定A矩阵是否为方阵。
通过判断A矩阵的行数和列数是否相等来判断A矩阵是否为方阵。
其次，在确定A矩阵为方阵后还须判别A是否可逆。
这里我们可以利用定理1，通过判别A的行列式是否为0，来确定A是否存在逆矩阵。
最后，运用求伴随矩阵、初等行变换及直接三角分解法的计算方法来编写代码。
2.2三种求逆算法及解释分析
以下是运用Matlab7.0编写的求逆矩阵的三个函数，以及对这三个函数的分析及比较：
1、利用求伴随矩阵的方法编写的求逆算法（代码如下）
M函数函数名：
[F]=solvenif1（A）
function[F]=solvenif1（A）
n=rank（A）;[x,y]=size（A）;
ifx~=y%判断A矩阵是否为方阵
disp（'请注意：
因为该矩阵非方阵，所以该矩阵无逆.'）
F='A无逆矩阵';
return;
else
%disp（'请注意：
该矩阵为方阵.'）
d=det（A）;
ifd==0%判断A矩阵是否存在逆矩阵
disp（'请注意：
因为该矩阵非满秩，所以该矩阵无逆.'）
F='A无逆矩阵';
return;
else
%disp（'请注意：
因为该矩阵满秩，所以该矩阵逆存在.'）
fori=1:
n
forj=1:
n
B=A;
B（i,:
）=[];B（:
j）=[];
C（j,i）=（（-1）^（i+j））*det（B）/d;
end
end
F=C;
disp（'该矩阵的逆F='）
end
end
分析：
利用求伴随矩阵的方法编写的求逆算法最主要的就是求元素的代数余子式。
故一定要把元素所在的行和列删除，且需要剩余部分组成的余子式输出。
解释：
B（i,:
）=[];B（:
j）=[];%是为了选取元素的代数余子式。
其中
（i，j=1,2,3﹒﹒﹒n）为的代数余子式，B为（n-1）阶方阵。
C（j,i）=（（-1）^（i+j））*det（B）/d;%是求A*与A的行列式d的比，结果即为（=/d）。
运行代码：
⑴
>>clear;A=[27-19;-113-37],[F]=solvenif1（A）
A=
27-19
-113-37
请注意：
因为该矩阵非方阵,所以该矩阵无逆.
F=
A无逆矩阵
⑵
>>clear;A=[23-1;-13-3;115-11],[F]=solvenif1（A）
A=
23-1
-13-3
115-11
请注意：
因为该矩阵非满秩,所以该矩阵无逆.
F=
A无逆矩阵
⑶
>>clear;A=[123;221;343];[F]=solvenif1（A）
该矩阵的逆F=
F=
1.00003.0000-2.0000
-1.5000-3.00002.5000
1.00001.0000-1.0000
>>A*F
ans=
100
010
001
注：
该算法对任何大型可逆方阵都适用
2、利用初等行变换求逆矩阵（代码如下）
M函数函数文件名：
[F]=solvenif2（A）
function[F]=solvenif2（A）
n=rank（A）;[x,y]=size（A）;
ifx~=y
disp（'请注意：
因为该矩阵非方阵，所以该矩阵无逆.'）
F='A无逆矩阵';
else
%disp（'请注意：
该矩阵为方阵.'）
d=det（A）;
ifd==0
disp（'请注意：
因为该矩阵非满秩，所以该矩阵无逆.'）
F='A无逆矩阵';
else
%disp（'请注意：
因为该矩阵满秩，所以该矩阵逆存在.'）
B=eye（x）;
A=[A,B];%把A矩阵和单位阵B按行合并
y=2*y;
fork=1:
x%该for循环是在选主元
max=abs（A（k,k））;
r=k;
forL=k+1:
x
ifmaxmax=abs（A（L,k））;
r=L;
end
end
t=A（k,:
）;
A（k,:
）=A（r,:
）;
A（r,:
）=t;
s=A（k,k）;
forj=1:
y
A（k,j）=A（k,j）/s;%求行乘数
end
fori=1:
x
ifi~=k
forj=k+1:
y
A（i,j）=A（i,j）-A（i,k）*A（k,j）;%行变换
end
end
end
end
F=A（:
x+1:
y）;
disp（'该矩阵的逆F='）
end
end
分析：
利用初等行变换求逆矩阵相对求伴随矩阵的方法就麻烦的多了。
对于一个大型的矩阵在按照固定的循环做初等行变换的过程中，可能会遇到恰好某一行的主对角线上的元素（主元）为0。
则此时就需要换行了。
那“换行的标准是什么，要跟哪行换？
”就要处理好。
遇到这种情况，这里我们可以把该行与下面主元的绝对值最大的元素换行。
依次循环，最终输出后（n+1）列到2n列的矩阵即为所求。
解释：
行乘数：
在每次行变换的过程中，下一行需减去上一行的倍。
这里的即为行乘数。
运行代码：
>>A=unifrnd（0,100,6,6）,[F]=solvenif2（A）
%A=unifrnd（0,100,6,6）生成[0,100]区间上的连续型均匀分布6行6列的随机数矩阵
A=
95.012945.646892.181341.027013.88911.5274
23.11391.850473.820789.365020.276574.6786
60.684382.140717.62665.789119.872244.5096
48.598244.470340.570635.286860.379293.1815
89.129961.543293.547081.316627.218846.5994
76.209779.193791.69040.986119.881441.8649
该矩阵的逆F=
F=
0.05740.02750.03650.0015-0.0622-0.0241
-0.0442-0.0253-0.0147-0.00570.05340.0156
-0.0139-0.0061-0.0211-0.00090.01530.0186
-0.0169-0.0076-0.0061-0.00370.0313-0.0060
-0.0364-0.0461-0.04710.02530.06130.0090
0.02720.03310.0299-0.0019-0.0513-0.0065
>>A*F
ans=
1.0000-0.00000.00000.0000-0.00000.0000
01.00000-0.00000.00000.0000
-0.0000-0.00001.0000-0.00000.00000.0000
-0.0000-0.0000-0.00001.0000-0.00000
-0.0000-0.0000-0.00000.00001.00000.0000
-0.0000000.0000-0.00001.0000
3、利用三角分解的方法（LU分解）编写的求逆算法（代码如下）
M函数函数名：
[F]=solvenif3（A）
function[L,U,F]=solvenif3（A）
%实现对矩阵A的LU分解，U为上三角矩阵，L为下三角矩阵，并输出A的逆矩阵
n=rank（A）;[x,y]=size（A）;
ifx~=y
disp（'请注意：
因为该矩阵非方阵,不能进行LU分解,所以该矩阵无逆.'）
L='L不存在';
U='U不存在';
F='A无逆矩阵';
else
%disp（'请注意：
该矩阵为方阵.'）
d=det（A）;
ifd==0
disp（'请注意：
因为该矩阵非满秩,不能进行LU分解,所以该矩阵无逆.'）
L='L不存在';
U='U不存在';
F='A无逆矩阵';
else
%disp（'请注意：
因为该矩阵满秩，所以该矩阵逆存在.'）
B=eye（n）;
A=[A,B];
L=zeros（n,n）;
U=zeros（n,2*n）;
fori=1:
n
L（i,i）=1;
end
fork=1:
n
forj=k:
2*n
U（k,j）=A（k,j）-sum（L（k,1:
k-1）.*U（1:
k-1,j）'）;
end
fori=k+1:
n
L（i,k）=（A（i,k）-sum（L（i,1:
k-1）.*U（1:
k-1,k）'））/U（k,k）;
end
end
C=U（:
n+1:
2*n）;
U=U（:
1:
n）;
F=U\C;
end
disp（'U为上三角矩阵，L为下三角矩阵，F为A的逆矩阵'）
end
注：
利用三角分解法进行LU分解具有唯一性。
上述的三角分解为直接分解，没有进行选主元的过程，故当A矩阵出现病态时，可能导致结果不正确。
对选主元的三角分解这里不在进行详细说明了。
有兴趣的读者请看附录一的选主元三角分解求逆矩阵的Matlab代码。
运行代码：
>>A=unifrnd（0,100,5,5）;[L,U,F]=solvenif3（A）
U为上三角矩阵，L为下三角矩阵，F为A的逆矩阵
L=
1.00000000
0.72671.0000000
0.88461.34861.000000
0.57300.96331.17341.00000
0.74281.13730.1601-0.63091.0000
U=
58.279222.595020.906956.782941.5375
041.561622.789438.15880.3159
0029.1057-95.770250.2685
000103.3674-81.5894
0000-13.9437
F=
0.0307-0.0596-0.07810.04930.0950
-0.0129-0.0129-0.06600.02840.0867
-0.00110.00650.0532-0.0075-0.0624
-0.00430.03600.0396-0.0260-0.0566
-0.00560.03810.0646-0.0452-0.0717
>>A*F
ans=
1.00000.00000.000000.0000
-0.00001.0000-0.00000.00000.0000
0.00000.00001.00000.00000.0000
-0.00000.00000.00001.00000.0000
0.0000-0.0000-0.00000.00001.0000
2.3三种算法的简单比较
运用Matlab7.0中的计时函数对上述三种求逆的算法进行比较：
随机生成一个20阶的方阵，分别运用这两种方法求逆，并计算运算时间。
>>A=unifrnd（0,100,20,20）
结果略
>>tic;[F]=solvenif1（A）;toc
该矩阵的逆F=
Elapsedtimeis0.047000seconds.
>>tic;[F]=solvenif2（A）;toc
该矩阵的逆F=
Elapsedtimeis0.015000seconds.
>>tic;[L,U,F]=solvenif3（A）;toc
U为上三角矩阵，L为下三角矩阵，F为A的逆矩阵
Elapsedtimeis0.190000seconds.
仅从计算时间上看，分别运用这三种算法对同一个20阶方阵求逆，初等行变换求逆的效率与另外的两种相比较好。
三、特殊矩阵的求逆方法
3.1对称正定矩阵求逆
设C为n阶对称正定矩阵，y、x为n维向量，其关系式为：
y=C·x（3-1）
确定了上的一个映像，如能写出逆关系：
x=B·y（3-2）
则B为C的逆阵，即B=
现将式（3-1）写成
......
......（3-3）
.......
......
因C对称正定，必有0,
升级会员 