14

My question is twofold:

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

What's the "correct" way to access an element in a sparse matrix?

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

My view on this topic: You wouldn't do anything different. var = S(row, col) is as efficient as it gets.

What's the "correct" way to add elements to a sparse matrix?

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

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; % This code does work, however, it is slow.

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.

Stewie Griffin
  • 14,889
  • 11
  • 39
  • 70
  • 3
    In reference to accessing elements of a sparse matrix, [there is no conversion to full](http://stackoverflow.com/a/23023268/2180721). – Oleg Apr 11 '14 at 22:01

2 Answers2

13

MATLAB stores sparse matrices in compressed column format. This means that when you perform an operations like A(2,2) (to get the element in at row 2, column 2) MATLAB first access the second column and then finds the element in row 2 (row indices in each column are stored in ascending order). You can think of it as:

 A2 = A(:,2);
 A2(2)

If you are only accessing a single element of sparse matrix doing var = S(r,c) is fine. But if you are looping over the elements of a sparse matrix, you probably want to access one column at a time, and then loop over the nonzero row indices via [i,~,x]=find(S(:,c)). Or use something like spfun.

You should avoid constructing a dense matrix A and then doing S = sparse(A), as this operations just squeezes out zeros. Instead, as you note, it's much more efficient to build a sparse matrix from scratch using triplet-form and a call to sparse(i,j,x,m,n). MATLAB has a nice page which describes how to efficiently construct sparse matrices.

The original paper describing the implementation of sparse matrices in MATLAB is quite a good read. It provides some more info on how the sparse matrix algorithms were originally implemented.

codehippo
  • 1,379
  • 1
  • 8
  • 19
  • Do you mean that `A2` is actually placed in memory (but not shown)? And if so, is A2 a sparse vector or a regular vector? – Dennis Jaheruddin Apr 14 '14 at 11:12
  • 1
    @DennisJaheruddin this is just a way to think about how the algorithm for selecting the A(2,2) would work. A2 is not created, it just represents the (sparse) second column of the sparse matrix A. – codehippo Apr 14 '14 at 13:22
  • @DennisJaheruddin If you actually read the paper on page 9, the column is created as "3.1.3. The sparse accumulator", and is dense. – Oleg Apr 14 '14 at 20:12
  • I'll give you the bounty, as I think the last paper covers it all. :-) thanks! I hope anyone interested cares to read the paper :-) however, I'll keep the question open for a while, in case others have any additional valuable input... – Stewie Griffin Apr 18 '14 at 14:28
6

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

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;

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

Sparse matrix insertion benchmark: single insertion operation

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:

Sparse matrix insertion benchmark: incremental insertions

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.

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.

Bottom line: which method should I use?

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

  • If you intend to insert elements one at a time, use direct insertion.
  • If you intend to insert a large (>10) number of elements at a time, rebuild the matrix from index vectors.
buzjwa
  • 2,632
  • 2
  • 24
  • 37
  • Your test leaves out two important aspects (that IMHO change the verdict): 1) pre-allocation of nzmax elements and 2) concatenation/query of indices,e.g. `i_new`. With the first point I mean that you can reserve some additional space for new values to avoid any copy-write of indices and values. My guess is that a bit of this pre-allocation effect can be seen from your test since you allocate approximately 1e4^2*0.0005 ~ 5e5 values and then zero-them out. However, the copy-write penalty kicks in heavily only after 1e4*(100 insertions). So, an explicit test is needed. – Oleg Apr 13 '14 at 22:11
  • With point two, I mean that a real use would require adding values incrementally, i.e. I would first query the new indices (or concatenate growing the old ones), then add some more and re-construct the matrix, and repeat the same on the next iteration. In general, I think smarter use of sparse matrices requires getting an estimate or a guess for the max number of values you might need to allocate, and then simply use direct assignment. – Oleg Apr 13 '14 at 22:11
  • You are making excellent points here. I expanded the tests and will update my answer once they're done. I agree with your last statement and hope to verify it by this benchmark. – buzjwa Apr 14 '14 at 07:14
  • @OlegKomarov See my updated answer. It turns out the initial `nzmax` didn't change much. Let me know what you think. – buzjwa Apr 14 '14 at 11:34
  • @Naveh I think you reached the wrong conclusion. You assume you need to do one insertion at a time (and rebuild the sparse matrix each time you do an insertion)? Why not, as the OP said, create i,j,s vectors that contain all of your insertions and then only create the sparse matrix once? If you do that, using index vectors will be much faster than direct insertions for your second test. – codehippo Apr 14 '14 at 13:35
  • I agree that building only once from index vectors is fastest. Incremental insertion is not an assumption, just one of the options I checked. If you only need to insert elements once then use index vectors (this is the conclusion from the first graph). If however you need to insert elements one at a time, use direct insertion (this is the conclusion from the second graph). – buzjwa Apr 14 '14 at 13:44
  • 1
    @codehippo, I agree with Naveh on this one. In my opinion, the answer is unambiguous. If you are inserting a chunk of elements, use indices. If you want to insert one element at a time, direct insertion is faster. Let's say your initial matrix is `S0`. Assuming you want to insert `k` elements, your final matrix is `Sk`. If you're only interested in `S0` and `Sk`, use indices. However, if you for some reason is interested in `S1`, `S2`, `S3` ..., use direct insertion. – Stewie Griffin Apr 14 '14 at 15:08
  • @RobertP. Consider the case where you preform a direct insertion near the beginning of the first column. My understanding is that this will force nearly all the nonzeros to be moved (since the matrix is stored in CSC). My guess is that this will cost the same as forming the matrix from scratch, even though you are only inserting a single element. – codehippo Apr 14 '14 at 15:15
  • With respect to my point 1), you are actually inserting already more than what you allocated. However, I wouldn't bother going too much in details. Creating a sparse matrix by columns, is the way to go for multiple insertions. Although intuition would suggest that direct insertion of a few elements should get slower as the insert gets closer to the first elements, meaning more elements afterwards will need to be re-written, and thus it would also depend on the size of the array being re-written, I did not manage to overthrow the result by re-building the sparse array. – Oleg Apr 14 '14 at 20:16
  • @codehippo My previous comment addresses your last point. In general rebuilding the sparse array by concatenating subscripts with the new element will be a worse alternative to direct insert, I guess due to the re-write of the indices. – Oleg Apr 14 '14 at 20:23