|
7 稀疏矩阵 稀疏矩阵是一种特殊类型的矩阵，即矩阵中包括较多的零元素。对于稀疏矩阵的这种特性，在MATLAB中可以只保存矩阵中非零元素及非零元素在矩阵中的位置。在用稀疏矩阵进行计算时，通过消去零元素可以减少计算的时间。 7.1 稀疏矩阵的存储方式 对一般矩阵而言，MATLAB保存矩阵内的每一个元素，矩阵中的零元素与其他元素一样，需要占用同样大小的内存空间。但对于稀疏矩阵，MATLAB仅存储稀疏矩阵中的非零元素及其对应的位置，其他空余位置只是在访问时以默认的零元素来填充。对于一个含有大量零元素的大型矩阵，采用这种方法可以大大地减少数据占据的内存空间。 MATLAB采用3个内部数组来保存元素为实数的稀疏矩阵。 稀疏矩阵也可用于存储复数。当稀疏矩阵用于存储复数数据时，需用第4个内部数组保存非零复数的虚部。一个复数非零，是指其实部或虚部至少其中一个不为零。 【例2-16】 稀疏矩阵与一般矩阵的内存占用对比。 >> M_full = magic(1100); % 创建一个1100´1100 矩阵 >> M_full(M_full > 50) = 0; % 将>50的元素设置为0 >> M_sparse = sparse(M_full); % 创建稀疏矩阵 >> whos Name Size Bytes Class Attributes M_full 1100x1100 9680000 double M_sparse 1100x1100 9608 double sparse 本例中，M_full和 M_sparse两个变量存储的实际上是同一个矩阵，但是二者因为采用的存储形式分别为一般矩阵和稀疏矩阵，所以占用的内存量却相差了近1000倍。因为MATLAB版本不同，操作系统不同（例如32位和64位），内部存储格式也有些变化，但总体上来说占用的内存空间比一般矩阵小很多。 7.2 稀疏矩阵的创建 MATLAB决不会自动地创建一个稀疏矩阵，这需要用户来决定是否使用稀疏矩阵。在创建一个矩阵前，用户需要根据此矩阵中是否包含较多的零元素，采用稀疏矩阵技术是否有利，来决定是否采用稀疏矩阵的形式。把矩阵中非零元素的个数除以所有元素的个数，就叫做矩阵的密度，密度越小的矩阵采用稀疏矩阵的格式越有利。 要将一般矩阵转换为稀疏矩阵，可以使用函数sparse，如s=sparse (A)，是指将矩阵A转换为稀疏矩阵。另外，使用函数full则可把稀疏矩阵转换为一般矩阵。 【例2-17】 一般矩阵与稀疏矩阵的转换示例。 >> A=[0 0 0 1;0 1 0 0;1 2 0 0;0 0 3 0] A = 0 0 0 1 0 1 0 0 1 2 0 0 0 0 3 0 >> s=sparse(A) s = (3,1) 1 (2,2) 1 (3,2) 2 (4,3) 3 (1,4) 1 >> B=full(s) B = 0 0 0 1 0 1 0 0 1 2 0 0 0 0 3 0 从本例的结果中可以看出所有s的非零元素列表及其对应的行列序号。所有非零元素保存在一列中，反映了数据的内部结构。 稀疏矩阵的创建一般有以下几种方式。 1．直接创建稀疏矩阵 使用函数sparse，可以用一组非零元素直接创建一个稀疏矩阵。该函数调用格式为： S=sparse(i,j,s,m,n) 其中i和j都为矢量，分别是指矩阵中非零元素的行号与列号，s是一个全部为非零元素矢量，元素在矩阵中排列的位置为(i,j)。m为输出的稀疏矩阵的行数，n为输出的稀疏矩阵的列数。 【例2-18】 稀疏矩阵的创建。 >> S=sparse([1 3 2 1 4],[3 1 4 1 4],[1 2 3 45],4,4) S = (1,1) 4 (3,1) 2 (1,3) 1 (2,4) 3 (4,4) 5 >> full(S) ans = 4 0 1 0 0 0 0 3 2 0 0 0 0 0 0 5 本例中通过sparse函数直接创建了稀疏矩阵S。sparse函数中的前两个输入变量[1 3 2 1 4]和[3 1 4 1 4]就是元素在矩阵中排列的位置，第3个输入变量[1 2 3 4 5]就是稀疏矩阵前面两个输入变量中的位置所对应的元素的值，而最后的两个输入变量4是指输出的稀疏矩阵的行数是4，输出的稀疏矩阵的列数同样也为4。通过full函数把稀疏矩阵转换为一般矩阵，这样就可以清楚地看出sparse函数输入和输出之间的关系。 需要指出的是：函数sparse还有一个变化形式，可以设置最大数目的非零元素。如有必要，可以在函数sparse中添加第6个输入参数，设置稀疏矩阵中非零元素的最大数目，这样以后要在矩阵中添加非零元素，就无需再修改矩阵的结构。具体的使用方法请查阅help文档。 2．从对角线元素中创建稀疏矩阵 要将一个矩阵的对角线元素保存在一个稀疏矩阵中，可以使用函数spdiags。其调用格式为： S=spdiags(B,d,m,n) 函数spdiags用于创建一个尺寸为m行n列的稀疏矩阵S，其非零元素来自矩阵B中的元素且按对角线排列，参数d指定矩阵B中用于生成稀疏矩阵B的对角线位置。矩阵的主对角线可以认为是第0条对角线，每向右上移动一条对角线编号加1，向左下移动一条对角线编号减1，也就是说B中的j列填充矢量d元素，j指定对角线。 【例2-19】 稀疏矩阵的创建。 >> B=[1 2 3;4 5 6;7 8 9;10 11 12] B = 1 2 3 4 5 6 7 8 9 10 11 12 >> d=[-3 0 2] d = -3 0 2 >> A=spdiags(B,d,7,4) A = (1,1) 2 (4,1) 1 (2,2) 5 (5,2) 4 (1,3) 9 (3,3) 8 (6,3) 7 (2,4) 12 (4,4) 11 (7,4) 10 >> full(A) ans = 2 0 9 0 0 5 0 12 0 0 8 0 1 0 0 11 0 4 0 0 0 0 7 0 0 0 0 10 本例生成了一个7行4列的稀疏矩阵A。B的第1列元素排列在主对角线以下的第3条对角线上，第2列元素排列在主对角线上，第3列中的非零元素排列在主对角线上方的第2条对角线上。这里需要注意B中的第三列并没有全部分布在第2条对角线上，而是最后两个元素9和12排列在该对角线上。 3．从外部文件中导入稀疏矩阵 用外部文件创建的文本文件，如果该文件中的数据按3或者4列排列，则可将这个文本文件载入内存，用于创建一个稀疏矩阵。 【例2-20】 稀疏矩阵的创建。 假如有这样一个文件uphill.dat（用户可以通过记事本打开、编辑其内容），文件内含有以下数据： 1 1 1.000000000000000 1 2 0.500000000000000 2 2 0.333333333333333 1 3 0.333333333333333 2 3 0.250000000000000 3 3 0.200000000000000 1 4 0.250000000000000 2 4 0.200000000000000 3 4 0.166666666666667 4 4 0.142857142857143 4 4 0.000000000000000 那么通过使用load命令可以将此数据文件载入MATLAB，然后对其进行操作。实际中用户可以在命令窗口输入： >> load uphill.dat % 用load命令将数据的文本文件uphill.dat载入工作空间 >> H = spconvert(uphill) H = (1,1) 1.0000 (1,2) 0.5000 (2,2) 0.3333 (1,3) 0.3333 (2,3) 0.2500 (3,3) 0.2000 (1,4) 0.2500 (2,4) 0.2000 (3,4) 0.1667 (4,4) 0.1429 > >> full(H)
ans =
1.0000 0.5000 0.3333 0.2500 0 0.3333 0.2500 0.2000 0 0 0.2000 0.1667 0 0 0 0.1429 本例首先使用load函数导入了一个3列数据的文本文件uphill.dat，用户可以通过在命令行中输入变量名uphill 可以查看数据uphill 中的具体内容来验证数据读取是否正确，然后调用spconvert将uphill 转换为相应的稀疏矩阵H。通过调用full函数可以直观地看出得到的稀疏矩阵。 MATLAB使用load函数来导入外部数据文件，具体的用法可以参阅第10章。 7.3 稀疏矩阵的运算 多数MATLAB自带的数学函数都可用于处理稀疏矩阵，此时可以将稀疏矩阵当做一般矩阵看待。另外MATLAB也提供有一些专门针对稀疏矩阵的函数。处理稀疏矩阵时，计算的复杂程度与稀疏矩阵中非零元素的数目成正比，也与矩阵的行列尺寸有关。比如稀疏矩阵的乘法、乘方、包含一定次数的线性方程组等，都是比较复杂的运算。 用函数处理稀疏矩阵时，计算结果要遵循以下一些原则。 （1）MATLAB函数处理一个矩阵时，不管这个矩阵是一般矩阵还是稀疏矩阵，其返回值为一个数值或矩阵。返回值都按一般矩阵方式进行保存，并不会根据接受的参数是稀疏矩阵，而将结果保存为稀疏矩阵。 （2）函数处理一个数值或矢量返回一个矩阵时，如果矩阵为零矩阵、元素全为1的矩阵、随机矩阵或单位矩阵，这些矩阵全为一般矩阵形式。对于零矩阵，有一种类似稀疏矩阵的存储方法，因为零矩阵中没有非零元素，所以不能将一个零矩阵转换为一个稀疏矩阵，但指令zeros(m,n)和sparse(m,n)是可用的。对于单位矩阵和随机矩阵，可以使用类似稀疏矩阵的操作指令，即speye和sprand。对于元素全为1的矩阵，则没有类似的操作指令。 （3）以矩阵为参数返回矩阵或矢量的一元函数，返回值的存储类型与参数的存储类型相同。例如矩阵S的cholesky分解，如果S为一般矩阵，结果也为一般矩阵；如果S为稀疏矩阵，结果也为稀疏矩阵。对于列向处理矩阵的函数，如求各列最大值的函数max，求各列之和的函数sum等，也都返回与参数相同的存储类型。如果参数是稀疏矩阵，即使返回的矩阵或矢量全为非零元素，也用稀疏方式表示。例外情况只有函数sparse和full，因为它们用于一般矩阵和稀疏矩阵之间的转换。 （4）对于有两个输入参数为矩阵的情况，如果输入的两个矩阵都为稀疏矩阵，则输出仍为稀疏矩阵；都为一般矩阵，结果也为一般矩阵。如果输入参数一个为稀疏矩阵，一个为一般矩阵，结果通常为一般矩阵，但在能够保证矩阵稀疏性不变的运算中，结果则为稀疏矩阵。 （5）使用方括号对矩阵进行组合时，如果组合的矩阵中有稀疏矩阵，结果则为稀疏矩阵。 （6）子矩阵在右边的赋值操作，返回值为右边子矩阵的储存类型，子矩阵在左边赋值不改变其储存类型。 【例2-21】 稀疏矩阵的组合。 >> A=[1 0 0;0 0 1;1 2 0] A = 1 0 0 0 0 1 1 2 0 >> B=sparse(A) B = (1,1) 1 (3,1) 1 (3,2) 2 (2,3) 1 >> C=[A(:,1),B(:,2)] C = (1,1) 1 (3,1) 1 (3,2) 2 本例将矩阵A的第1列和矩阵B的第2列组成了新的矩阵C，从结果可知，C为稀疏矩阵。 【例2-22】 稀疏矩阵子矩阵的赋值。 >> A=[1 0 0;0 0 1;1 2 0]; >> B=sparse(A); >> C=sparse(cat(1,full(B),A)) C = (1,1) 1 (3,1) 1 (4,1) 1 (6,1) 1 (3,2) 2 (6,2) 2 (2,3) 1 (5,3) 1 >> i=[1 2 3]; >> j=[1 2 3]; >> T=C(i,j) T = (1,1) 1 (3,1) 1 (3,2) 2 (2,3) 1 >> C(j,i)=full(T) % 将一般矩阵赋值给一稀疏矩阵，仍返回稀疏矩阵 C = (1,1) 1 (3,1) 1 (4,1) 1 (6,1) 1 (3,2) 2 (6,2) 2 (2,3) 1 (5,3) 1 7.4 稀疏矩阵的交换与重新排序 稀疏矩阵S的行交换与列交换可以用以下两种方法表示。 （1）对于交换矩阵P，对稀疏矩阵S进行行交换可表示为P*S，进行列交换可表示为P*S’。 （2）对于一个交换矢量p，p为一般矢量包含1到n个自然数的一个排列。对稀疏矩阵进行行交换，可以表示为S(p,:)。S(:,p)为列交换形式。对于矩阵S的某一列进行行交换，可以表示为S(p,n)，如S(p,1)为对第1列进行交换。 【例2-23】 稀疏矩阵S的交换。 >> p=[1 3 2 4]; >> S=eye(4,4) S = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 >> P=S(p,:) P = 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 >> V=S(p,2) V = 0 0 1 0 矩阵P的第1行为S的第1行，第2行为S的第3行，等等。即对矩阵S的行，按照矢量p指定的顺序进行调整。 >> S1=speye(4,4) S1 = (1,1) 1 (2,2) 1 (3,3) 1 (4,4) 1 >> P1=S1(p,:) % 对于稀疏矩阵行列的交换，返回的形式仍为稀疏矩阵 P1 = (1,1) 1 (3,2) 1 (2,3) 1 (4,4) 1 对于稀疏矩阵S1进行行列的交换，返回的P1仍为稀疏矩阵。对稀疏矩阵的列重新排序，有时可以使矩阵分解的速度更快，最简单的矩阵排序是根据矩阵中非零元素的个数进行的，这种方法对于元素极不规则的矩阵很有效，特别适用于非零元素在行或列中数目变化较大的矩阵。MATLAB提供有一个非常简单的函数colperm，可以实现这种排序方法。此函数的M文件仅有以下几行： function p=colperm(S) if ~ismatrix(S) % 判断输入变量是否是一个矩阵 error(message('MATLAB:colperm:invalidInput')); % 不满足条件的话返回错误信息 end [~,p] = sort(full(sum(S ~= 0, 1))); 程序的第5行，实现了以下4个功能。 （1）调用S ~= 0判断矩阵中各元素是否为0，若不为0则返回逻辑值1。 （2）函数sum求上一步创建的矩阵各列的和，也即为各列中非零元素的个数。 （3）函数full将上一步创建的矢量转换为一般矢量的格式。 （4）使用函数sort对上一步操作创建的矢量元素进行升序排序，函数sort的第2个输出参数p，即为对矩阵S中各列中非零元素的个数进行重新排序的交换矢量。 【例2-24】 对下面的矩阵A，先用函数colperm获取一个交换矢量p，然后根据矢量p对矩阵A的列，按照非零元素的个数升序排序。 >> A=[0 1 2 3;3 2 1 0;0 0 2 0;1 0 0 2] A = 0 1 2 3 3 2 1 0 0 0 2 0 1 0 0 2 >> p=colperm(A) p = 1 2 4 3 >> B=A(:,p) B = 0 1 3 2 3 2 0 1 0 0 0 2 1 0 2 0 结果显示，矩阵B就是将A的各列按照非零元素的个数升序排序的结果。 7.5 稀疏矩阵视图 MATLAB提供有spy函数，用于观察稀疏矩阵非零元素的分布视图。本小节举例来说明spy函数的用法。 【例2-25】 稀疏矩阵视图示例。本例采用spy函数绘制Buckminster Fuller网格球顶的60×60邻接矩阵视图。这个矩阵还可用来表示碳60模型和足球。 >>B = bucky; >>spy(B) 得到的结果如图2-2所示。图中显示了稀疏矩阵B的非零元素分布视图。
|
|
|
