使用稀疏矩阵时的最佳实践
Ste*_*fin 12 performance benchmarking matlab sparse-matrix
这个问题是基于这个问题的讨论.我之前一直在使用稀疏矩阵,我相信我使用它们的方式很有效.
我的问题是双重的:
在下面,A = full(S)其中S是稀疏矩阵.
访问稀疏矩阵中元素的"正确"方法是什么？
也就是说,稀疏等价物var = A(row, col)是什么？
我对这个话题的看法:你不会做任何不同的事情.var = S(row, col)尽可能高效.我在这方面受到了以下解释的挑战:
按照你所说的方式访问第2行和第2列的元素,S(2,2)与添加新元素相同:
var = S(2,2)=>A = full(S)=>var = A(2,2)=>S = sparse(A) => 4.
这句话真的可以说对吗？
将元素添加到稀疏矩阵的"正确"方法是什么？
也就是说,稀疏等价物会A(row, col) = var是什么？(假设A(row, col) == 0开始)
众所周知,A(row, col) = var对于大型稀疏矩阵来说,简单地做就很慢.从文档:
如果要更改此矩阵中的值,可能会尝试使用相同的索引:
B(3,1)= 42; %此代码确实有效,但速度很慢.
我对这个主题的看法:当使用稀疏矩阵时,你经常从向量开始并使用它们以这种方式创建矩阵:S = sparse(i,j,s,m,n).当然,您也可以像这样创建它:S = sparse(A)或者sprand(m,n,density)类似的东西.
如果从第一种方式开始,您只需执行以下操作:
i = [i; new_i];
j = [j; new_j];
s = [s; new_s];
S = sparse(i,j,s,m,n);
如果你开始没有向量,你会做同样的事情,但find首先使用:
[i, j, s] = find(S);
i = [i; new_i];
j = [j; new_j];
s = [s; new_s];
S = sparse(i,j,s,m,n);
现在你当然会有这些向量,如果你多次进行这个操作,可以重用它们.然而,最好一次添加所有新元素,而不是在循环中执行上述操作,因为增长矢量很慢.在这种情况下,new_i,new_j和new_s将对应于该新的元素向量.
cod*_*ppo 12
MATLAB以压缩列格式存储稀疏矩阵.这意味着当您执行类似A(2,2)的操作(以获取第2行第2列中的元素)时,MATLAB首先访问第二列,然后在第2行中找到元素(每列中的行索引都存储按升序排列).您可以将其视为:
A2 = A(:,2);
A2(2)
如果您只访问稀疏矩阵的单个元素,那么var = S(r,c)就可以了.但是如果你循环遍历稀疏矩阵的元素,你可能希望一次访问一列,然后通过遍历非零行索引[i,~,x]=find(S(:,c)).或者使用类似的东西spfun.
你应该避免构造一个密集的矩阵A然后做S = sparse(A),因为这个操作只是挤出零.相反,正如您所注意到的,使用三元组形式和调用来从头开始构建稀疏矩阵要高效得多sparse(i,j,x,m,n).MATLAB有一个很好的页面,描述了如何有效地构造稀疏矩阵.
描述在MATLAB中实现稀疏矩阵的原始论文是一个很好的阅读.它提供了有关如何最初实现稀疏矩阵算法的更多信息.
编辑:根据Oleg的建议修改答案(见评论).
这是我问题第二部分的基准.为了测试直接插入,矩阵初始化为空变化nzmax.为了测试从索引向量重建,这是无关紧要的,因为矩阵是在每次调用时从头开始构建的.测试这两种方法进行单次插入操作(具有不同数量的元素),或者进行增量插入,一次一个值(直到相同数量的元素).由于计算应变,我将每个测试用例的重复次数从1000降低到100.我相信这在统计上仍然可行.
Ssize = 10000;
NumIterations = 100;
NumInsertions = round(logspace(0, 4, 10));
NumInitialNZ = round(logspace(1, 4, 4));
NumTests = numel(NumInsertions) * numel(NumInitialNZ);
TimeDirect = zeros(numel(NumInsertions), numel(NumInitialNZ));
TimeIndices = zeros(numel(NumInsertions), 1);
%% Single insertion operation (non-incremental)
% Method A: Direct insertion
for iInitialNZ = 1:numel(NumInitialNZ)
disp(['Running with initial nzmax = ' num2str(NumInitialNZ(iInitialNZ))]);
for iInsertions = 1:numel(NumInsertions)
tSum = 0;
for jj = 1:NumIterations
S = spalloc(Ssize, Ssize, NumInitialNZ(iInitialNZ));
r = randi(Ssize, NumInsertions(iInsertions), 1);
c = randi(Ssize, NumInsertions(iInsertions), 1);
tic
S(r,c) = 1;
tSum = tSum + toc;
end
disp([num2str(NumInsertions(iInsertions)) ' direct insertions: ' num2str(tSum) ' seconds']);
TimeDirect(iInsertions, iInitialNZ) = tSum;
end
end
% Method B: Rebuilding from index vectors
for iInsertions = 1:numel(NumInsertions)
tSum = 0;
for jj = 1:NumIterations
i = []; j = []; s = [];
r = randi(Ssize, NumInsertions(iInsertions), 1);
c = randi(Ssize, NumInsertions(iInsertions), 1);
s_ones = ones(NumInsertions(iInsertions), 1);
tic
i_new = [i; r];
j_new = [j; c];
s_new = [s; s_ones];
S = sparse(i_new, j_new ,s_new , Ssize, Ssize);
tSum = tSum + toc;
end
disp([num2str(NumInsertions(iInsertions)) ' indexed insertions: ' num2str(tSum) ' seconds']);
TimeIndices(iInsertions) = tSum;
end
SingleOperation.TimeDirect = TimeDirect;
SingleOperation.TimeIndices = TimeIndices;
%% Incremental insertion
for iInitialNZ = 1:numel(NumInitialNZ)
disp(['Running with initial nzmax = ' num2str(NumInitialNZ(iInitialNZ))]);
% Method A: Direct insertion
for iInsertions = 1:numel(NumInsertions)
tSum = 0;
for jj = 1:NumIterations
S = spalloc(Ssize, Ssize, NumInitialNZ(iInitialNZ));
r = randi(Ssize, NumInsertions(iInsertions), 1);
c = randi(Ssize, NumInsertions(iInsertions), 1);
tic
for ii = 1:NumInsertions(iInsertions)
S(r(ii),c(ii)) = 1;
end
tSum = tSum + toc;
end
disp([num2str(NumInsertions(iInsertions)) ' direct insertions: ' num2str(tSum) ' seconds']);
TimeDirect(iInsertions, iInitialNZ) = tSum;
end
end
% Method B: Rebuilding from index vectors
for iInsertions = 1:numel(NumInsertions)
tSum = 0;
for jj = 1:NumIterations
i = []; j = []; s = [];
r = randi(Ssize, NumInsertions(iInsertions), 1);
c = randi(Ssize, NumInsertions(iInsertions), 1);
tic
for ii = 1:NumInsertions(iInsertions)
i = [i; r(ii)];
j = [j; c(ii)];
s = [s; 1];
S = sparse(i, j ,s , Ssize, Ssize);
end
tSum = tSum + toc;
end
disp([num2str(NumInsertions(iInsertions)) ' indexed insertions: ' num2str(tSum) ' seconds']);
TimeIndices(iInsertions) = tSum;
end
IncremenalInsertion.TimeDirect = TimeDirect;
IncremenalInsertion.TimeIndices = TimeIndices;
%% Plot results
% Single insertion
figure;
loglog(NumInsertions, SingleOperation.TimeIndices);
cellLegend = {'Using index vectors'};
hold all;
for iInitialNZ = 1:numel(NumInitialNZ)
loglog(NumInsertions, SingleOperation.TimeDirect(:, iInitialNZ));
cellLegend = [cellLegend; {['Direct insertion, initial nzmax = ' num2str(NumInitialNZ(iInitialNZ))]}];
end
hold off;
title('Benchmark for single insertion operation');
xlabel('Number of insertions'); ylabel('Runtime for 100 operations [sec]');
legend(cellLegend, 'Location', 'NorthWest');
grid on;
% Incremental insertions
figure;
loglog(NumInsertions, IncremenalInsertion.TimeIndices);
cellLegend = {'Using index vectors'};
hold all;
for iInitialNZ = 1:numel(NumInitialNZ)
loglog(NumInsertions, IncremenalInsertion.TimeDirect(:, iInitialNZ));
cellLegend = [cellLegend; {['Direct insertion, initial nzmax = ' num2str(NumInitialNZ(iInitialNZ))]}];
end
hold off;
title('Benchmark for incremental insertions');
xlabel('Number of insertions'); ylabel('Runtime for 100 operations [sec]');
legend(cellLegend, 'Location', 'NorthWest');
grid on;
我在MATLAB R2012a中运行了这个.执行单个插入操作的结果总结在此图中:
这表明如果只进行一次操作,使用直接插入比使用索引向量要慢得多.使用索引向量的情况下的增长可能是因为向量本身的增长或者是由较长的稀疏矩阵构造,我不确定哪个.nzmax用于构建矩阵的初始化似乎对它们的增长没有影响.
此图表总结了进行增量插入的结果:
在这里我们看到了相反的趋势:使用索引向量较慢,因为逐步增长它们并在每一步重建稀疏矩阵的开销.理解这一点的一种方法是查看上一个图中的第一点:对于插入单个元素,使用直接插入而不是使用索引向量进行重建更有效.在incrementlal的情况下,这个单个插入是重复完成的,因此根据MATLAB的建议使用直接插入而不是索引向量变得可行.
这种理解还表明,如果我们一次逐步添加100个元素,那么有效的选择就是使用索引向量而不是直接插入,因为第一个图表显示此方法对于此大小的插入更快.在这两个体系之间是一个你可能应该试验看哪种方法更有效的领域,尽管结果可能表明这些方法之间的区别是可忽略不计的.
底线:我应该使用哪种方法？
我的结论是,这取决于您的预期插入操作的性质.
- 如果您打算一次插入一个元素,请使用直接插入.
- 如果要一次插入大(> 10)个元素,请从索引向量重建矩阵.