explain sparse convolution algorithm

There are (at least) two functions to convolve two sparse arrays on the file exchange: sconv2 and the convn method of n-D sparse array class. Both use a similar algorithm, which I have re-written for 1D data:

    % 
function C = SparseConv1(A, B)
assert( isvector(A) );
A = A(:);
assert( isvector(B) );
B = B(:);
[idxA,~,avals] = find(A(:));
[idxB,~,bvals] = find(B(:));
C = avals(:)*bvals(:).';
idxArep = repmat(idxA(:),1,numel(idxB));
idxBrep = repmat(idxB(:).',numel(idxA),1); 
        % sparse adds together elements that have duplicate subscripts
        % but don't understand how this gives the convolution.
C = sparse(idxArep(:) + idxBrep(:) -1, ones(numel(idxArep),1), C(:));

Question: Can someone derive the algorithm from the traditional sum of product of offset vectors formula?

Or provide a reference to the derivation?

edit: 6/11/17 I corrected the expressions for idxArep and idxBrep

edit: 6/23/17 I added another simplified version as an answer that shows the algorithm more clearly

[David Goodmanson]

Hello Robert, I was taking a look at this and it does not appear to produce the same results as conv:

A = [1 2 3];   
B = [3 5 7 9];
x = SparseConv1(A,B)
y = conv(A,B)

x =
   (1,1)        3
   (2,1)       16
   (3,1)       45
   (4,1)       21
   (5,1)       32
   (6,1)       27

y =
     3    11    26    38    39    27

[Bruno Luong]

OP made a mistake in his transcription, it should be:

idxBrep = repmat(idxB(:).',numel(idxA),1);
idxArep = repmat(idxA(:),1,numel(idxB));

[Robert Alvarez]

I corrected the code.

[Robert Alvarez]

I think one way to derive the algorithm is to observe that we can expand the A vector as a sum of basis vectors {u_k} of length numel(A) with a '1' at location k and zeros elsewhere multiplied by the scalars A(k).

A = Sum[k=1..numel(A)] A(k)*u_k

With a similar expansion for B.

The u_k vectors are orthonormal.

Answer 1

[David Goodmanson]

Hi Robert,

Since the posted code is oriented toward sparse matrices, it uses ‘find’ to eliminate the calculation for zeros contained A and B, of which there are probably a lot. The code below parallels the posted code but for simplicity it takes all possible values of A and B and does not bother with find. This means that idxA = 1:length(A), similarly for B. It includes the correction, thanks to Bruno Luong.

Convolution has the form

D(j) = Sum{i}  A(i)*B(j+1-i).

This is equivalent to

D(j) = Sum{i,k}  A(i)*B(k)   ---  with the condition k = j+1-i  -->  j = i+k-1.

(The 1 is there so that when 1=1 and k=1, then j=1 and not 2).

The algorithm needs every possible product of an element of A with an element of B, and the outer product M = A*B.’ is a matrix of this form, where

M(i,k) = A(i)*B(k).

Each matrix index pair (i,k) corresponds to j = i+k-1, and the values of j look like

[1 2 3 
 2 3   ...
 3

You just need to construct a matrix like that in order to grab all the right entries of M for each j, and sum over those entries. At this point it’s probably easiest to do an example.

>> A = [1 2 3].';      B = [3 5 7 9].';
>> M = A*B.'

M =
     3     5     7     9
     6    10    14    18
     9    15    21    27

>> nA = length(A); nB = length(B);
>> indA = repmat((1:nA)',1,nB)
>> indB = repmat(1:nB,nA,1)
>> Z = indA+indB-1

indA =
     1     1     1     1
     2     2     2     2
     3     3     3     3

indB =
     1     2     3     4
     1     2     3     4
     1     2     3     4

Z =
     1     2     3     4
     2     3     4     5
     3     4     5     6

Summing entries in M by eye, for constant j in Z, gives, 3, 11, 26 etc. compared to

>> conv(A,B).'
   ans =   3    11    26    38    39    27

The last command in the original code reads out the matrices (in this notation) M, indA and indB column by column (so everything still matches up), creates indices (j,1) for a sparse column vector and does the sum as you mentioned, where sparse matrix elements with the same indices are added.

The two repmat lines can be replaced more compactly by

[indA indB] = ndgrid(1:nA,1:nB)

[Robert Alvarez]

Thanks for the response, David. I am looking for a more formal derivation, particularly for the part where you go to example. From numerical examples, it looks like the algorithm works but it would be useful (and interesting) to have a mathematical proof.

[Robert Alvarez]

I think the explanation needs to be able to show, for example, why the previous formulas for idxArep and idxBrep did not work when the input vectors are of unequal length.

[David Goodmanson]

By the previous formulas do you mean the ones that got corrected? If so, that starts to get into mistakes in coding (those happen) at least as much as proving what all the indices are doing.

If you want a strict mathematical proof I don't think that the sparse matrix way of doing it will get there anyway. That's because after you make the big column vector of index values (I called them j) for the final result,

Vecj =  idxArep(:) + idxBrep(:) -1,

the calculation depends on the sparse algorithm finding a value j, recognizing that it is the same as a previous value j and adding the results into C(j). The method does not tell you where in Vecj the values of j occurred. It is not that hard to work out other indices that tell you those locations in Vecj where a given j occurs, but in that case you would not be using the sparse command any more.

[Robert Alvarez]

"By the previous formulas do you mean the ones that got corrected? If so, that starts to get into mistakes in coding (those happen) at least as much as proving what all the indices are doing."

I think that finding the reason that the current versions of idxArep and idxBrep work and the previous ones did not gives a lot of insight into the algorithm.

Yes, we can see from numerical examples like yours that this is picking off the correct values of C but why/how? An algebraic formula derived from the fundamental definition of a convolution is required to truly understand the algorithm IMO.

[Robert Alvarez]

" It is not that hard to work out other indices that tell you those locations in Vecj where a given j occurs, but in that case you would not be using the sparse command any more."

I think that a useful explanation has to model the behavior of the sparse function regarding adding duplicate entries. And doing this is crucial to understanding the algorithm.

[David Goodmanson]

I totally agree the comparing previous vs. current versions of idxArep & idxBrep provides insight.

I made your code as it is is now into a script. Since the example A and B do not contain zeros, then idxA is simply 1:length(A), same for B. The key is that all three matrices C, idxA and idzB are going to be read out columnwise so their shapes have to match in order for everything to match up in the steps after that. If you run the code you can see the problem. Small examples are important. Once you resolve that issue, we could discuss the sparse duplicate issue if you would like.

A = [1 2 3]
B = [3 5 7 9]
A = A(:);
B = B(:);
[idxA,~,avals] = find(A(:))
[idxB,~,bvals] = find(B(:))
C = avals(:)*bvals(:).'

idxArep = repmat(idxA(:),1,numel(idxB))
idxBrep = repmat(idxB(:).',numel(idxA),1)
idxArep_old = repmat(idxA(:).',numel(idxB),1)
idxBrep_old = repmat(idxB(:),1,numel(idxA))

[Robert Alvarez]

I found a bug with SparseConv1. The following test case gives a result that differs from conv: x = [0 1 2 3 4]; h = [5 6 0]; SparseConv1(x,h) outputs [0 5 16 27 38 24]. conv(x,h) outputs [0 5 16 27 38 24 0]. Note the zero at the end. I need to understand the algorithm better to fix this bug.

[David Goodmanson]

good catch. This is not exactly a bug because you end up with a sparse matrix, all of whose nonzero elements are correct. If you do some test cases with both leading and trailing zeros in x and h, then convert from sparse to full, you will see that the leading zeros come out correctly, the nonzero entries occur in the right locations, and only the trailing zeros are missing.

When 'full' converts from sparse it is able to put in the leading zeros before the first nonzero element because the answer has to begin with element 1. But 'full' has no way of knowing how long the output of conv is supposed to be, which is nA+nB-1. If you put that info into sparse,

C = sparse(idxArep(:)+idxBrep(:)-1, ones(numel(idxArep),1),C(:),nA+nB-1,1)

then full is able to come up with the right number of trailing zeros.

I think the zeros are a distraction, and the best way to approach the sparse conv algorithm is to work it out from beginning to end without treating nonzero entries in h and x as special. In the not-special case you don't need to use 'find', idxA is simply 1:nA, and avals = A; similarly for B. Then come back later and see what happens with special treatment for zeros in A and B.

Sparse is being used in two ways, as a method to save space by not saving zeros explicitly, and for its feature of adding duplicate entries. The second one is the important one. As I mentioned in a previous comment, finding a way to do that using extra indices and without resorting to sparse is how to get to a real proof, and it's a good topic for discussion.

Answer 2

[Robert Alvarez]

A simplified version that shows the algorithm more clearly

    % Simplified code version 2

    function C = Conv1x(A,B)
    idxA = find(A);
    idxB = find(B);

    idxArep = repmat(idxA(:),1,numel(idxB));
    idxBrep = repmat(idxB(:).',numel(idxA),1); 

% indexes for output convolution      
  idxs4conv =  idxArep(:) + idxBrep(:) -1; 

       % compute convolution from indexes and values
    C = zeros(numel(A) + numel(B) - 1, 1);
    idxs4conv_uniq = unique(idxs4conv);
    for k = 1:numel(idxs4conv_uniq);
       idx_uniq = idxs4conv_uniq(k);
       idxs_repeated = find(idxs4conv==idx_uniq);
       for ki = 1:numel(idxs_repeated)
           kA = idxArep(idxs_repeated(ki));
           kB = idxBrep(idxs_repeated(ki));
           C(idx_uniq) = C(idx_uniq) + A(kA)*B(kB);
       end
    end

[David Goodmanson]

Looks like you have scoped this out in its entirety. Good one! There's a lot to be said for working out the details, compared to most people who just plug into an algorithm like this one and are not motivated to do so.

Earlier in this thread I maintained that using 'find' and the like does not lead to a mathematically rigorous process. But having considered further I think that there is nothing wrong with that method.