Efficient multiplication by large structured matrix of ones and zeros

I am trying to do a relatively simple operation, with a prototypical example:

N = 20;
B = de2bi((1:2^N)-1,N);

c = rand(2^N,1);

result = sum(B.*c,1);

Now this works well, and seems relatively efficient, but when N grows larger, the B matrix will be 2^N x N, which gets prohibitive quite fast. I could make it a sparse matrix which would save half the memory, but that just buys one extra N before memory constraints are reached again. I'm wondering if there are any clever and efficient ways to perform this operation without necessarily saving the B matrix given that it is so structured. Thanks!

Answer 1

[Matt J]

It may seem counter-intuitive that this version is faster than the versions in my earlier comment, since it does twice as much summation. I imagine it's because it's better multi-threaded and requires less memory allocation.

N = 24; % not so small
c = rand(2^N,1);

tic;
    chat = reshape(c,repmat(2,1,N));
    e=1:N; 
    
    clear result2
    for i = 1:N
        tmp=sum(chat,setdiff(e,i));
        result2(N+1-i) = tmp(2); %#ok<SAGROW> 
    end
 toc

Elapsed time is 0.260331 seconds.

[John D'Errico]

Nice. The funny thing is, when I first looked at this question, I was thinking that it would not be possible to see a gain over the brute force matrix multiply. My gut was wrong of course.

[Adam Shaw]

Fantastic solution, thanks! Note that I think you just need result2(i) = tmp(2) - it seems flipped from the brute force matrix multiplication to me.

Answer 2

[John D'Errico]

Are there more efficient ways to do this? Well, perhaps, arguably so. Why you want to do it, is I supose your problem.

THINK about what you are doing here. What does that multiply do? It forms linear combinations of the vector C, with the bits of the numbers 0:2^N-1.

For example...

N = 5; % small, so we can see what happens.
B = de2bi((1:2^N)-1,N);
c = rand(2^N,1);
format long g
result = sum(B.*c,1)

result = 1×5

9.65859767786247 7.1472229443395 9.89441940297802 8.23946581063777 8.39175263434243

Now using a different code, can we reproduce that? One idea is to convert c into a multi-dimensional array. Then sum only the appropriate elements, permuting the elements of that array.

chat = reshape(c,repmat(2,1,N));
result2 = zeros(1,N);
for i = 1:N
  chati = reshape(permute(chat,[i,setdiff(1:N,i)]),2,[]);
  result2(i) = sum(chati(2,:));
end
result2

result2 = 1×5

9.65859767786247 7.1472229443395 9.89441940297802 8.23946581063777 8.39175263434243

So I never computed B. I did a fair amount of work to perform the computations though. Is there a gain when N is large? It may start to gain when N grows moderately large.

[Matt J]

No need to permute:

N = 24; % not so small
c = rand(2^N,1);

tic
  B = de2bi((1:2^N)-1,N);
  result = sum(B.*c,1);
toc

Elapsed time is 8.992176 seconds.

tic;
    chat = reshape(c,repmat(2,1,N));
    result2 = zeros(1,N);
    for i = 1:N
      chati = reshape(permute(chat,[i,setdiff(1:N,i)]),2,[]);
      result2(i) = sum(chati(2,:));
    end
toc

Elapsed time is 2.891254 seconds.

tic;

    chat = reshape(c,repmat(2,1,N));
    inds0=repmat({':'},1,N);
    
    result2 = zeros(1,N);
    for i = 1:N
      inds=inds0; 
      inds{i}=2;  
      result2(i) = sum(chat(inds{:}),'all');
    end
toc

Elapsed time is 1.551022 seconds.