MATLAB: How to record all cases where two columns sum to a third

vectorization

I have a matrix and I would like to record all instances of triplets of columns where two columns sum to the third. So let us say the matrix is
2 0 1 2 1
1 1 0 0 1
Then I would like to produce a new matrix like so:
1 1 5
2 3 2
4 5 3
What this is saying is that column 1 can be expressed as the sum of columns 2 and 4, or as the sum of columns 3 and 5. In addition, column 5 can be expressed as the sum of columns 2 and 3.
What I have tried so far is to use loops:
D=[2,0,1,2,1;1,1,0,0,1];
N=size(D,2);
eq=[];
for ii=1:N
  if ii==1
    for jj=(ii+1):(N-1)
      for kk=(jj+1):N
        if isequal(D(:,ii),D(:,jj)+D(:,kk))
          eq=[eq,[ii;jj;kk]];
        end
      end
    end
  elseif ii==2
    for jj=[ii-1,(ii+1):(N-1)]
      for kk=(max(ii,jj)+1):N
        if isequal(D(:,ii),D(:,jj)+D(:,kk))
          eq=[eq,[ii;jj;kk]];
        end
      end
    end
  elseif ii>2 && ii<(N-1)
    for jj=[1:(ii-1),(ii+1):(N-1)]
      for kk=[(jj+1):(ii-1),(max(ii,jj)+1):N]
        if isequal(D(:,ii),D(:,jj)+D(:,kk))
          eq=[eq,[ii;jj;kk]];
        end
      end
    end
  elseif ii==(N-1)
    for jj=1:(ii-1)
      for kk=[(jj+1):(ii-1),N]
        if isequal(D(:,ii),D(:,jj)+D(:,kk))
          eq=[eq,[ii;jj;kk]];
        end
      end
    end
  elseif ii==N
    for jj=1:(ii-2)
      for kk=jj:(ii-1)
        if isequal(D(:,ii),D(:,jj)+D(:,kk))
          eq=[eq,[ii;jj;kk]];
        end
      end
    end
  end          
end
The code works, but it is convoluted and also very slow if my matrix D has lots of columns. I have a hunch that this could be vectorized, but I cannot figure out how. Does anyone know of a more elegant solution?

Best Answer

A solution that does not involve any loop:
m = [2 0 1 2 1
     1 1 0 0 1];
[col1, col2] = ndgrid(1:size(m, 2));   %cartesian product of column indices
mask = triu(true(size(col1)), 1);      %mask to only keep columns where col2 > col1
col1 = col1(mask);
col2 = col2(mask);
[issum, colsum] = ismember((m(:, col1) + m(:, col2)).', m.', 'rows');  %note the transpose because of the 'rows' option
result = [colsum(issum), col1(issum), col2(issum)].'  %use sortrows before the transpose if the order is important
Related Question