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              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]

0     0    0     1

0     1    0     0

1     2    0     0

0     0    3     0

>> s=sparse(A)

(3,1)        1

(2,2)        1

(3,2)        2

(4,3)        3

(1,4)        1

>> B=full(s)

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)

(1,1)        4

(3,1)        2

(1,3)        1

(2,4)        3

(4,4)        5

>> full(S)

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]

1     2    3

4     5    6

7     8    9

10    11   12

>> d=[-3 0 2]

-3     0    2

>> A=spdiags(B,d,7,4)

(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)

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.0000000

1    2    0.0000000

2    2    0.3333333

1    3    0.3333333

2    3    0.0000000

3    3    0.0000000

1    4    0.0000000

2    4    0.0000000

3    4    0.6666667

4    4    0.2857143

4    4    0.0000000

那么通过使用load命令可以将此数据文件载入MATLAB,然后对其进行操作。实际中用户可以在命令窗口输入:

>> load uphill.dat  %  用load命令将数据的文本文件uphill.dat载入工作空间

>> H = spconvert(uphill)

(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)

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]

1     0    0

0     0    1

1     2    0

>> B=sparse(A)

(1,1)        1

(3,1)        1

(3,2)        2

(2,3)        1

>> C=[A(:,1),B(:,2)]

(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))

(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)

(1,1)        1

(3,1)        1

(3,2)       2

(2,3)        1

>> C(j,i)=full(T)      %  将一般矩阵赋值给一稀疏矩阵,仍返回稀疏矩阵

(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的交换。

>> S=eye(4,4)

1     0    0     0

0     1    0     0

0     0    1     0

0     0    0     1

>> P=S(p,:)

1     0    0     0

0     0    1     0

0     1    0     0

0     0    0     1

>> V=S(p,2)

矩阵P的第1行为S的第1行,第2行为S的第3行,等等。即对矩阵S的行,按照矢量p指定的顺序进行调整。

>> S1=speye(4,4)

(1,1)        1

(2,2)        1

(3,3)        1

(4,4)        1

>> P1=S1(p,:)     %  对于稀疏矩阵行列的交换,返回的形式仍为稀疏矩阵

(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'));   % 不满足条件的话返回错误信息

[~,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的列,按照非零元素的个数升序排序。

0     1    2     3

3     2    1     0

0     0    2     0

1    0     0     2

>> p=colperm(A)

1     2    4     3

>> B=A(:,p)

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的非零元素分布视图。