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