处理稀疏矩阵时的最佳实践 [英] Best practice when working with sparse matrices

问题描述

此问题基于此问题中的讨论.我之前一直在处理稀疏矩阵，我的信念是我一直在使用稀疏矩阵.

This question is based on a discussion in this question. I have been working with sparse matrices earlier and my belief is that the way I've been working with these is efficient.

我的问题有两个:

在下面，A = full(S)，其中S是稀疏矩阵.

In the below, A = full(S) where S is a sparse matrix.

也就是说，与var = A(row, col)等效的稀疏是什么?

That is, what would the sparse equivalent to var = A(row, col) be?

我对这个主题的看法:您将不会做任何不同的事情. var = S(row, col)一样高效.对此我的解释受到了挑战:

My view on this topic: You wouldn't do anything different. var = S(row, col) is as efficient as it gets. I have been challenged on this with the explanation:

以您所说的S(2,2)的方式访问第2行和第2列上的元素 与添加新元素相同:var = S(2,2) => A = full(S) => var = A(2,2) => S = sparse(A) => 4.

Accessing element on row 2 and column 2 the way you said, S(2,2) does the same as adding a new element: var = S(2,2) => A = full(S) => var = A(2,2) => S = sparse(A) => 4.

这个说法真的正确吗?

也就是说，A(row, col) = var的稀疏等效项是什么? (假设A(row, col) == 0开始)

That is, what would the sparse equivalent of A(row, col) = var be? (Assuming A(row, col) == 0 to begin with)

众所周知，对于大型稀疏矩阵，简单地执行A(row, col) = var会很慢.从文档:

It is known that simply doing A(row, col) = var is slow for large sparse matrices. From the documentation:

如果您想更改此矩阵中的值，可能会很受诱惑 使用相同的索引:

If you wanted to change a value in this matrix, you might be tempted to use the same indexing:

B(3,1)= 42; ％该代码确实有效，但是速度很慢.

B(3,1) = 42; % This code does work, however, it is slow.

我在此主题上的观点:在处理稀疏矩阵时，通常从向量开始，并使用它们以这种方式创建矩阵:S = sparse(i,j,s,m,n).当然，您也可以这样创建它:S = sparse(A)或sprand(m,n,density)或类似的东西.

My view on this topic: When working with sparse matrices, you often start with the vectors and use them to create the matrix this way: S = sparse(i,j,s,m,n). Of course, you could also have created it like this: S = sparse(A) or sprand(m,n,density) or something similar.

如果您以第一种方式开始，只需执行以下操作即可:

If you start of the first way, you would simply do:

i = [i; new_i];
j = [j; new_j];
s = [s; new_s];
S = sparse(i,j,s,m,n); 

如果一开始就没有向量，您将做同样的事情，但是要使用

If you started out not having the vectors, you would do the same thing, but use find first:

[i, j, s] = find(S);
i = [i; new_i];
j = [j; new_j];
s = [s; new_s];
S = sparse(i,j,s,m,n); 

现在，您当然将拥有向量，并且如果您多次执行此操作，则可以重用它们.但是，最好一次添加所有新元素，而不是循环执行上述操作，因为

Now you would of course have the vectors, and can reuse them if you're doing this operation several times. It would however be better to add all new elements at once, and not do the above in a loop, because growing vectors are slow. In this case, new_i, new_j and new_s will be vectors corresponding to the new elements.

推荐答案

编辑:答案根据Oleg的建议进行了修改(请参见评论).

EDIT: Answer modified according to suggestions by Oleg (see comments).

这是我对问题第二部分的基准.为了测试直接插入，使用变化的nzmax将矩阵初始化为空.对于测试从索引向量进行的重建，这是不相关的，因为矩阵是在每次调用时从头开始构建的.对这两种方法进行了测试(一次插入操作(元素数量不同))或一次递增插入(一次插入一个值)(最多元素数量相同).由于计算量大，我将每个测试用例的重复次数从1000减少到100.我相信这在统计上还是可行的.

Here is my benchmark for the second part of your question. For testing direct insertion, the matrices are initialized empty with a varying nzmax. For testing rebuilding from index vectors this is irrelevant as the matrix is built from scratch at every call. The two methods were tested for doing a single insertion operation (of a varying number of elements), or for doing incremental insertions, one value at a time (up to the same numbers of elements). Due to the computational strain I lowered the number of repetitions from 1000 to 100 for each test case. I believe this is still statistically viable.

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中运行了它.下图总结了执行一次插入操作的结果:

I ran this in MATLAB R2012a. The results for doing a single insertion operations are summarized in this graph:

这表明，如果仅执行一次操作，则使用直接插入比使用索引向量要慢得多.使用索引向量的情况下的增长可能是由于向量本身的增长，也可能是由于更长的稀疏矩阵构造引起的，我不确定是哪个.最初用于构建矩阵的nzmax似乎对其生长没有影响.

This shows that using direct insertion is much slower than using index vectors, if only a single operation is done. The growth in the case of using index vectors can be either because of growing the vectors themselves or from the lengthier sparse matrix construction, I'm not sure which. The initial nzmax used to construct the matrices seems to have no effect on their growth.

此图总结了进行增量插入的结果:

The results for doing incremental insertions are summarized in this graph:

在这里，我们看到了相反的趋势:由于逐步增加索引向量并在每个步骤重建稀疏矩阵的开销，使用索引向量的速度较慢.理解这一点的一种方法是查看上图中的第一点:对于插入单个元素，使用直接插入比使用索引向量进行重建更有效.在增量情况下，这种单次插入是重复进行的，因此违反MATLAB的建议，使用直接插入而不是索引向量变得可行.

Here we see the opposite trend: using index vectors is slower, because of the overhead of incrementally growing them and rebuilding the sparse matrix at every step. A way to understand this is to look at the first point in the previous graph: for insertion of a single element, it is more effective to use direct insertion rather than rebuilding using the index vectors. In the incrementlal case, this single insertion is done repetitively, and so it becomes viable to use direct insertion rather than index vectors, against MATLAB's suggestion.

这种理解还表明，如果我们一次递增地添加100个元素，那么有效的选择将是使用索引向量而不是直接插入，因为第一张图显示了该方法对于插入的速度更快.这个大小.在这两种方案之间的一个区域中，您可能应该进行实验以查看哪种方法更有效，尽管结果可能表明那里的方法之间的差异可忽略不计.

This understanding also suggests that were we to incrementally add, say, 100 elements at a time, the efficient choice would then be to use index vectors rather than direct insertion, as the first graph shows this method to be faster for insertions of this size. In between these two regimes is an area where you should probably experiment to see which method is more effective, though probably the results will show that the difference between the methods is neglibile there.

底线:我应该使用哪种方法?

我的结论是，这取决于您要进行的插入操作的性质.

My conclusion is that this is dependant on the nature of your intended insertion operations.

• 如果您打算一次插入一个元素，请使用直接插入.
• 如果您打算一次插入大量(> 10)的元素，请从索引向量重建矩阵.

这篇关于处理稀疏矩阵时的最佳实践的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！