第6矩阵分析75

1、 第六章第六章 matlab使用详解使用详解2012学年选修课学年选修课第6章 矩阵分析 matlab内置了大量的数值计算函数，这些函数封装了常用的数值计算功能。利用这些数值计算函数，能够从烦琐的编程工作中解放出来，集中精力解决问题。本课程将matlab数值计算分为矩阵分析、函数分析和数据分析等初等数值计算内容进行讲解。教学重点：本讲的矩阵分析主要讨论以下问题： 矩阵基本运算，如加、减、乘、除四则运算； 矩阵特征量，如行列式、秩等； 矩阵分解； 矩阵函数； 稀疏矩阵。6.1 matlab数值计算中的矩阵分析 矩阵分析无论是在数学理论还是实际工程问题中都具有重要的应用，例如，线性方程组的解与矩阵

2、除法、矩阵的特征量（如行列式、逆、条件数、秩等）、矩阵分解相关；matlab的最初雏形是为了解决大规模矩阵运算而编写的一系列函数模块。矩阵作为matlab的基本数据结构，一直是matlab的核心，是matlab基本的运算单元，其大部分的内建函数也都支持矩阵作为输入变量，用户在编写自用程序时也应当尽量使用矩阵作为输入变量。6.2 矩阵基本运算 矩阵的加、减、乘、除四则运算、幂运算、比较运算和逻辑运算等代数运算是matlab数值计算最基础的部分。 这里可以粗略地将矩阵运算分为两类，即普通数值运算（四则运算、幂运算）和关系运算（比较运算、逻辑运算。 本节最后将特别介绍矩阵的按位运算。6.2.1 矩阵

3、的加、减 矩阵的加、减运算定义为相应元素的加减。对矩阵a、b，其和（差） ，c也为 矩阵，且cmn=amn+bmn。 矩阵的加、减运算要求参与运算的矩阵具有相同的大小，或者其中之一为标量，例如 矩阵a与标量 的和（差） ， 为 矩阵，且 。cabmnmnxcaxcmnmnmncax应用举例： 例：已知矩阵 a=1:4;5:8， b=0:3;1:4， c=1:3;2:4， x3， 求a+b，a+c，a+x。在matlab中执行结果如下： a+bans = 1 3 5 7 6 8 10 12 a+c ? error using = plusmatrix dimensions must agree.

4、 a+3ans = 4 5 6 7 8 9 10 11注意：矩阵加减运算要求参与运算的矩阵具有相同的大小，即行数、列数相等，否则系统会拒绝运算，并给出出错信息。一种例外的情况是参与运算的矩阵之一为标量。6.2.2 矩阵的线性代数乘法*应用举例： 已知矩阵a=1,2,1;-1,3,1，矩阵b=-1,2;2,3;4,6，求a*b和b*a。 在matlab中执行结果如下： a*b ans = 7 14 11 13 b*a ans = -3 4 1 -1 13 5 -2 26 10 注意：矩阵乘法要求 被乘矩阵的列数与乘矩阵的行数相等，否则系统给出错误信息。可以看到a*b并不等于b*a。6.2.3 矩

5、阵的按位乘法 matlab定义了另一类矩阵乘法，称之为矩阵的按位乘法。由ca*b表示，其中a、b为m x n矩阵，按位乘法所得结果c也为m x n矩阵。说明： 参与按位乘法运算的矩阵应具有相同的大小，或者至少有一个为标量。 矩阵和标量的线性代数乘法与按位乘法是等效的。 矩阵按位乘法是可交换的，即a*bb*a。应用举例：已知a=1,2;3,4; b=-1,1;1,1; c=1,2,3;求a.*b，b.*a，a.*c。matlab执行结果如下： a.*bans = -1 2 3 4 b.*aans = -1 2 3 4 a.*c? error using = timesmatrix dimensi

6、ons must agree. 参与按位乘法运算的矩阵必须具有相同的行数和列数6.2.4 矩阵线性代数除法 矩阵除法是矩阵乘法的逆运算，matlab也定义了两类矩阵除法。第一类是矩阵的线性代数除法，对应于矩阵线性代数乘法的逆运算。矩阵线性代数除法又有两种算子，即右除算子和左除算子，如表所示。矩阵线性代数除法运 算 符名 称说 明/右除b/a对应线性方程x*ab的解左除ab对应线性方程a*xb的解 matlab在进行矩阵线性代数除法运算时，对右除a/b,要求a的列数与b的列数相等；对左除ab,要求a的行数与b的行数相等。6.2.4 矩阵线性代数除法*应用举例：例：已知 a=rand(3), b=

7、1,2,3t,分别求线性方程组a*xb和x*abt的解;在matlab中执行结果如下： a=rand(3)a = 0.9501 0.4860 0.4565 0.2311 0.8913 0.0185 0.6068 0.7621 0.8214 b=1,2,3; x=abx = -1.3480 2.5460 2.2860验证是否a*x等于b a*xans = 1.0000 2.0000 3.0000显然a*xb；ab对应线性方程a*xb的解 b/a对应线性方程 x*ab的解应用举例： x*abt的解为x=b/a matlab命令窗口中输入 x=b/a执行结果如下：x = -1.7100 -0.774

8、0 4.6200 验证是否x*abt 在matlab中输入 x*a执行结果如下：ans = 1 2 3显然x*abt 注意：左除式要求参与运算的矩阵具有相同的行数，右除式要求参与运算的矩阵具有相同的列数。 与矩阵按位乘法相似，matlab也定义了矩阵按位除法。矩阵按位除法也有两种算子，即左除算子和右除算子。 矩阵的按位除法要求参与运算的矩阵具有相同的大小，或至少有一个为标量。6.2.5 矩阵按位除法矩阵按位除法运 算 符名 称说 明./点右除c=a./b,则cmnamn/bmn.点左除c=a.b,则cmnbmn/amn应用举例： 例：已知a=1,3;2,1，b=2,6;2,1，求a./b和a.

9、b在matlab中执行结果如下： a./bans = 0.5000 0.5000 1.0000 1.0000 a.bans = 2 2 1 1c=a./b,则cmnamn/bmnc=a.b,则cmnbmn/amn6.2.6 矩阵的幂 矩阵的幂与矩阵乘法具有紧密的联系，matlab也定义了两类矩阵幂运算。第一类与矩阵线性代数乘法相对应，由 表示，其中 为n阶方阵。 注意：这里特别要注意的是参与幂运算的矩阵必须是方阵。 因为：矩阵乘法要求 被乘矩阵的列数与乘矩阵的行数相等，由于被乘矩阵与乘矩阵为同一矩阵因此要求矩阵a为方阵。cannca aa 应用举例： 已知a=1,2;-1,2，b=1,2;-1

10、,2;3,5，求a2和b2 在matlab命令窗口中输入 a=1,2;-1,2; b=1,2;-1,2;3,5;执行结果如下： a2ans = -1 6 -3 2 b2? error using = mpower matrix must be square.6.2.7 矩阵按位幂 与矩阵乘法相对应，matlab定义了另外一类矩阵幂运算，即矩阵的按位幂运算。矩阵的按位幂运算由c=a.r表示，其中a为m x n矩阵，按位幂运算的结果c也为m x n矩阵，且cmnarmn。 矩阵按位幂运算不要求参与运算的矩阵是方阵。应用举例： 已知a=1,-1;2,2;1,3，求a.2和a.3。 在matlab命令

11、窗口中输入： a=1,-1;2,2;1,3； a.2 a.3执行结果如下： a.2 a.3ans = ans= 1 1 1 -1 4 4 8 8 1 9 1 27 显然a.2不等于a2矩阵按位幂运算不要求参与运算的矩阵是方阵。6.2.8 矩阵按位运算 按位运算是matlab为矩阵设计的一种简洁、高效、安全的运算模式，能够提高代码执行的高效和安全程度。前面介绍的矩阵加减、按位乘除、按位幂都是按位运算符。按位运算符一般有一个（.）作为前导符。.*./ .表：常见的矩阵按位运算符+加-减.*按位乘./按位右除.按位左除.按位幂6.3 矩阵特征量 线性代数中有一些矩阵特征量用于刻画矩阵某方面的性质，包

12、括矩阵的行列式、秩、条件数、以及矩阵的逆等。6.3.1 矩阵的行列式 关于矩阵行列式的概念，这里不作赘述，如有疑问，请参考任何一本线性代数方面的书籍。如n阶矩阵的行列式不等于0，即|a|0时，称矩阵a非奇异，如果限定线性方程组的系数矩阵为方阵，当a非奇异，则线性方程有惟一解。 对n 阶方阵，matlab调用函数det(a)得到矩阵行列式|a|。 注意：只能求方阵的行列式，否则系统报错应用举例： 例：已知amagic(3),求|a|。在matlab命令窗口中输入如下代码： a=magic(3)；a = 8 1 6 3 5 7 4 9 2 det(a)执行结果如下：ans = -360应用举例：

13、例：矩阵a1,2,3;3,4,5求其行列. a=1,2,3;3,4,5a = 1 2 3 3 4 5 det(a)? error using = detmatrix must be square6.3.2 矩阵的逆 若线形方程组ax=b且系数矩阵a非奇异，即|a|不等于0，则线性方程组有惟一解，该惟一解为xa1b，其中a1为a的逆矩阵。对非奇异矩阵a，其逆矩阵a1是满足以下条件的矩阵：aa1= a1a=i （i为单位矩阵）。 matlab调用函数inv(a)求a的逆矩阵a1。 当a奇异时即|a| =0，inv（a）依然能被执行，matlab只会给出警告信息。应用举例： 已知a=1,2,3;1,

14、4,9;1,8,27b=12;16;8求a1 验证aa1= a1a=i ， 求axb的解。在matlab命令窗口中输入如下代码 a=1,2,3;1,4,9;1,8,27； b=12;16;8； cinv(a)执行结果如下：c = 3.0000 -2.5000 0.5000 -1.5000 2.0000 -0.5000 0.3333 -0.5000 0.1667 a*cans = 1.0000 0 0 0.0000 1.0000 0 0.0000 0 1.0000应用举例： c*aans = 1.0000 -0.0000 -0.0000 0 1.0000 0 0.0000 0 1.0000 c*

15、bans = -0.0000 10.0000 -2.6667应用举例：例：b1,2,3;1,2,3;1,2,3,求逆b1b = 1 2 3 1 2 3 1 2 3 inv(b)warning: matrix is singular to working precision.ans = inf inf inf inf inf inf inf inf inf 解为无穷det(b)=06.3.3 矩阵的秩 矩阵线性无关的行数与列数称为矩阵的秩。 matlab中的函数rank用于求矩阵的秩，其调用格式为：r=rank(a)。例如： t=1 1 1;2 2 3t = 1 1 1 2 2 3 r=rank

16、(t)r = 26.4 矩阵分解 矩阵分解是矩阵理论的重要内容，在信号处理、自动控制等众多领域中有着非常广泛的应用。矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，以便于理论分析或数值计算。目前比较常用的矩阵分解方法有：特征分解（evd）、schur（舒尔）分解、cholesky（柯利）分解、三角分解(lu)分解、正交分解（qr）分解、奇异值分解(svd)分解等。6.4 矩阵分解表列出了一些常用的矩阵分解及其对应的matlab实现函数 对n阶方阵a，为标量，v是非零的n维列向量，且满足avv。则称为方阵a的特征值，v是a相对于特征值的特征向量。 在matlab中

17、，计算矩阵a的特征值和特征向量的函数是eig(a)，常用的调用格式有2种：(1) e=eig(a) 矩阵a的全部特征值，构成向量e。(2) v,d=eig(a) 矩阵a的全部特征值，构成对角阵d，并且a的特征向量构成v的列向量。6.4.1 特征分解应用举例： 例：求3阶范德蒙德矩阵a的特征分解。 a=vander(1:3)a = 1 1 1 4 2 1 9 3 1 v,d=eig(a)v = -0.2738 -0.3487 0.2014 -0.5006 0.1162 -0.7710 -0.8213 0.9300 0.6042d = 5.8284 0 0 0 -2.0000 0 0 0 0.17

18、165.8284,-2.0000,0.1716为a的特征值求解线性方程组： 利用矩阵的lu(三角分解)、qr（正交分解）和cholesky（柯利）分解法对求解求解大型方程组非常有用。 优点是运算速度快、可以节省磁盘空间、节省内存6.4.3 cholesky分解 如果a为对称正定矩阵，则cholesky分解可将矩阵a分解为上三角矩阵和其转置矩阵的乘积，即a=rt*r，其中r上三角矩阵。 matlab提供函数chol用于正定矩阵的cholesky分解。chol常用的调用格式有以下两种： r = chol(x) r,p = chol(x) 如果x为正定矩阵，则返回上三角矩阵r，此时p=0，表示函数调

19、用成功；如果x非正定，则前一种调用会产生错误信息，后一种调用不会产生错误，而是将p设为正整数，因此，可以通过查询p的状态检查cholesky分解是否成功，也可以依此判断x的正定性。应用举例： 例：已知a=1,1,1,1;1,2,3,4;1,3,6,10;1,4,10,20,求矩阵a的cholesky分解。在matlab命令窗口中输入： a=1,1,1,1;1,2,3,4;1,3,6,10;1,4,10,20; r,p=chol(a)r = 1 1 1 1 0 1 2 3 0 0 1 3 0 0 0 1p = 0r上为三角矩阵, p=0，表示函数调用成功，也说明矩阵a是正定矩阵应用举例*：已知a

20、16 4 8;4 5 -4;8 -4 22,b=28 5 26t 通过cholesky分解法求ax=b的解因为a=rt*r，因此a*x=b变成rt*r*x=b x=r(rtb) a=16 4 8;4 5 -4;8 -4 22 b=28 5 26 r=chol(a) x=r(rb) lu分解法，可以把任意矩阵分解为下三角矩阵的基本变换形式和上三角矩阵的乘积。即a=lu(或palu) l为下三角矩阵，u为上三角矩阵,p为置换矩阵。 在matlab中函数lu用于矩阵的lu分解，其基本的调用格式为: l,u=lu(a); l,u,p=lu(a)。6.4.4 lu(三角)分解矩阵的lu分解最常用于求解线

21、性方程组axb。首先做lu分解palu，线性方程组axb转换为luxpb，求解过程分两步进行：1）首先求解线性方程组lypb，可得yl(pb);2) 接着求原方程组的解uxy，得xuy。6.4.4 lu(三角)分解*应用举例*： 例：利用lu分解求线性方程组1,2,3;4,5,6;7,8,9x=2;3;4的解。axb转换为luxpb，x=ul(pb)在matlab命令窗口中输入如下代码： a=1,2,3;4,5,6;7,8,9; b=2;3;4; l,u,p=lu(a)； y=l(p*b)； x=uy6.4.5 qr（正交）分解 对于任何长方矩阵a，都可以进行qr分解，其中q为正交矩阵，r为上

22、三角矩阵，即a=qr。 在matlab中常用的调用格式为：q,r=qr(a); 因为a=qr 因此a*x=b变成q*r*x=b x=r(qb)应用举例*：已知a=16 4 8;4 5 -4;8 -4 22,b=28 5 26t通过qr分解法求axb的解a*x=b变成q*r*x=b x=r(qb) a=16 4 8;4 5 -4;8 -4 22 b=28 5 26 q,r=qr(a) x=r(qb)6.5 矩阵函数 矩阵函数是矩阵理论的重要概念，在信号处理、系统控制等领域有着重要的应用。如果将矩阵看作一个线性系统，那么矩阵函数可以看作系统的级联、合成。利用矩阵函数的概念可以得到很多工程应用中有用

23、的工具。6.5.1 常用矩阵函数注意：1）矩阵函数仅支持方阵作为输入参数。2）sqrtm(a)与sqrt(a)不同，b=sqrt(a)是对a中的每个元素求平方根，那么就有b.2=a b=sqrtm(a)是求矩阵a的平方根，那么就有b2=a 应用举例：例：已知a=1,1;0,4求sqrt(a)和sqrtm(a)。在matlab命令窗口中输入如下代码： a=1,1;0,4a = 1 1 0 4 sqrt(a)ans = 1 1 0 2 sqrtm(a)ans = 1.0000 0.3333 0 2.00006.6 稀疏矩阵 实际工程中的数据处理任务面临大容量数据的挑战，当涉及大型矩阵的数值计算时，

24、一个重要的问题是存储和执行效率的问题。稀疏矩阵的概念，正是为了解决这一问题而提出的。从数学性质上看，稀疏矩阵与一般的矩阵没有差别，但在数据存储和执行算法上有着很大的不同。本节在讲述稀疏矩阵时，经常与全(满)矩阵作对比，使大家对稀疏矩阵的概念和使用方法有一个更加透彻的理解。 本节将一般的矩阵称为全矩阵（full matrix），以区别于稀疏矩阵。6.6.1 稀疏矩阵与全矩阵 稀疏矩阵是这样一类矩阵，其元素仅有少数不为0，而大量的元素为0。稀疏矩阵的这种性质使得matlab可以对其采用不同于全矩阵的存储方式和执行算法以提高效率。 matlab利用二维数组存储全矩阵，对零元、非零元不作区分，统一采用

25、浮点数；但在存储稀疏矩阵时只存储非零元及其对应的索引值（整型）。显然，这种存储方式能够大大提高稀疏矩阵的存储效率。应用举例： 例：利用函数sprand创建1000 x 1000稀疏矩阵。在matlab命令窗口中输入如下代码： a=sprand(1000,1000,0.1); %sprand创建稀疏矩阵，前两个参数分别为稀疏矩阵的行数和列数，最后一个参数为非0元所占比例。 full_a=full(a);%full函数将稀疏矩阵转化为全矩阵 whos name size bytes class a 1000 x1000 1144964 double array (sparse)full_a 100

26、0 x1000 8000000 double array 可以看出稀疏矩阵所需存储空间仅是1144964个字节，是对应全矩阵的1144964/8000000=0.1431倍，由此看出稀疏矩阵可以提高效率。 稀疏矩阵的执行效率。 对上例中的两个矩阵，a和full_(a)分别执行乘数运算。 tic;2*a;tocelapsed time is 0.016000 seconds. tic;2*full_a;tocelapsed time is 0.047000 seconds.说明稀疏矩阵提高了运行效率。 %tic表示计时开始，toc表示计时结束，时间差即为中间代码执行时间。应用举例：6.6.2 创

27、建稀疏矩阵 除了通过将全矩阵转换为稀疏矩阵外，matlab还提供了一系列函数用于创建稀疏矩阵，如表所示。稀疏矩阵创建函数函 数 名说 明sparse生成一般的稀疏矩阵speye生成单位稀疏矩阵sprand生成均匀分布随机稀疏矩阵sprandn生成正态分布随机稀疏矩阵sprandsym生成对称随机稀疏矩阵spdiags生成对角稀疏矩阵 具体的： 函数a=sparse(s)将矩阵s转化为稀疏存储方式的矩阵a。 sparse(m,n)：生成一个mn的所有元素都是0的稀疏矩阵。 函数speye(m,n)生成一个mn的稀疏存储单位矩阵 函数aspdiags(b,d,m,n)创建mn的稀疏矩阵，a的第d个

28、对角线为b的列。6.6.2 创建稀疏矩阵此外，还有一些和稀疏矩阵操作有关的函数。例如u,v,s=find(a)：返回矩阵a中非0元素的下标和元素。full(a)：返回和稀疏存储矩阵a对应的完全存储方式矩阵。6.6.2 创建稀疏矩阵应用举例： 例：创建全矩阵a，并将a转化为稀疏矩阵。在matlab命令窗口中输入如下代码： a=eye(4)； a(1,3)=2a = 1 0 2 0 0 1 0 0 0 0 1 0 0 0 0 1 s=sparse(a)s = (1,1) 1 (2,2) 1 (1,3) 2 (3,3) 1 (4,4) 16.6.3 稀疏矩阵操作 一般地，能用于全矩阵的操作函数对稀疏

29、矩阵同样有效，并且具有相似的操作规则，现总结如下： 用于矩阵拼接的函数，如cat，horzcat、vertcat、repmat，若输入参数中有一个为稀疏矩阵，则返回结果为稀疏矩阵； 矩阵变形函数，如ctranspose、flipdim、fliplr、flipud、reshape、rot90、transpose，这些函数都是单输入函数，若输入为稀疏矩阵，则返回结果也为稀疏矩阵； 矩阵结构信息函数，如isempty、isscalar、isvector、length、ndims、numel、size； 矩阵数据类型信息函数，如ischar、isfloat、isinteger、islogical、is

30、number、isreal，这些函数对稀疏矩阵输入返回稀疏矩阵。6.6.3 稀疏矩阵操作应用举例*： 例：1)创建一个单位稀疏矩阵，2)并访问以下下标a(1,1),a(1,3),3)将矩阵a旋转90度。应用举例： 例：1)创建一个单位稀疏矩阵，2)并访问以下下标a(1,1),a(1,3),3)将矩阵a旋转90度。在matlab命令窗口中输入代码： a=speye(3)结果显示为：a = (1,1) 1 (2,2) 1 (3,3) 1 a(1,1)ans = 1 a(1,3)ans = 0 rot90(a)ans = (3,1) 1 (2,2) 1 (1,3) 1应用举例： 例：创建100*10

31、0正态分布随机稀疏矩阵a，求其非零元、非零元个数、并将其非零元替换为1。 a=sprandn(100,100,0.2); whos name size bytes class a 100 x100 21932 double array (sparse)grand total is 1794 elements using 21932 bytes vnz=nonzeros(a）;%返回稀疏矩阵a的非零元组成的列向量 whos name size bytes class a 100 x100 21932 double array (sparse) vnz 1794x1 14352 double ar

32、ray应用举例： nz=nnz(a)%求稀疏矩阵非零元个数nz = 1794 b=spones(a)； %返回与a具有相同稀疏结构的稀疏矩阵，但是非零元均为1。 vnz=nonzeros(b); %返回稀疏矩阵b的非零元组成的列向量 all(vnz=1) %检查替换结果ans = 16.6.4 稀疏矩阵的运算 全矩阵的四则运算对稀疏矩阵都是有效的，但是返回结果有可能是稀疏矩阵或者是全矩阵，这要视具体情况而定，现总结如下： 对于单个稀疏矩阵的输入，大部分函数都返回稀疏矩阵，但也有一部分函数返回全矩阵； 对于多个矩阵输入，若其中有一个及一个以上的全矩阵，则大部分的函数都将返回全矩阵； 对于二元运算，如矩阵的加减、乘、除，只要其中有一个为全矩阵，则结果返回全矩阵；稀疏矩阵与标量的加减运算返回全矩阵； 稀疏矩阵的数乘为稀疏矩阵； 稀疏矩阵的幂仍然是稀疏矩阵。6.6.4 稀疏矩阵的运算应用举例： 例：创建4阶对角稀疏矩阵a，并求其行列式，逆和秩。 d=2,4,3,8; a=spdiags(d,0,4,4)a