稀疏矩阵运算中的索引
由于稀疏矩阵是以压缩稀疏列格式存储的,因此为稀疏矩阵进行索引的相关成本与为满矩阵进行索引的相关成本不同。在只需更改稀疏矩阵中的若干元素时,这类成本可忽略不计,因此,在这类情况下,通常使用常规数组索引来重新分配值:
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.
因此,最好使用构造函数(例如 sparse 或 spdiags 函数)一次构造所有稀疏矩阵。
例如,假定需要稀疏形式的坐标矩阵 C:
C=(4000−10400−10040−101010141−14)
使用 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 如何存储稀疏矩阵的详细信息,请参阅 John R. Gilbert、Cleve Moler 和 Robert Schreiber 合著的 Sparse Matrices In Matlab:Design and Implementation, (SIAM Journal on Matrix Analysis and Applications, 13:1, 333–356 (1992))。