Why is this so slow?

Hi, All -- the seemingly simple script shown below is exceedingly slow. It takes several minutes on a fast workstation. I can't figure out why it takes so long. What am I missing? Thanks, in advance.

N = 10^6;
a = 0;
b = 4;
r = a + (b-a)*rand(N,1);
c = r.*sinh(r);
k = (1:N)';
s(N,1) = 0;
for j = 1:N
    s(j,1) = sum(c(1:j));
end
s = (b-a)*s./k;

Answer 1

[madhan ravi]

Instead of loop why not simply use cumsum() ?

[Moe Szyslak]

I was not aware of cumsum() -- that's blazing fast! That said, I'm still not sure why the loop is so very slow....

[Rik]

If you want to implement it with a loop, you can still gain a lot of speed:

N=10^4;
c=rand(N,1);

tic
s=[];s(N,1)=0;
for n=1:N
    s(n)=sum(c(1:n));
end
t1=toc;
tic
s=[];s(N,1)=0;
s(1)=c(1);
for n=2:N
    s(n)=s(n-1)+c(n);
end
t2=toc;
tic
s=cumsum(c);
t3=toc;
clc
fprintf('fast loop is %.0f times faster\ncumsum is %.0f times faster\n',...
    t1/t2,t1/t3)

That is over 100 times faster on my machine, with cumsum increasing the speed 10 times over that.

[Moe Szyslak]

That's much cleaner. I see similar speedups on my machine. Thanks again.

Answer 2

[Rik]

Loops with many elements are slow. You also forgot to properly pre-allocate s.

N = 10^5;
a = 0;
b = 4;
r = a + (b-a)*rand(N,1);
c = r.*sinh(r);
k = (1:N)';
%s(N,1) = 0;
s=zeros(N,1);
for j = 1:N
    s(j,1) = sum(c(1:j));
end
s = (b-a)*s./k;

clc
s2 = (b-a)*cumsum(c)./k;
fprintf('largest error = %.2e\n',max(abs(s-s2)))

[madhan ravi]

Yes Rik , I once asked the question about the syntax variable(N,1) = 0 is the same as variable = zeros(N,1)

[Moe Szyslak]

I thought that x(n,1)=0 was equivalent to x=zeros(n,1)? I also seem to recall reading somewhere that the former is actually a little faster than the latter, but maybe I'm wrong?

[Rik]

It works if the variable doesn't exist yet, but if it already exists (e.g. because you aborted a 10^6 run to do a 10^5 instead), this will cause a dimension error with this code.

I didn't remember the conclusion from your question, so I looked it up.

%replace this
s(N,1) = 0;
%with this
s=[];s(N,1) = 0;

Because apparently that is faster.

Except in this case, where it is faster to create the array with cumsum instead.

[madhan ravi]

Yes Rik