《数字信号处理》教学参考资料（MATLAB 5手册）第9章 稀疏矩阵

China-bub.com 下载 第9章稀疏矩阵 在许多问题中提到了含有大量0元素的矩阵，这样的矩阵称为稀疏矩阵。比如求解普通或 者部分微分方程的数值解。为了节省存储空间和计算时间，MATLAB考虑到矩阵的稀疏性， 在对它运算时有特殊的命令。 9.1矩阵为什么稀疏 一个稀疏矩阵中有许多元素等于零，这便于矩阵的计算和保存。如果MATLAB把一个矩 阵当作稀疏矩阵，那么只需在m×3的矩阵中存储m个非零项。第1列是行下标，第2列是列下 标，第3列是非零元素值，不必保存零元素。如果存储一个浮点数要8个字节，存储每个下标 要4个字节，那么整个矩阵在内存中存储需要16×m个字节。 ■例9.1 A=eye(1000): 得到一个1000×1000的单位矩阵，存储它需要8Mb空间。如果使用命令： B=speye(1000); 用一个1000×3的矩阵来代表，每行包含有一个行下标、列下标和元素本身。现在只需 16Kb的空间就可以存储1000×1000的单位矩阵，它只需要满单位矩阵的0.2%存储空间。对于 许多的广义矩阵也可这样来作。 ■ 稀疏矩阵的计算速度更快，因为MATLAB只对非零元素进行操作，这是稀疏矩阵的第二 个突出的优点。 ■例9.2 假设矩阵A、B和例9.1中的矩阵一样。计算2*A需要一百万次的浮点运算，而计算2*B只 需要2000次浮点运算。 ■ 因为MATLAB不能自动创建稀疏矩阵，所以要用特殊的命令来得到稀疏矩阵，在下一节 中将给出这些命令。前面章节中的算术和逻辑运算都适用于稀疏矩阵。 9.2创建和转换稀疏矩阵 在MATLAB中，用命令sparse来创建一个稀疏矩阵。 命令集87 创建稀疏矩阵 sparse(A) 由非零元素和下标建立稀疏矩阵A。如果A已是一个稀疏 矩阵，则返回A本身。 sparse(m,n) 生成一个m×n的所有元素都是0的稀疏矩阵

下载 第9章 稀 疏 矩 阵 在许多问题中提到了含有大量 0元素的矩阵，这样的矩阵称为稀疏矩阵。比如求解普通或 者部分微分方程的数值解。为了节省存储空间和计算时间， M AT L A B考虑到矩阵的稀疏性， 在对它运算时有特殊的命令。 9.1 矩阵为什么稀疏 一个稀疏矩阵中有许多元素等于零，这便于矩阵的计算和保存。如果 M AT L A B把一个矩 阵当作稀疏矩阵，那么只需在 m×3的矩阵中存储 m个非零项。第1列是行下标，第 2列是列下 标，第3列是非零元素值，不必保存零元素。如果存储一个浮点数要 8个字节，存储每个下标 要4个字节，那么整个矩阵在内存中存储需要 1 6×m个字节。 ■ 例9 . 1 A = e y e ( 1 0 0 0 ) ; 得到一个1 0 0 0×1 0 0 0的单位矩阵，存储它需要8 Mb空间。如果使用命令： B = s p e y e ( 1 0 0 0 ) ; 用一个 1 0 0 0×3的矩阵来代表，每行包含有一个行下标、列下标和元素本身。现在只需 1 6 K b的空间就可以存储1 0 0 0×1 0 0 0的单位矩阵，它只需要满单位矩阵的 0 . 2 %存储空间。对于 许多的广义矩阵也可这样来作。 稀疏矩阵的计算速度更快，因为 M AT L A B只对非零元素进行操作，这是稀疏矩阵的第二 个突出的优点。 ■ 例9 . 2 假设矩阵A、B和例9 . 1中的矩阵一样。计算 2*A需要一百万次的浮点运算，而计算 2*B只 需要2 0 0 0次浮点运算。 因为M AT L A B不能自动创建稀疏矩阵，所以要用特殊的命令来得到稀疏矩阵，在下一节 中将给出这些命令。前面章节中的算术和逻辑运算都适用于稀疏矩阵。 9.2 创建和转换稀疏矩阵 在M AT L A B中，用命令s p a r s e来创建一个稀疏矩阵。 命令集8 7 创建稀疏矩阵 s p a r s e ( A ) 由非零元素和下标建立稀疏矩阵 A。如果A已是一个稀疏 矩阵，则返回A本身。 s p a r s e ( m , n ) 生成一个m×n的所有元素都是0的稀疏矩阵。 ■ ■

124 China-pub.com MATLAB5手册 下载 sparse (u,v,a) 生成一个由长度相同的向量u,v和a定义的稀疏矩阵。其 中u和v是整数向量，a是一个实数或者复数向量。(4，v) 对应值a,如果a中有零元素，则将这个元素排除在外。 稀疏矩阵的大小为max(m×max(v)。 sparse (u,v,a,m,n) 生成一个mXn的稀疏矩阵，(w,)对应值a。向量u,v和a 必须长度相同。 sparse (u,v,a,m,n, 生成一个m×n的含有nzmax个非零元素的稀疏矩阵。（山， nzmax)） )对应值a。nzmax的值必须大于或者等于向量u和v的长度。 find(x) 返回向量x中非零元素的下标。如果x=X是一个矩阵，那 么X的向量就作为一个长向量来考虑。 [u,v]=find(A) 返回矩阵A中非零元素的下标。 [u,v,s]=find(A) 返回矩阵A中非零元素的下标。用向量中元素的值及u和v中 相应的下标，实际上就是向量、v和s作为命令sparsel的参数。 spconvert(D) 将一个有三列的矩阵转换成一个稀疏矩阵。D中的第1列作 为行的下标，第2列作为列的下标，最后一列作为元素值。 而且可以使用命令f将稀疏矩阵转换成一个满矩阵。 命令集88 转换成满矩阵 full(S) 将稀疏矩阵S转换成一个满矩阵。 ■例9.3 (a)创建一个5×5的单位矩阵： A=eye(5) 将矩阵A转换成稀疏矩阵B: B sparse(A) B= (1,1) 1 (2,2) 1 (3,3) 1 (4,4) 1 (5,5) 1 (b)假设MATLAB中给出如下的向量： ind1=[123342]; 1nd2=[121453]; number=[012305]; 这样就有了行向量，但是也可使用列向量。运行命matrix=sparse(indl,ind2,number), 结果为： Smatrix (3,1) 2

s p a r s e ( u , v , a ) 生成一个由长度相同的向量 u，v和a定义的稀疏矩阵。其 中u和v是整数向量， a是一个实数或者复数向量。 (ui, vi) 对应值 ai，如果 a中有零元素，则将这个元素排除在外。 稀疏矩阵的大小为m a x (u)×m a x (v)。 s p a r s e ( u , v , a , m , n ) 生成一个m×n的稀疏矩阵，(ui, vi)对应值ai。向量u，v和a 必须长度相同。 s p a r s e ( u , v , a , m , n , 生成一个m×n的含有n z m a x个非零元素的稀疏矩阵。(ui, n z m a x ) vi)对应值ai。n z m a x的值必须大于或者等于向量u和v的长度。 f i n d ( x ) 返回向量x中非零元素的下标。如果 x=X是一个矩阵，那 么X的向量就作为一个长向量来考虑。 [ u , v ] = f i n d ( A ) 返回矩阵A中非零元素的下标。 [ u , v , s ] = f i n d ( A ) 返回矩阵A中非零元素的下标。用向量s中元素的值及u和v中 相应的下标，实际上就是向量u、v和s作为命令s p a r s e的参数。 s p c o n v e r t ( D ) 将一个有三列的矩阵转换成一个稀疏矩阵。 D中的第1列作 为行的下标，第2列作为列的下标，最后一列作为元素值。 而且可以使用命令f u l l将稀疏矩阵转换成一个满矩阵。 命令集8 8 转换成满矩阵 f u l l ( S ) 将稀疏矩阵S转换成一个满矩阵。 ■ 例9 . 3 (a) 创建一个5×5的单位矩阵： A = e y e ( 5 ) 将矩阵A转换成稀疏矩阵B： (b) 假设M AT L A B中给出如下的向量： 这样就有了行向量，但是也可使用列向量。运行命令S m a t r i x = s p a r s e ( i n d 1 , i n d 2 , n u m b e r )， 结果为： 1 2 4 M ATLAB 5 手册 下载

China-bub.coM 第9章稀硫矩阵 125 下载 (2,2) 1 (2,3) (3,4) 3 其中有去掉了两个零元素。将这个矩阵转换成满矩阵，输入： Fullmatrix=full(Smatrix) 得到的结果为： Fullmatrix 0 0 0 0 0 1 5 0 0 2 0 0 3 0 0 0 0 0 注意，稀疏矩阵和得到的满矩阵的大小是分别是由ind1和ind2中最大元素值确定的，即 使相应的值是零，并且在列出的稀疏矩阵中去掉这个值。 输入命令whos可得到： Name Size Bytes Class A 5x5 200 double array 5x5 84 sparse array Fullmatrix 4x5 160 double array Smatrix 4x5 96 sparse array ind1 1x6 48 double array ind2 1x6 48 double array number 1x6 48 double array Grand total is 74 elements using 684 bytes 可以看出虽然两个矩阵的大小相同，但是其中稀疏矩阵需要的存储空间更小些。 (c)在处理稀疏矩阵时find命令很有用。命令对于稀疏矩阵或者满矩阵都返回相同的结果。 返回得到的三个向量直接用来重新创建一个稀疏矩阵。令Smatrix定义在(b)中，运行命令： [ind1,ind2,number]find(Smatrix); Smaller sparse(ind1,ind2,number) 得到的结果为： Smaller (3,1) 2 (2,2) 1 (2,3) 5 (3,4) 3 用下面命令得到的矩阵和(b)中得到的矩阵是不一样的： Fullsmall full(Smaller) Fullsmall 0 0 0 O 0 0 2 ■

其中有去掉了两个零元素。将这个矩阵转换成满矩阵，输入： F u l l m a t r i x = f u l l ( S m a t r i x ) 得到的结果为： 注意，稀疏矩阵和得到的满矩阵的大小是分别是由 i n d 1和i n d 2中最大元素值确定的，即 使相应的值是零，并且在列出的稀疏矩阵中去掉这个值。 输入命令w h o s可得到： 可以看出虽然两个矩阵的大小相同，但是其中稀疏矩阵需要的存储空间更小些。 (c) 在处理稀疏矩阵时f i n d命令很有用。命令对于稀疏矩阵或者满矩阵都返回相同的结果。 返回得到的三个向量直接用来重新创建一个稀疏矩阵。令 S m a t r i x定义在( b )中，运行命令： 得到的结果为： 用下面命令得到的矩阵和( b )中得到的矩阵是不一样的： 第9章 稀 疏 矩 阵 1 2 5 下载 ■

126 MATLAB5手册 China-pub.com 下载 9.3稀疏矩阵运算 MATLAB中对满矩阵的运算和函数同样可用在稀疏矩阵中。结果是稀疏矩阵还是满矩阵， 这取决于运算符或者函数及下列的操作数： ·当函数用一个矩阵作为输入参数，输出参数为一个标量或者一个给定大小的向量时，输 出参数的格式总是返回一个满阵形式，如命令size。 ·当函数用一个标量或者一个向量作为输入参数，输出参数为一个矩阵时，输出参数的格式 也总是返回一个满矩阵，如命令eye。还有一些特殊的命令可以得到稀疏矩阵，如命令speye。 ·对于单参数的其他函数来说，通常返回的结果和参数的形式是一样的，如diag。 ·对于双参数的运算或者函数来说，如果两个参数的形式一样，那么也返回同样形式的结 果。在两个参数形式不一样的情况下，除非运算的需要，均以满矩阵的形式给出结果。 ·两个矩阵的组和[AB],如果A或B中至少有一个是满矩阵，则得到的结果就是满矩阵。 ·表达式右边的冒号是要求一个参数的运算符，遵守这些运算规则。 ·表达式左边的冒号不改变矩阵的形式。 ■例9.4 假设有： A=eye(5);B=sparse(A);h=[1;2;0;4;5]; 这是一个5×5的单位满矩阵和相应的稀疏矩阵。 (a)C=5*B,结果为： C= (1,1) 5 (2,2) 5 (3,3) 6 (4,4) (5,5) 6 这是一个稀疏矩阵。 (b)D=A+B,给出的结果为： D= 2 0 0 0 2 0 0 0 0 0 2 0 0 0 0 0 2 0 0 O 0 2 这是一个满矩阵。 (c)=B\h,结果为： X= 1 2 4 5 ■ 这是一个满向量

9.3 稀疏矩阵运算 M AT L A B中对满矩阵的运算和函数同样可用在稀疏矩阵中。结果是稀疏矩阵还是满矩阵， 这取决于运算符或者函数及下列的操作数： • 当函数用一个矩阵作为输入参数，输出参数为一个标量或者一个给定大小的向量时，输 出参数的格式总是返回一个满阵形式，如命令 s i z e。 • 当函数用一个标量或者一个向量作为输入参数，输出参数为一个矩阵时，输出参数的格式 也总是返回一个满矩阵，如命令e y e。还有一些特殊的命令可以得到稀疏矩阵，如命令s p e y e。 • 对于单参数的其他函数来说，通常返回的结果和参数的形式是一样的，如 d i a g。 • 对于双参数的运算或者函数来说，如果两个参数的形式一样，那么也返回同样形式的结 果。在两个参数形式不一样的情况下，除非运算的需要，均以满矩阵的形式给出结果。 • 两个矩阵的组和[A B]，如果A或B中至少有一个是满矩阵，则得到的结果就是满矩阵。 • 表达式右边的冒号是要求一个参数的运算符，遵守这些运算规则。 • 表达式左边的冒号不改变矩阵的形式。 ■ 例9 . 4 假设有： 这是一个5×5的单位满矩阵和相应的稀疏矩阵。 (a) C = 5*B，结果为： 这是一个稀疏矩阵。 (b) D = A + B，给出的结果为： 这是一个满矩阵。 (c) x = B \ h，结果为： 这是一个满向量。 1 2 6 M ATLAB 5 手册 下载 ■

China-pub.com 第9章稀硫矩阵 127 下载 有许多命令可以对非零元素进行操作。 命令集89 矩阵的非零元素 nnz(A) 求矩阵A中非零元素的个数。它既可求满矩阵也可求稀疏矩阵。 spy(A) 画出稀疏矩阵A中非零元素的分布。也可用在满矩阵中，在 这种情况下，只给出非零元素的分布。 sPy(A,cstr,s1ze)用指定的颜色cstr(见表13-I)和在size规定的范围内画出稀疏 矩阵A中非零元素的分布。 nonzeros(A) 按照列的顺序找出矩阵A中非零的元素。 spones(A) 把矩阵A中的非零元素全换为1。 spalloc(m,n, 产生一个mXn阶只有nzmax个非零元素的稀疏矩阵。这样可以 nzmax) 有效地减少存储空间和提高运算速度。 nzmax(A) 给出为矩阵A中非零元素分配的内存数。不一定和nnz(A)得 到的数相同：参见sparse或者spalloc。 issparse(A) 如果矩阵A是稀疏矩阵，则返回1：否则返回0。 spfun(fcn,A) 用A中所有非零元素对函数fcn求值，如果函数不是对稀疏矩 阵定义的，同样也可以求值。 spfun(A) 求稀疏矩阵A的结构秩。对于所有的矩阵来说，都有 sprank(A≥rank(A)。 ■例9.5 用下面的命令定义稀疏矩阵： A=sparse(diag(ones(5,1),1))+sparse(diag(ones(5,1),-1)); 现在创建一个大矩阵： Big=kron(A,A) 这个矩阵Big是什么样子呢？Kronecker张量积给出一个大矩阵，它的元素是矩阵A的元素 之间可能的乘积。因为参量都是稀疏矩阵，所以得到的矩阵也是一个稀疏矩阵。可以用命令 whos和issparse来确认一下。 查看矩阵Big的结构图，可输入spy(Big),结构如图9-I所示。 : :8 ∷ ∷： 10 25、30 图91用spy命令显示的矩阵结构图

有许多命令可以对非零元素进行操作。 命令集8 9 矩阵的非零元素 n n z ( A ) 求矩阵A中非零元素的个数。它既可求满矩阵也可求稀疏矩阵。 s p y ( A ) 画出稀疏矩阵 A中非零元素的分布。也可用在满矩阵中，在 这种情况下，只给出非零元素的分布。 s p y ( A , c s t r , s i z e ) 用指定的颜色 c s t r(见表1 3 - 1 )和在s i z e规定的范围内画出稀疏 矩阵A中非零元素的分布。 n o n z e r o s ( A ) 按照列的顺序找出矩阵A中非零的元素。 s p o n e s ( A ) 把矩阵A中的非零元素全换为1。 s p a l l o c ( m , n , 产生一个m×n阶只有n z m a x个非零元素的稀疏矩阵。这样可以 n z m a x ) 有效地减少存储空间和提高运算速度。 n z m a x ( A ) 给出为矩阵A中非零元素分配的内存数。不一定和 n n z ( A )得 到的数相同；参见s p a r s e或者s p a l l o c。 i s s p a r s e ( A ) 如果矩阵A是稀疏矩阵，则返回1；否则返回0。 s p f u n ( f c n , A ) 用A中所有非零元素对函数 f c n求值，如果函数不是对稀疏矩 阵定义的，同样也可以求值。 s p f u n ( A ) 求稀疏矩阵 A的结构秩。对于所有的矩阵来说，都有 s p r a n k ( A≥r a n k ( A )。 ■ 例9 . 5 用下面的命令定义稀疏矩阵： 现在创建一个大矩阵： Big=kron(A, A) 这个矩阵B i g是什么样子呢？K r o n e c k e r张量积给出一个大矩阵，它的元素是矩阵 A的元素 之间可能的乘积。因为参量都是稀疏矩阵，所以得到的矩阵也是一个稀疏矩阵。可以用命令 w h o s和i s s p a r s e来确认一下。 查看矩阵B i g的结构图，可输入s p y ( B i g )，结构如图9 - 1所示。 图9-1 用s p y命令显示的矩阵结构图 第9章 稀 疏 矩 阵 1 2 7 下载

128 MATLAB5手册 China-pub.com 下载 可以看出Bg是一个块双对角矩阵。 9.4 稀疏矩阵的特例 MATLAB中有四个基本稀疏矩阵，它们是单位矩阵、随机矩阵、对称随机矩阵和对角矩阵。 命令集90 单位稀疏矩阵 speye(n) 生成n×n的单位稀疏矩阵。 speye(m,n) 生成m×n的单位稀疏矩阵。 命令speyel(A)得到的结果和sparse(eye(A))是一样的，但是没有涉及到满阵的存储。 命令集91 随机稀疏矩阵 sprand(A) 生成与A有相同结构的随机稀疏矩阵，且元素服从均匀分布。 sprand(m,n,dens) 生成一个m×n的服从均匀分布的随机稀疏矩阵，有es×mX n个非零元素，0≤des≤l。参数des是非零元素的分布密度。 sprand (m,n,dens, 生成一个近似的条件数为1/rc、大小为m×n的随机稀疏矩 rc) 阵。如果rc=rc是一个长度为l≤I(min(m,n)的向量，那么 矩阵将r心作为它1个奇异值的第一个，其他的奇异值为0。 sprandn(A) 生成与A有相同结构的随机稀疏矩阵，且元素服从正态分布。 sprandn (m,n,dens, 生成一个m×n的服从正态分布的随机稀疏矩阵，和sprand rc) 一样。 sprandsym(S) 生成一个随机对称稀疏矩阵。它的下三角及主对角线部分 与S的结构相同，矩阵元素服从正态分布。 sprandsym(n,dens) 生成一个m×n的随机对称稀疏矩阵。矩阵元素服从正态 分布，分布密度为dens。 sprandsym(n,dens, 生成一个近似条件数为1rc的随机对称稀疏矩阵。元素以0 rc) 对称分布，但不是正态分布。如果c=rc是一个向量，则矩 阵有特征值rC,。也就是说，如果rc是一个正向量，则矩阵 是正定矩阵。 sprandsym(n,dens, 生成一个正定矩阵。如果=1，则矩阵是由一正定对称矩 rc,k) 阵经随机Jacobij旋转得到的，其条件数正好等于l/rc:如 果k=2,则矩阵为外积的换位和，其条件数近似等于1/c。 sprandsym(S,dens, 生成一个与矩阵S结构相同的稀疏矩阵，近似条件数为l/rc。 rC,3) 参数dens被忽略，但是这个参数在这个位置以便函数能确认 最后两个参数的正确与否。 ■例9.6 (a)假设有矩阵A:

可以看出B i g是一个块双对角矩阵。 9.4 稀疏矩阵的特例 MATLAB中有四个基本稀疏矩阵，它们是单位矩阵、随机矩阵、对称随机矩阵和对角矩阵。 命令集9 0 单位稀疏矩阵 s p e y e ( n ) 生成n×n的单位稀疏矩阵。 s p e y e ( m , n ) 生成m×n的单位稀疏矩阵。 命令speye(A) 得到的结果和s p a r s e ( e y e ( A ) )是一样的，但是没有涉及到满阵的存储。 命令集9 1 随机稀疏矩阵 s p r a n d ( A ) 生成与A有相同结构的随机稀疏矩阵，且元素服从均匀分布。 s p r a n d ( m , n , d e n s ) 生成一个m×n的服从均匀分布的随机稀疏矩阵，有d e n s×m× n个非零元素，0≤d e n s≤1。参数d e n s是非零元素的分布密度。 s p r a n d ( m , n , d e n s , 生成一个近似的条件数为1 /rc、大小为m×n的随机稀疏矩 r c ) 阵。如果rc=rc是一个长度为l≤l ( m i n (m, n) )的向量，那么 矩阵将rci作为它l个奇异值的第一个，其他的奇异值为 0。 s p r a n d n ( A ) 生成与A有相同结构的随机稀疏矩阵，且元素服从正态分布。 s p r a n d n ( m , n , d e n s , 生成一个m×n的服从正态分布的随机稀疏矩阵，和sprand r c ) 一样。 s p r a n d s y m ( S ) 生成一个随机对称稀疏矩阵。它的下三角及主对角线部分 与S的结构相同，矩阵元素服从正态分布。 s p r a n d s y m ( n , d e n s ) 生成一个 m×n的随机对称稀疏矩阵。矩阵元素服从正态 分布，分布密度为d e n s。 s p r a n d s y m ( n , d e n s , 生成一个近似条件数为1 /rc的随机对称稀疏矩阵。元素以0 r c ) 对称分布，但不是正态分布。如果 rc=rc是一个向量，则矩 阵有特征值rci。也就是说，如果 rc是一个正向量，则矩阵 是正定矩阵。 s p r a n d s y m ( n , d e n s , 生成一个正定矩阵。如果k= 1，则矩阵是由一正定对称矩 r c , k ) 阵经随机J a c o b i旋转得到的，其条件数正好等于 1 /rc；如 果k= 2，则矩阵为外积的换位和，其条件数近似等于 1 /rc。 s p r a n d s y m ( S , d e n s , 生成一个与矩阵S结构相同的稀疏矩阵，近似条件数为1 /rc。 r c , 3 ) 参数d e n s被忽略，但是这个参数在这个位置以便函数能确认 最后两个参数的正确与否。 ■ 例9 . 6 (a) 假设有矩阵A： 1 2 8 M ATLAB 5 手册 下载 ■

China-pub.com 第9章稀硫矩阵 129 下载 /0100 1000 A= 0100 0010/ 输入Random=sprandn(A),可得到随机稀疏矩阵： Random (2,1) -0.4326 (1,2) -1.6656 (3,2) 0.1253 (4,3) 0.2877 矩阵中随机数的位置和矩阵A中非零元素的位置相同。 (b)对于(a)中的矩阵A,输入： B=sprandsym(A) 结果为： B= (2,1) -1.1465 (1,2) -1.1465 (3,2) 1.1909 (2,3) 1.1909 (4,3) 1.1892 (3,4) 1.1892 这是一个用矩阵A的下三角及主对角线部分创建的对称矩阵，在非零元素的位置用随机数 作为元素值。 ■ 用命令spdiags可以取出对角线元素，并创建带状对角矩阵。假设矩阵A的大小为mXn, 在p个对角线上有非零元素。B的大小为min(mXn)Xp,它的列是矩阵A的对角线。向量d的长 度为P,其整型分量给定了A的对角元： d0主对角线上的对角线 命令集92 对角稀疏矩阵 [B,d]=spdiags(A) 求出A中所有的对角元，对角元保存在矩阵B中，它们的 下标保存在向量d中。 spdiags(A,d) 生成一个矩阵，这个矩阵包含有矩阵中向量规定的对角元。 spdiags(B,d,A) 生成矩阵A,用矩阵B中的列替换d定义的对角元。 A=spdiags(B,d,m,n)用保存在由d定义的B中的对角元创建稀疏矩阵A。 例11.4给出了如何使用spdiags命令来解普通微分方程组。 9.5系数阵为稀疏矩阵的线性方程组 在许多实际应用中要保留稀疏矩阵的结构，但是在计算过程中的中间结果会减弱它的稀 疏性，如LU分解。这就会导致增加浮点运算次数和存储空间。为了避免这种情况发生，在

输入R a n d o m = s p r a n d n ( A )，可得到随机稀疏矩阵： 矩阵中随机数的位置和矩阵 A中非零元素的位置相同。 (b) 对于( a )中的矩阵A，输入： B = s p r a n d s y m ( A ) 结果为： 这是一个用矩阵A的下三角及主对角线部分创建的对称矩阵，在非零元素的位置用随机数 作为元素值。 用命令s p d i a g s可以取出对角线元素，并创建带状对角矩阵。假设矩阵 A的大小为m×n， 在p个对角线上有非零元素。 B的大小为m i n (m×n)×p，它的列是矩阵A的对角线。向量d的长 度为p，其整型分量给定了A的对角元： di0 主对角线上的对角线 命令集9 2 对角稀疏矩阵 [ B , d ] = s p d i a g s ( A ) 求出A中所有的对角元，对角元保存在矩阵 B中，它们的 下标保存在向量d中。 s p d i a g s ( A , d ) 生成一个矩阵，这个矩阵包含有矩阵A中向量d规定的对角元。 s p d i a g s ( B , d , A ) 生成矩阵A，用矩阵B中的列替换d定义的对角元。 A = s p d i a g s ( B , d , m , n ) 用保存在由d定义的B中的对角元创建稀疏矩阵A。 例11 . 4给出了如何使用s p d i a g s命令来解普通微分方程组。 9.5 系数阵为稀疏矩阵的线性方程组 在许多实际应用中要保留稀疏矩阵的结构，但是在计算过程中的中间结果会减弱它的稀 疏性，如 L U分解。这就会导致增加浮点运算次数和存储空间。为了避免这种情况发生，在 第9章 稀 疏 矩 阵 1 2 9 下载 ■

130 China-pub.com MATLAB5手册 下载 MATLAB中用命令对矩阵进行重新安排。这些命令都列在下面的命令集93中。通过help命令 可以得到每个命令更多的帮助信息，也可见helpdesk。 命令集93 矩阵变换 colmmd(A) 返回一个变换向量，使得矩阵A列的秩为最小。 symmmd (A) 返回使对称矩阵秩为最小的变换。 symrcm(A) 矩阵A的Cuthill--McKee逆变换。矩阵A的非零元素在主对角线附近。 colperm(A) 返回一个矩阵A的列变换的向量。列按非零元素升序排列。有时这 是LU因式分解前有用的变换：1uA(:, )。如果A是一个对称矩 阵，对行和列进行排序，这有利Cholesky分解：dola5,) randperm(n) 给出正数1,2，，的随机排列，可以用来创建随机变换矩阵。 dmperm(A) 对矩阵A进行Dulmage-Mendelsohn分解，输入help dmper可 得更多信息。 ■例9.7 创建一个秩为4的变换矩阵，可输入： i =[1 2 3 4];aa ones(1,4);perm randperm(4) P=sparse(i,perm,aa) 一旦运行perm=randperm(4),就会得到： perm 3 给出的变换矩阵为： p= (2,1) 1 (4,2) 1 (3,3) 1 (1,4) 1 如果矩阵A为： 7 4 3 5 2 A= 63 1 8 1 1 2 输入命令： RowChange P*A,ColChange A*P 运行结果为： RowChange 5 6 3 1 5 2 4 2 1 1 2 7 4 3 4

M AT L A B中用命令对矩阵进行重新安排。这些命令都列在下面的命令集 9 3中。通过h e l p命令 可以得到每个命令更多的帮助信息，也可见 h e l p d e s k。 命令集9 3 矩阵变换 c o l m m d ( A ) 返回一个变换向量，使得矩阵 A列的秩为最小。 s y m m m d ( A ) 返回使对称矩阵秩为最小的变换。 s y m r c m ( A ) 矩阵A的C u t h i l l - M c K e e逆变换。矩阵A的非零元素在主对角线附近。 c o l p e r m ( A ) 返回一个矩阵A的列变换的向量。列按非零元素升序排列。有时这 是L U因式分解前有用的变换：lu(A(:, j))。如果A是一个对称矩 阵，对行和列进行排序，这有利于C h o l e s k y分解：chol(A(j, j)) 。 r a n d p e r m ( n ) 给出正数1，2，. . .，n的随机排列，可以用来创建随机变换矩阵。 d m p e r m ( A ) 对矩阵A进行D u l m a g e - M e n d e l s o h n分解，输入help dmperm可 得更多信息。 ■ 例9 . 7 创建一个秩为4的变换矩阵，可输入： 一旦运行p e r m = r a n d p e r m ( 4 )，就会得到： 给出的变换矩阵为： 如果矩阵A为： 输入命令： 运行结果为： 1 3 0 M ATLAB 5 手册 下载

China-pub.com 第9章稀硫矩阵 131 下载 ColChange 4 又 3 2 4 1 6 3 2 1 8 1 ■ 有两个不完全因式分解命令，它们是用来在解大线性方程组前进行预处理的。用 helpdesk命令可得更多信息。 命令集94 不完全因式分解 cho1inc(A,opt)进行不完全Cholesky分解，变量opt取下列值之一： dropto1指定不完全分解的舍入误差，O给出完全分解。 michol如果nichol=l,则从对角线上抽取出被去掉的元素。 rdiag用sqrt(dropto1*norm(x(:,j)))代替上三角分 解因子中的零元素，为零元素所在的列。 [L,0,P]= 返回矩阵X的不完全分解得到的三个矩阵L、U和P,变量opt取 luinc(X,opt) 下列值之一： droptol指定分解的舍入误差。 milu 改变分解以便从上三角角分解因子中抽取被去掉的列元素。 udiag 用dropto1值代替上三角角分解因子中的对角线上的 零元素。 thresh中心极限。 解稀疏线性方程组既可用左除运算符解，也可用一些特殊命令来解。 命令集95 稀疏矩阵和线性方程组 spparms(keystr,op)设置稀疏矩阵算法的参数，用help spparms可得详细信息。 spaugment (A,c) 根据[c*!A,A可创建稀疏矩阵，这是二次线性方程组的最 小二乘问题。参见7.7节。 symbfact (A) 给出稀疏矩阵的Cholesky和LU因式分解的符号分解因子。 用help symbfad可得详细信息。 稀疏矩阵的范数计算和普通满矩阵的范数计算有一个重要的区别。稀疏矩阵的欧几里德 范数不能直接求得。如果稀疏矩阵是一个小矩阵，则用norm(ful1(A))来计算它的范数： 但是对于大矩阵来说，这样计算是不可能的。然而MATLAB可以计算出欧几里德范数的近似 值，在计算条件数时也是一样。 命令集96 稀疏矩阵的近似欧几里德范数和条件数 normest (A) 计算A的近似欧几里德范数，相对误差为10-6。 normest (A,tol) 计算A的近似欧几里德范数，设置相对误差tOl,而不用缺省时 的10-6

有两个不完全因式分解命令，它们是用来在解大线性方程组前进行预处理的。用 h e l p d e s k命令可得更多信息。 命令集9 4 不完全因式分解 c h o l i n c ( A , o p t ) 进行不完全C h o l e s k y分解，变量o p t取下列值之一： d r o p t o l 指定不完全分解的舍入误差， 0给出完全分解。 m i c h o l 如果m i c h o l = 1，则从对角线上抽取出被去掉的元素。 r d i a g 用s q r t ( d r o p t o l*n o r m ( X ( : , j ) ) )代替上三角分 解因子中的零元素，j为零元素所在的列。 [ L , U , P ]= 返回矩阵X的不完全分解得到的三个矩阵 L、U和P，变量o p t取 l u i n c ( X , o p t ) 下列值之一： d r o p t o l 指定分解的舍入误差。 m i l u 改变分解以便从上三角角分解因子中抽取被去掉的列元素。 u d i a g 用d r o p t o l值代替上三角角分解因子中的对角线上的 零元素。 t h r e s h 中心极限。 解稀疏线性方程组既可用左除运算符解，也可用一些特殊命令来解。 命令集9 5 稀疏矩阵和线性方程组 s p p a r m s ( k e y s t r , o p ) 设置稀疏矩阵算法的参数，用 help spparms可得详细信息。 s p a u g m e n t ( A , c ) 根据[ c*l A; A’ 0 ] 创建稀疏矩阵，这是二次线性方程组的最 小二乘问题。参见7 . 7节。 s y m b f a c t ( A ) 给出稀疏矩阵的 C h o l e s k y和L U因式分解的符号分解因子。 用help symbfact可得详细信息。 稀疏矩阵的范数计算和普通满矩阵的范数计算有一个重要的区别。稀疏矩阵的欧几里德 范数不能直接求得。如果稀疏矩阵是一个小矩阵，则用 n o r m ( f u l l ( A ) )来计算它的范数； 但是对于大矩阵来说，这样计算是不可能的。然而 M AT L A B可以计算出欧几里德范数的近似 值，在计算条件数时也是一样。 命令集9 6 稀疏矩阵的近似欧几里德范数和条件数 n o r m e s t ( A ) 计算A的近似欧几里德范数，相对误差为 1 0－6。 n o r m e s t ( A , t o l ) 计算A的近似欧几里德范数，设置相对误差 t o l，而不用缺省时 的1 0－6。 第9章 稀 疏 矩 阵 1 3 1 下载 ■

132 China-pub.CoM MATLAB5手册 下载 [nrm,nit]= 计算近似rm范数，还给出计算范数迭代的次数nit。 normest(A) condest(A) 求矩阵A条件数的1-范数中的下界估计值。 [c,v]= 求矩阵A的1-范数中条件数的下界估计值c和向量v,使得 condest (A,tr) IAv=(Al·Iv)/c。如果给定r,则给出计算的过程。=1， 给出每步过程：t=一I,给出商c/rcond(A)。 ■例9.8 假设给出： Sprs speye(4);Sprs(4,1)=19;Sprs(3,2)=4; 用normApprox=normest(Sprs)计算出： normApprox 19.0525 用theNorm=norm(full(Sprs)得： theNorm 19.0525 为了找到它们之间的差别，计算difference=theNorm-normApprox,得： difference= 8.5577e-09 在许多应用中，normest计算得到的近似值是一个很好的近似欧几里德范数，它的计算 步数要比norm要少得多：可参见7.6节。 ◆ 用etree命令来找到稀疏对称矩阵的消元树，用向量「来描述消元树，还可用 etreeplot命令画出来。元素f是矩阵的上三角Cholesky分解因子中i行上第l非元素的列下 标。如果有非零元素，则=0。消元树可以这样来建立： 节点是f的孩子，或者如果0，则节点是树的根节点。 命令集97 矩阵的消元树 etree(A) 求A的消元树向量f,这个命令有可选参数：输入help etree获取帮助。 etreeplot(A) 画出向量定义的消元树图形。 treeplot(p,c,d) 画出指针向量p的树图形，参数c和d分别指定节点的颜色 和分支数。etreeplot可以调用这个命令。 treelayout 显示树的结构，treeplot可以调用这个命令。 ■例9.9 假设有对称稀疏矩阵B:

[ n r m , n i t ] = 计算近似n r m范数，还给出计算范数迭代的次数 n i t。 n o r m e s t ( A ) c o n d e s t ( A ) 求矩阵A条件数的1 -范数中的下界估计值。 [ c , v ]= 求矩阵A的1 -范数中条件数的下界估计值 c和向量v，使得 c o n d e s t ( A , t r ) | |Av| | = ( | |A| |·| |v| | ) / c。如果给定 t r，则给出计算的过程。 t r= 1， 给出每步过程；t r=－1，给出商c / r c o n d ( A )。 ■ 例9 . 8 假设给出： 用n o r m A p p r o x = n o r m e s t ( S p r s )计算出： 用t h e N o r m = n o r m ( f u l l ( S p r s ) )得： 为了找到它们之间的差别，计算 d i f f e r e n c e = t h e N o r m - n o r m A p p r o x，得： 在许多应用中，n o r m e s t计算得到的近似值是一个很好的近似欧几里德范数，它的计算 步数要比n o r m要少得多；可参见7 . 6节。 用 e t r e e 命令来找到稀疏对称矩阵的消元树，用向量 f 来描述消元树，还可用 e t r e e p l o t命令画出来。元素fi是矩阵的上三角C h o l e s k y分解因子中i行上第1非零元素的列下 标。如果有非零元素，则fi= 0。消元树可以这样来建立： 节点i是fi的孩子，或者如果fi= 0，则节点i是树的根节点。 命令集9 7 矩阵的消元树 e t r e e ( A ) 求A的消元树向量 f，这个命令有可选参数；输入 h e l p e t r e e获取帮助。 e t r e e p l o t ( A ) 画出向量f定义的消元树图形。 t r e e p l o t ( p , c , d ) 画出指针向量 p的树图形，参数 c和d分别指定节点的颜色 和分支数。e t r e e p l o t可以调用这个命令。 t r e e l a y o u t 显示树的结构，t r e e p l o t可以调用这个命令。 ■ 例9 . 9 假设有对称稀疏矩阵B： 1 3 2 M ATLAB 5 手册 下载 ■