Speed up this bottleneck line of code

Hello.

I have a big chunk of code, here is a part of it (arrays a, b and c, indeces ind1 and ind2 are actually different, its just an artificially constructed example):

n = 10000;
m = 5000;
k = 20;

a = rand(n,m);
b = rand(n,m);
c = rand(m+k,n)*1.1-0.1;

ind1 = repmat([ones(1,m/k),0]==1,1,k)';
ind2 = repmat([0,ones(1,m/k)]==1,1,k)';

tic;
d = m * sum(a .* b .* max(0, c(ind1,:) - c(ind2,:))',2);
toc

I am not satisfied with it's speed. After profiling my code i found that string

d = m * sum(a .* b .* max(0, c(ind1,:) - c(ind2,:))',2);

is a bottleneck of my code. Indeces ind1 and ind2 are almost entirely true.

Now im not sure now to improve from here, i've tried to google how to speed up .* operation or speed up indexing of matrix c but with no luck. Any advice would be great.

Thanks!

Answer 1

[Matt J]

If you can construct c in transposed form, there are some savings to be had,

n = 10000;
m = 5000;
k = 20;


a = rand(n,m);
b = rand(n,m);
c = rand(m+k,n)*1.1-0.1;

tic;
s=m/k+1;
[ind1,ind2]=deal((1:m+k)');
ind1(s:s:end)=[];
ind2(1:s:end)=[];
toc

Elapsed time is 0.004463 seconds.

tic;
d1 = m * sum(a .* b .* max(0, c(ind1,:)-c(ind2,:))',2);
toc

Elapsed time is 0.555621 seconds.

ct=c';
I=[(1:m)',(1:m)'];
S=repmat(m*[1,-1],m,1);

tic;
J=[ind1,ind2]; 
A=sparse(J,I,S,m+k,m);
d2 = sum( a.*b.*max(0,ct*A) ,2);
toc

Elapsed time is 0.347437 seconds.

percentDifference = norm(d1-d2)/norm(d1)*100

percentDifference = 4.0847e-14

[Alex]

Thanks!

It helps to cut off about 1/3 of time from this line!

Any suggestions about .* operator? Is it optimizable?

[Matt J]

No, all the basic matrix operators are very well optimized by MathWorks.

Answer 2

[Walter Roberson]

n = 10000;
m = 5000;
k = 20;

tic
a = rand(n,m);
b = rand(n,m);
c = rand(m+k,n)*1.1-0.1;
toc

Elapsed time is 1.545849 seconds.

ind1 = repmat([ones(1,m/k),0]==1,1,k)';
ind2 = repmat([0,ones(1,m/k)]==1,1,k)';

tic;
d = m * sum(a .* b .* max(0, c(ind1,:) - c(ind2,:))',2);
toc

Elapsed time is 0.560316 seconds.

Your bottleneck is in your generation of a, b, and c -- not in your calculation of d.

[Walter Roberson]

You would possibly be slightly faster if you were to transpose a, b, and c, so that your logical indexing is by column instead of by row. Selecting entire columns allows using compact instructions to copy blocks of memory; selecting by row requires examination of the indices at each point, element-by-element, because memory storage order proceeds row by row.

[Walter Roberson]

ind1 is equivalent to selecting 1:m/k . ind2 is equivalent to selecting 2:m/k+1. It might potentially be faster to use the numeric indices instead of logical indices.

[Alex]

What i mean by

"arrays a, b and c, indeces ind1 and ind2 are actually different, its just an artificially constructed example"

is that i dont post this part of code here, this code is actually different and optimized. These a, b, c, ind1 and ind2 are generated just to demonstrate calculations in row that is a bottleneck.