Fast multiplication: binary matrix with double matrix

Hello everyone,

I am trying to speed up my Matlab code at the moment. The most time consuming part of the code is the multiplication of two matrices A*B, where A is binary (only 0 or 1 entries) and B is a double matrix.

The size of the matrices isn't that large, it's only time consuming because its in the inner loop of some iteration and thus is performed 100k upto a million times. (The matrix B is changing in each iteration, but A stays the same.) So each bit of performance speed-up could help me out here.

At the moment, I am just converting A from binary to double and use the standard matrix multiplication A*B. However, I wonder if there is a faster way to do it. since A is binary, A*B is not a 'real multiplication' but just an addition of some elements of B (defined by the non-zero pattern of A). Anyone has a clue how to do so? And neither A nor B are sparse, if that is important.

[Jos (10584)]

matrix muliplication is highly optimized. I doubt that, for instance, summing selected elements of B would be faster. What are the sizes of A and B? And do you need to create B as a separate matrix, or can you use A to create B?

[Florian]

A and B are normally around 512 x 512. I need to create both as a separate matrix. A is created once in the beginning and B is overwritten in each iteration.

I wonder if maybe using mex and writing the addition in C would be faster. However, I would take me some time to get into all of this. I wanted to ask here first.

[Jos (10584)]

So, just to check, both A and B are square and of the same size? A is a logical matrix, and B is a double matrix created (overwritten) by another function inside a loop?

Writing a mex-function to do this would require the same structure as matrix multiplication (looping and summing) so I doubt you can get it any faster.

What about decreasing the size of A by removing replicated rows?

[Florian]

A and B don't need to be square, but size(A,2)==size(B,1). And yes, the loop is like

for k=1:numIterations
  X = A*B;
  %some other calculations on X
  B = B+X;
end

A does not have replicated rows.

The structure in the mex file would be the same, but instead of having a multiplication inside the loop (element of A times element of B), there would be an if statement (if element of A is unequal 0) and an addition (add element of B to the result). The if statement could be moved to the outer loop I guess. So, I am not sure if that could save some time alltogether. But if there is no better approach, I might try to write this.

[Bruno Luong]

I doubt you can beat MATLAB matrix multiplication (highly optimized) with MEX file. The time depends on little on the number of operations, but memory copying, etc ... count.

Answer 1

[Walter Roberson]

It turns out to be faster to leave A as logical when you do the matrix multiplication.

Using logical indexing, or doing find() and adding only those elements, is much slower unless the occupancy fraction drops to below 10% as a rough estimate.

[Florian]

Are you sure it's faster using a logical matrix? On my system it seems to be faster with double. Using Matlab R2018. Did I miss something in my example below?

At = randn(512)>0;
Bt = randn(512);
tic;
for k=1:10000
  Bt = At*Bt;
end
toc
At = double(At);
tic;
for k=1:10000
  Bt = At*Bt;
end
toc
Elapsed time is 29.873778 seconds.
Elapsed time is 19.203822 seconds.

[Walter Roberson]

If you want addition of some elements selected by A, then you should be using A.*B rather than A*B

[Florian]

I think you got my question wrong. I want to compute the product A*B. However, because A is binary, this product can be interpreted as "the addition of elements of B". My question was, if this fact can be used somehow to speed up A*B.

Also, how does A.*B sum up elements of B? Isn't it just setting some of the elements to zero?.

[Walter Roberson]

Ah, I guess I must have mis-interpreted the question.

Answer 2

[Bruno Luong]

Here is two ways, it won't be faster than A*B as other has warned you

n = 512;
A = sprand(n,n,0.1) > 0 ;
B = rand(n,n);

m = size(A,1);
p = size(B,2);
[i,j] = find(A);
C = zeros(m,p);
tic
for k = 1:p
    C(:,k) = accumarray(i,reshape(B(j,k),[],1),[m,1]);
end
toc

m = size(A,1);
p = size(B,2);
[i,j] = find(A);
k = 1:p;
[I,K] = ndgrid(i,k);
tic
C = accumarray([I(:),K(:)],reshape(B(j,k),[],1));
toc

tic
C = A*B;
toc