Main Content

访问稀疏矩阵

非零元素

以下多条命令可以提供有关稀疏矩阵的非零元素的概要信息:

  • nnz 返回稀疏矩阵中的非零元素数。

  • nonzeros 返回包含稀疏矩阵的所有非零元素的列向量。

  • nzmax 返回为稀疏矩阵的非零项分配的存储空间量。

要尝试上述中的一些命令,请加载提供的稀疏矩阵 west0479,该矩阵是 Harwell-Boeing 集合之一。

load west0479
whos
  Name            Size             Bytes  Class     Attributes

  west0479      479x479            34032  double    sparse   

该矩阵为八个阶段的化工精馏塔建模。

尝试以下命令。

nnz(west0479)
ans =

        1887
format short e
west0479
west0479 =

  (25,1)      1.0000e+00
  (31,1)     -3.7648e-02
  (87,1)     -3.4424e-01
  (26,2)      1.0000e+00
  (31,2)     -2.4523e-02
  (88,2)     -3.7371e-01
  (27,3)      1.0000e+00
  (31,3)     -3.6613e-02
  (89,3)     -8.3694e-01
  (28,4)      1.3000e+02
     .
     .
     .
nonzeros(west0479)
ans =

   1.0000e+00
  -3.7648e-02
  -3.4424e-01
   1.0000e+00
  -2.4523e-02
  -3.7371e-01
   1.0000e+00
  -3.6613e-02
  -8.3694e-01
   1.3000e+02
    .
    .
    .

注意

使用 Ctrl+C 随时停止列出 nonzeros

请注意,最初在默认情况下,nnznzmax 的值相同。即,非零元素数等于为非零元素分配的存储位置数。但是,如果将其他的数组元素置零,MATLAB® 不会动态释放内存。将某些矩阵元素的值更改为零时会更改 nnz 的值,但不会更改 nzmax 的值。

但是,可以根据需要将尽可能多的非零元素添加到矩阵中。不受 nzmax 原始值的限制。

索引和值

对于任何矩阵,无论是满矩阵还是稀疏矩阵,find 函数都会返回非零元素的索引和值。其语法是

[i,j,s] = find(S);

find 返回向量 i 中的非零值的行索引、向量 j 中的列索引以及向量 s 中的自身非零值。下面的示例使用 find 查找稀疏矩阵中的非零索引和值。sparse 函数同时使用 find 输出和矩阵大小重新创建矩阵。

S1 = west0479;
[i,j,s] = find(S1);
[m,n] = size(S1);
S2 = sparse(i,j,s,m,n);

稀疏矩阵运算中的索引

由于稀疏矩阵是以压缩稀疏列格式存储的,因此为稀疏矩阵进行索引的相关成本与为满矩阵进行索引的相关成本不同。在只需更改稀疏矩阵中的若干元素时,这类成本可忽略不计,因此,在这类情况下,通常使用常规数组索引来重新分配值:

B = speye(4);
[i,j,s] = find(B);
[i,j,s]
ans =

     1     1     1
     2     2     1
     3     3     1
     4     4     1
B(3,1) = 42;
[i,j,s] = find(B);
[i,j,s]
ans =

     1     1     1
     3     1    42
     2     2     1
     3     3     1
     4     4     1
在存储新矩阵时,为使 42 位于 (3,1) 位置,MATLAB 会在非零值向量和下标向量中插入额外的一行,然后移动 (3,1) 后面的所有矩阵值。

如果线性索引超过 2^48-1(即当前矩阵中允许的元素数上限),使用线性索引在大型稀疏矩阵中访问或指定元素将失败。

S = spalloc(2^30,2^30,2);
S(end) = 1
Maximum variable size allowed by the program is exceeded.

要访问其线性索引大于 intmax 的元素,请使用数组索引:

S(2^30,2^30) = 1
S =

         (1073741824,1073741824)              1

尽管在稀疏矩阵中进行索引以更改单个元素的成本可忽略不计,但该成本在循环环境下会增加,而且在大型矩阵中该操作可能会使执行速度变得很慢。因此,在需要更改大量稀疏矩阵元素的情况下,最好使用向量化方法而不要使用循环方法来执行该操作。例如,考虑稀疏单位矩阵:

n = 10000;
A = 4*speye(n);
以循环方式更改 A 的元素慢于类似的向量化运算:
tic
A(1:n-1,n) = -1; 
A(n,1:n-1) = -1; 
toc
Elapsed time is 0.003344 seconds.
tic
for k = 1:n-1
  C(k,n) = -1; 
  C(n,k) = -1; 
end
toc
Elapsed time is 0.448069 seconds.
由于 MATLAB 以压缩稀疏列格式存储稀疏矩阵,因此,在循环的每个遍历期间,它都需要移动 A 中的多个条目。

如果为稀疏矩阵预分配内存,然后以类似的逐个元素的方式填充,会使对稀疏数组进行索引产生大量开销:

S1 = spalloc(1000,1000,100000);
tic;
for n = 1:100000
    i = ceil(1000*rand(1,1));
    j = ceil(1000*rand(1,1));
    S1(i,j) = rand(1,1);
end
toc
Elapsed time is 2.577527 seconds.

构建索引和值向量则无需为稀疏数组进行索引,因此这种方法的速度快得多:

i = ceil(1000*rand(100000,1));
j = ceil(1000*rand(100000,1));
v = zeros(size(i));
for n = 1:100000
    v(n) = rand(1,1);
end

tic;
S2 = sparse(i,j,v,1000,1000);
toc
Elapsed time is 0.017676 seconds.

因此,最好使用构造函数(例如 sparsespdiags 函数)一次构造所有稀疏矩阵。

例如,假定需要稀疏形式的坐标矩阵 C

C=(4000104001004010101014114)

使用 sparse 函数,以及行下标、列下标和值组成的三联对组,直接构造该五列矩阵:

i = [1 5 2 5 3 5 4 5 1 2 3 4 5]';
j = [1 1 2 2 3 3 4 4 5 5 5 5 5]';
s = [4 1 4 1 4 1 4 1 -1 -1 -1 -1 4]';
C = sparse(i,j,s)
C =

   (1,1)        4
   (5,1)        1
   (2,2)        4
   (5,2)        1
   (3,3)        4
   (5,3)        1
   (4,4)        4
   (5,4)        1
   (1,5)       -1
   (2,5)       -1
   (3,5)       -1
   (4,5)       -1
   (5,5)        4
输出中值的顺序反映了底层的按列存储。有关 MATLAB 如何存储稀疏矩阵的详细信息,请参阅约翰·R·吉尔伯特、克里夫·莫勒尔和罗伯特·施赖伯合著的 Sparse Matrices In MATLAB:Design and Implementation, (SIAM Journal on Matrix Analysis and Applications, 13:1, 333–356 (1992))。

可视化稀疏矩阵

以图的形式查看非零元素在稀疏矩阵内的分布通常很有用。MATLAB spy 函数生成稀疏结构的模板视图,其中图表中的每点代表一个非零数组元素的位置。

例如:

加载提供的稀疏矩阵 west0479,该矩阵是 Harwell-Boeing 集合之一。

load west0479

查看稀疏结构体。

spy(west0479)

Figure contains an axes object. The axes object with xlabel nz = 1887 contains a line object which displays its values using only markers.

另请参阅

相关主题