How to avoid large matrix multiplication in Matlab

Question

There are two large matrices in my code, which have the same number of columns and different number of rows. Like A(20000X4000) and B(30000X4000). Both are 0-1 sparse.

I should check each row of A with all rows of B and count the number of common 1s. For example if A(1,:)=[0 1 0 1 1] and B([1 2],:)=[1 1 1 1 1;0 0 0 1 1], I need to get the result as 3 and 2.

Assume there is a large 0-1 matrix C(50000X4000) and its rows are labeled as either type A or type B. I should compare all rows of A and B together and enumerate 1s. If the number of 1s in each row of A and B are greater than some bounds, then I used those rows of A and B for the rest of calculation. So, I do not even need to store A and B, all I need is a list of pairs of indices of rows. Something like [(3,2),(3,5),...] which shows I should use the third row of A and the second row of B, also the third of A and the fifth of B and so on.

The first thing that came into my mind was A*B', it gives the correct result, but practically it is very costly and in some cases impossible to do this multiplication.

I converted the matrices into single data type, and it became a bit faster. Sparse did not help.

The task seems easy, just counting common 1s of each row of A and all rows of B, but not easy to implement. Considering that the code should do this task like 1000 times, then it is practically impossible.

Any idea how to do enumerate the common ones without multiplication? (btw, loops also did not help).

Thanks.

Answer 1

If multiplication of the whole matrices cannot be done, one idea is to process a vertical stripe at a time. For each stripe you compute the desired result, and accumulate it with that of the preceding stripes:

A = double(rand(5,300)<.5); %// random data
B = double(rand(4,300)<.5); %// random data

S = 10; %// stripe size
result = zeros(size(A,1),size(B,1)); %// initialize to 0
for s = 1:10:size(A,2) %// each vertical stripe, of width S
    ind = s+(0:S-1);
    result = result +  A(:,ind)*(B(:,ind)).';
end

Check:

>> result

result =

    73    72    62    72
    84    70    79    71
    83    84    76    77
    77    80    77    74
    71    71    70    74

>> A*B.'

ans =

    73    72    62    72
    84    70    79    71
    83    84    76    77
    77    80    77    74
    71    71    70    74

Answer 2

The solution that you tried is optimal or near optimal.

When I try this, it takes less than a minute:

A = round(rand(30000,4000));
B = round(rand(20000,4000));
tic,A*B';toc;

If you really need to do this thousands of times, there are only two scenarios that I can imagine:

1. You don't need to do it often, in that case just let it run and it will be done tomorrow
2. You want to do it often and speed matters, finding a much faster solution is simply not going to happen. Unless you have some very usefull information about the matrices that you will be multiplying.

If you find that this sample multiplication much more than a minute (say more than 10 minutes), you probably are using memory inefficiently. In that case try to get some more ram.

Answer 3

I do not know if this is really any better than what you have, because it does still have a for loop in it, but if someone can figure out how to remove that for loop you should be good to go.

%  create temp data
A = rand(20000,4000) < 0.5;
B = rand(30000,4000) < 0.5;
counts = zeros(size(A,1),size(B,1),'uint8');
for i = 1:size(A,1)
    counts(i,:) = sum(bsxfun(@eq,A(i,:),B),2);
end

Either way, the process is going to take a long time because you are comparing 30000 rows with 4000 elements each, 20000 times, or approximately 2.4e+12 comparisons. That is a huge task, and will definitely take a long time. Possibly try using parallel computing if you need it to be faster.

Answer 4

I did some benchmarking; on my machine (i7-3770 @ 3.40GHz), multiplying full matrices of size 30000 x 4000 and 4000 x 20000 takes about 55 seconds regardless of content, the same that Dennis Jaheruddin found. But using sparse matrices can make the computation faster, depending on sparseness. If I define the degree of sparseness r as the ratio between the number of 1s to elements of the matrix, I get the following results:

r      time / s 
0.001   0.07
0.002   0.3
0.005   2.1
0.01    8.3
0.02   25

Here is the code used to determine these numbers:

m = 20000;
n = 4000;
o = 30000;

r = 0.001;

N = round(r * m * n);
A = sparse(randi(m, N, 1), randi(n, N, 1), 1, m, n);

N = round(r * n * o);
B = sparse(randi(o, N, 1), randi(n, N, 1), 1, o, n);

tic
C = A * B';
toc