Thread Subject: Help speeding up/replacing loops

I've written some code, however it takes forever (several days) to run for larger values. Can anyone with more experience in matlab help me with the issue?? I understand vectorization is the way forward, but I'm unaware of how to do this for most of my loops.
Many thanks

Current Code:

function [] = the021testgraphnew(a1,an,aacc,x,N1,N2,g)
N = N1-N2;
ncut = N/10;
b = 0;
count = 0;
'creates the vector (1,2,...,ncut)';
X = 1:1:ncut;
'creates a vector with g random values';
C(1:g) = (rand(1)*pi);
'loops everything for each value of a (aka mu), this is time consuming but I cannot see another way around it';
for a = a1:aacc:an
    b = b + 1;
    A(b) = a;
    count = count + 1;
    'resets the expectation to 0';
    E = 0;
    'calculates the values of the transient';
    for j = 1:N2
        x = a*x*(1-x);
    end
    'calculates the values to actually be used and the expectation';
    for j = 1:N
        x = a*x*(1-x);
        B(j) = x;
        E = E + x;
    end
    E = E/N;
    'calculates the p''s and q''s';
    for f = 1:g
        p(1) = B(1)*cos(C(f));
        q(1) = B(1)*sin(C(f));
        for n = 2:N
                p(n) = p(n-1) + B(n)*cos(n*C(f));
                q(n) = q(n-1) + B(n)*sin(n*C(f));
        end
        'calculates the m(n)s for the range of c';
        for n = 1:ncut
            m = 0;
            for j = 1:(N-n)
                m = m + (p(n+j)-p(j))^2 + (q(n+j)-q(j))^2;
            end
            m = m/(N-n);
            Vosc = (E^2)*((1-cos(n*C(f)))/(1-cos(C(f))));
            D(n) = m - Vosc;
        end
        'creates a vector containing the values of k for each value of c';
        K(f) = corr(X,D);
    end
    'calculates the value of k to be plotted and stores it';
    k = median(K);
    KC(b) = k;
    'outputs how far through the program you''ve gotten so far';
    percents = ((count*aacc)/(an-a1))*100
end
'plots k against mu';
plot(A,KC)

I haven't looked over your code thoroughly, but one thing stands out immediately. Are you creating strings for comments? If so, you would be much better off to do this:

% A comment in the code.

instead of:

'A comment in the code';


MATLAB will ignore the true comment, but process the string, which requires time and resources.

"James Wright" <jameswright1001@yahoo.co.uk> wrote in message
news:h18ge7$gta$1@fred.mathworks.com...
> I've written some code, however it takes forever (several days) to run for
> larger values. Can anyone with more experience in matlab help me with the
> issue?? I understand vectorization is the way forward, but I'm unaware of
> how to do this for most of my loops.
> Many thanks
>
> Current Code:
>
> function [] = the021testgraphnew(a1,an,aacc,x,N1,N2,g)
> N = N1-N2;
> ncut = N/10;
> b = 0;
> count = 0;
> 'creates the vector (1,2,...,ncut)';
> X = 1:1:ncut;
> 'creates a vector with g random values';
> C(1:g) = (rand(1)*pi);
> 'loops everything for each value of a (aka mu), this is time consuming but
> I cannot see another way around it';
> for a = a1:aacc:an
> b = b + 1;
> A(b) = a;

MLINT should be flagging this line, and many of the other lines where you
create matrices inside of loops, and telling you to preallocate. Try
preallocating your matrices outside the loops and see if that improves the
performance.

MLINT:
http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_env/f9-11863.html

Preallocation:
http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_prog/f8-784135.html

*snip*

--
Steve Lord
slord@mathworks.com

"Steven Lord" <slord@mathworks.com> wrote in message <h18irl$1um$1@fred.mathworks.com>...
>
> "James Wright" <jameswright1001@yahoo.co.uk> wrote in message
> news:h18ge7$gta$1@fred.mathworks.com...
> > I've written some code, however it takes forever (several days) to run for
> > larger values. Can anyone with more experience in matlab help me with the
> > issue?? I understand vectorization is the way forward, but I'm unaware of
> > how to do this for most of my loops.
> > Many thanks
> >
> > Current Code:
> >
> > function [] = the021testgraphnew(a1,an,aacc,x,N1,N2,g)
> > N = N1-N2;
> > ncut = N/10;
> > b = 0;
> > count = 0;
> > 'creates the vector (1,2,...,ncut)';
> > X = 1:1:ncut;
> > 'creates a vector with g random values';
> > C(1:g) = (rand(1)*pi);
> > 'loops everything for each value of a (aka mu), this is time consuming but
> > I cannot see another way around it';
> > for a = a1:aacc:an
> > b = b + 1;
> > A(b) = a;
>
> MLINT should be flagging this line, and many of the other lines where you
> create matrices inside of loops, and telling you to preallocate. Try
> preallocating your matrices outside the loops and see if that improves the
> performance.
>
> MLINT:
> http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_env/f9-11863.html
>
> Preallocation:
> http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_prog/f8-784135.html
>
> *snip*
>
> --
> Steve Lord
> slord@mathworks.com
>
MLINT does flag them up and say that I should take them out of the loops, but that's what I'm unsure how to do. I've done it for simple things like C(1:g) = (rand(1)*pi); but I'm unsure how to do it for loops like:

p(1) = B(1)*cos(C(f));
q(1) = B(1)*sin(C(f));
for n = 2:N
      p(n) = p(n-1) + B(n)*cos(n*C(f));
      q(n) = q(n-1) + B(n)*sin(n*C(f));
end

where previous values in the vector are used to calculate the new one

Is it even possible to replace loops like:

p(1) = B(1)*cos(C(f));
q(1) = B(1)*sin(C(f));
for n = 2:N
      p(n) = p(n-1) + B(n)*cos(n*C(f));
      q(n) = q(n-1) + B(n)*sin(n*C(f));
end


??
I'm starting to think that maybe it's not

"James Wright" <jameswright1001@yahoo.co.uk> wrote in message <h18lud$rpl$1@fred.mathworks.com>...
> Is it even possible to replace loops like:
>
> p(1) = B(1)*cos(C(f));
> q(1) = B(1)*sin(C(f));
> for n = 2:N
> p(n) = p(n-1) + B(n)*cos(n*C(f));
> q(n) = q(n-1) + B(n)*sin(n*C(f));
> end
>
>
> ??
> I'm starting to think that maybe it's not



P = cumsum(B(1:N).*cos((1:N)*C(f)));
Q = cumsum(B(1:N).*sin((1:N)*C(f)));

"Matt Fig" <spamanon@yahoo.com> wrote in message <h18p3l$3u1$1@fred.mathworks.com>...
>
>
>
> P = cumsum(B(1:N).*cos((1:N)*C(f)));
> Q = cumsum(B(1:N).*sin((1:N)*C(f)));

Thank you so much, I'd been playing around with the cumsum trying to get it to work, but couldn't see how. Could you also help with vectorizing this?

for j = 1:N
    x = a*x*(1-x);
    B(j) = x;
    E = E + x;
end


Obviously I could rewrite this as,

B(1) = a*x*(1-x);
E = B(1);
for j = 2:N
     B(j) = a*(B(j-1))*(1-(B(j-1)));
     E = E + B(j);
end


But to vectorize it would I have to use the cumprod command?

"James Wright" <jameswright1001@yahoo.co.uk> wrote in message
news:h18kfh$j71$1@fred.mathworks.com...
> "Steven Lord" <slord@mathworks.com> wrote in message
> <h18irl$1um$1@fred.mathworks.com>...

*snip*

> MLINT does flag them up and say that I should take them out of the loops,
> but that's what I'm unsure how to do. I've done it for simple things like
> C(1:g) = (rand(1)*pi); but I'm unsure how to do it for loops like:
>
> p(1) = B(1)*cos(C(f));
> q(1) = B(1)*sin(C(f));
> for n = 2:N
> p(n) = p(n-1) + B(n)*cos(n*C(f));
> q(n) = q(n-1) + B(n)*sin(n*C(f));
> end
>
> where previous values in the vector are used to calculate the new one

You may not need to vectorize this at all if you preallocate. To
preallocate, just put values into the array before the loop until it's the
size that it will need to be in the end.

p = zeros(1, N);
q = zeros(1, N);
p(1) = B(1)*cos(C(f));
q(1) = B(1)*sin(C(f));
for n = 2:N
     p(n) = p(n-1) + B(n)*cos(n*C(f));
     q(n) = q(n-1) + B(n)*sin(n*C(f));
end

Now instead of growing the p and q arrays inside the loop (making it 2
elements long, then 3, then 4, etc.) you're just filling them in. That
means no memory growing inside the loop, so it'll be faster.

If preallocation doesn't improve performance enough, then you can continue
on to vectorization using CUMSUM as Matt suggested.

--
Steve Lord
slord@mathworks.com

Thanks for all the help so far, really appreciate it! This is an update of what my code looks like now, which I think is quite a bit better than before, but there's still a few points I could do with some help on.

function [] = the021test(a1,an,aacc,x,N1,N2,g)
N = N1-N2;
ncut = N/10;
b = 0;
count = 0;
kc = ((an - a1)/aacc);
%allocating vectors with the correct size, for memory reasons
B = zeros(1,N);
K = zeros(1,g);
m = zeros(1,ncut);
KC = zeros(1,kc);
A = zeros(1,kc);
%creates the vector (1,2,...,ncut)
X = 1:1:ncut;
%creates a vector with g random values
C(1:g) = (rand(1)*pi);
%loops everything for each value of a (aka mu), this is time consuming but
%I cannot see another way around it
for a = a1:aacc:an
    b = b + 1;
    A(b) = a;
    count = count + 1;
    %calculates the values of the transient
    for j = 1:N2
        x = a*x*(1-x);
    end
    %calculates the values to actually be used and the expectation
    B(1) = a*x*(1-x);
    for j = 2:N
        B(j) = a*(B(j-1))*(1-(B(j-1)));
    end
    E = mean(B);
    %calculates the p's and q's
    for f = 1:g
        p = cumsum(B(1:N).*cos((1:N)*C(f)));
        q = cumsum(B(1:N).*sin((1:N)*C(f)));
        %calculates the m(n)s for the range of c
        for n = 1:ncut
            for j = 1:(N-n)
                m(n) = m(n) + (p(n+j)-p(j))^2 + (q(n+j)-q(j))^2;
            end
            m(n) = (m(n))/(N-n);
        end
        %calculates the oscillatory terms
        Vosc(1:n) = (E^2)*((1-cos((1:ncut)*C(f)))/(1-cos(C(f))));
        %Works out the new meansquare displacement with the oscillating
        %term taken out
        D(1:ncut) = m(1:ncut) - Vosc(1:ncut);
        %creates a vector containing the values of k for each value of c
        K(f) = corr(X,D);
    end
    %calculates the value of k to be plotted and stores it
    k = median(K);
    KC(b) = k;
    %outputs how far through the program you've gotten so far
    percents = ((count*aacc)/(an-a1))*100
end
%plots k against mu
plot(A,KC)



I'm wondering how I would vectorize

B(1) = a*x*(1-x);
for j = 2:N
     B(j) = a*(B(j-1))*(1-(B(j-1)));
end

Sorry if it seems like I'm asking for help on something that seems the same, but the problem I'm having is because B(j-1) is used twice. Is there any way to sort this out?
Thanks

"James Wright" <jameswright1001@yahoo.co.uk> writes:

>Thanks for all the help so far, really appreciate it! This is an update of what my code looks like now, which I think is quite a bit better than before, but there's still a few points I could do with some help on.

>I'm wondering how I would vectorize
I think you should worry most about the loops within loops within loops.

>B(1) = a*x*(1-x);
>for j = 2:N
> B(j) = a*(B(j-1))*(1-(B(j-1)));
>end

Try (I think)

Bshifted=circshift(B(1:N),[1,1]);
B(1:N)=a.*Bshifted.*(1-Bshifted);
B(1) = a*x*(1-x);

tpl@eng.cam.ac.uk (Tim Love) wrote in message <h1apo1$kjc$1@gemini.csx.cam.ac.uk>...
> I think you should worry most about the loops within loops within loops.
>
> >B(1) = a*x*(1-x);
> >for j = 2:N
> > B(j) = a*(B(j-1))*(1-(B(j-1)));
> >end
>
> Try (I think)
>
> Bshifted=circshift(B(1:N),[1,1]);
> B(1:N)=a.*Bshifted.*(1-Bshifted);
> B(1) = a*x*(1-x);
This doesn't seem to give the same values for the elements of B

I don't see an easy vectorization of this snippet. However, I am wondering why you have not taken the good advice to pre-allocate your arrays before the loops? This can lead to speed increases which rival those of vectorization in some cases! For example, put this in a file and run the file:


clear
N = 10000;

tic
B(1) = 1; % No Pre-allocation
for ii = 2:N
    B(ii) = ii + ii.^2;
end
toc

tic
B2 = ones(1,N); % Pre-allocation
for jj = 2:N
   B2(jj) = jj + jj.^2;
end
toc

isequal(B,B2)

"James Wright" <jameswright1001@yahoo.co.uk> writes:


>This doesn't seem to give the same values for the elements of B
Sorry, I posted bogus code. I think I'd better leave it to others
to post something that works.

"Matt Fig" <spamanon@yahoo.com> wrote in message <h1b1j5$kks$1@fred.mathworks.com>...
> I don't see an easy vectorization of this snippet. However, I am wondering why you have not taken the good advice to pre-allocate your arrays before the loops? This can lead to speed increases which rival those of vectorization in some cases! For example, put this in a file and run the file:
>
>
> clear
> N = 10000;
>
> tic
> B(1) = 1; % No Pre-allocation
> for ii = 2:N
> B(ii) = ii + ii.^2;
> end
> toc
>
> tic
> B2 = ones(1,N); % Pre-allocation
> for jj = 2:N
> B2(jj) = jj + jj.^2;
> end
> toc
>
> isequal(B,B2)
I already did the pre-allocations, they're just at the very top of the program. I'm foolishly assuming I can pre-allocate them all at the start, would it actually be better to pre-allocate right before the loop?

Maybe finding a way to vectorize that isn't worth the time it'll take. The main time consumer is most probably this loop anyway.

    for f = 1:g
        p = cumsum(B(1:N).*cos((1:N)*C(f)));
        q = cumsum(B(1:N).*sin((1:N)*C(f)));
        %calculates the m(n)s for the range of c
        for n = 1:ncut
            m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
            m(n) = (m(n))/(N-n);
        end
        %calculates the oscillatory terms
        Vosc(1:n) = (E^2)*((1-cos((1:ncut)*C(f)))/(1-cos(C(f))));
        %Works out the new meansquare displacement with the oscillating
        %term taken out
        D(1:ncut) = m(1:ncut) - Vosc(1:ncut);
        %creates a vector containing the values of k for each value of c
        K(f) = corr(X,D);
    end


However, is it possible to vectorize this further? As the loops already contain vectorizations. I imagine this would also make things extremely complicated.

Thanks for all the help so far!

"James Wright" <jameswright1001@yahoo.co.uk> wrote in message
> I already did the pre-allocations, they're just at the very top of the program. I'm foolishly assuming I can pre-allocate them all at the start, would it actually be better to pre-allocate right before the loop?

I apologize, I didn't see that you had done this. As long as the sizes of the variables aren't changing every time through the big outer loop, you are fine to do it there.
Any more help I could give would depend on the input parameters. If you post everything needed to run the code the way you normally would, I might be able to offer more.

"Matt Fig" <spamanon@yahoo.com> wrote in message <h1b5er$8td$1@fred.mathworks.com>...
> I apologize, I didn't see that you had done this. As long as the sizes of the variables aren't changing every time through the big outer loop, you are fine to do it there.
> Any more help I could give would depend on the input parameters. If you post everything needed to run the code the way you normally would, I might be able to offer more.


Ok, the entirety of my program is:

function [] = the021testgraphvectorized(a1,an,aacc,x,N1,N2,g)
N = N1-N2;
ncut = N/10;
b = 0;
count = 0;
kc = ((an - a1)/aacc);
%allocating vectors with the correct size, for memory reasons
B = zeros(1,N);
K = zeros(1,g);
D = zeros(1,ncut);
m = zeros(1,ncut);
KC = zeros(1,kc);
A = zeros(1,kc);
%creates the vector (1,2,...,ncut)
X = 1:1:ncut;
%creates a vector with g random values
C(1:g) = (rand(1)*pi);
%loops everything for each value of a (aka mu), this is time consuming but
%I cannot see another way around it
for a = a1:aacc:an
    b = b + 1;
    A(b) = a;
    count = count + 1;
    %calculates the values of the transient
    %these are not used in the main calculation
    for j = 1:N2
        x = a*x*(1-x);
    end
    %calculates the values of x to actually be used and stores them
    B(1) = a*x*(1-x);
    for j = 2:N
        B(j) = a*(B(j-1))*(1-(B(j-1)));
    end
    %calculates the expectation
    E = mean(B);
    %calculates the p's and q's
    for f = 1:g
        p = cumsum(B(1:N).*cos((1:N)*C(f)));
        q = cumsum(B(1:N).*sin((1:N)*C(f)));
        %calculates the m(n)s for the range of c, m(n) is a vector containing values
        %of the meansquare displacement
        for n = 1:ncut
            m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
            m(n) = (m(n))/(N-n);
        end
        %calculates the oscillatory terms
        Vosc(1:n) = (E^2)*((1-cos((1:ncut)*C(f)))/(1-cos(C(f))));
        %Works out the new meansquare displacement with the oscillating
        %term taken out
        D(1:ncut) = m(1:ncut) - Vosc(1:ncut);
        %creates a vector containing the values of k for each value of c
        K(f) = corr(X,D);
    end
    %takes the median value of K as the chosen value for k to be plotted and stores it
    k = median(K);
    KC(b) = k;
    %outputs how far through the program you've gotten
    percents = ((count*aacc)/(an-a1))*100
end
%plots k against mu
plot(A,KC)

when running it I'm using the inputs (3.5,4,0.001,0.1,1050000,990000,150)


where it says "K(f) = corr(X,D)" in the code that calls:

function [k] = corr(x,y)
covxy = covar(x,y);
varx = covar(x,x);
vary = covar(y,y);
k = covxy/(sqrt(varx*vary));


which in turn is calling:

function [cov] = covar(x,y)
q = length(x);
xbar = mean(x);
ybar = mean(y);
cov = sum((x(1:q)-xbar).*(y(1:q)-ybar));
cov = cov/q;


I did try to eliminate the need for these programs by using matlabs "corrcoef" command which created a matrix with the correlation, but this seemed to muck things up somewhere along the line. If this is what you meant, hope it helps. I understand that there's a lot of code here and it's not the easiest to understand (due to my poor coding skills) and I can't thank you enough the help.

By far your major time hog is this:

for n = 1:ncut
     m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
     m(n) = (m(n))/(N-n);
end

Much of the time is spent indexing into p and q. As such I see no way of dramatically speeding this up without creating large intermediate arrays (because your vectors are large to begin with). Now maybe someone else sees it, but I don't. I think your only choice is to go to MEX. Good luck.

"Matt Fig" <spamanon@yahoo.com> wrote in message <h1bc89$83t$1@fred.mathworks.com>...
> By far your major time hog is this:
>
> for n = 1:ncut
> m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
> m(n) = (m(n))/(N-n);
> end
>
> Much of the time is spent indexing into p and q. As such I see no way of dramatically speeding this up without creating large intermediate arrays (because your vectors are large to begin with). Now maybe someone else sees it, but I don't. I think your only choice is to go to MEX. Good luck.
Thanks

On Jun 17, 2:48 pm, "James Wright" <jameswright1...@yahoo.co.uk>
wrote:
> "Matt Fig" <spama...@yahoo.com> wrote in message <h1bc89$83...@fred.mathworks.com>...
> > By far your major time hog is this:
>
> > for n = 1:ncut
> >      m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
> >      m(n) = (m(n))/(N-n);
> > end
>
> > Much of the time is spent indexing into p and q.  As such I see no way of dramatically speeding this up without creating large intermediate arrays (because your vectors are large to begin with).  Now maybe someone else sees it, but I don't.  I think your only choice is to go to MEX.  Good luck.
>
> Thanks

Hi, James.

I was looking at your code and I just noticed the following lines:

%creates a vector with g random values
C(1:g) = (rand(1)*pi);

Are you aware that all of the entries of C are identical? Is that
your intention? If so, you could eliminate needless indexing into C
in the remainder of your code by replacing C with a scalar. If not,
perhaps you meant to write

C = rand(1,g) * pi;

--Peter

Peter <petersamsimon2@hotmail.com> wrote in message <5a16f286-988c-4681-92d0-33f9a19062a8@d25g2000prn.googlegroups.com>...
> On Jun 17, 2:48?pm, "James Wright" <jameswright1...@yahoo.co.uk>
> wrote:
> > "Matt Fig" <spama...@yahoo.com> wrote in message <h1bc89$83...@fred.mathworks.com>...
> > > By far your major time hog is this:
> >
> > > for n = 1:ncut
> > > ? ? ?m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
> > > ? ? ?m(n) = (m(n))/(N-n);
> > > end
> >
> > > Much of the time is spent indexing into p and q. ?As such I see no way of dramatically speeding this up without creating large intermediate arrays (because your vectors are large to begin with). ?Now maybe someone else sees it, but I don't. ?I think your only choice is to go to MEX. ?Good luck.
> >
> > Thanks
>
> Hi, James.
>
> I was looking at your code and I just noticed the following lines:
>
> %creates a vector with g random values
> C(1:g) = (rand(1)*pi);
>
> Are you aware that all of the entries of C are identical? Is that
> your intention? If so, you could eliminate needless indexing into C
> in the remainder of your code by replacing C with a scalar. If not,
> perhaps you meant to write
>
> C = rand(1,g) * pi;
>
> --Peter
Thanks, I wasn't aware that they were all the same

On Jun 18, 6:38 am, "James Wright" <jameswright1...@yahoo.co.uk>
wrote:
> Peter <petersamsim...@hotmail.com> wrote in message <5a16f286-988c-4681-92d0-33f9a1906...@d25g2000prn.googlegroups.com>...
> > On Jun 17, 2:48?pm, "James Wright" <jameswright1...@yahoo.co.uk>
> > wrote:
> > > "Matt Fig" <spama...@yahoo.com> wrote in message <h1bc89$83...@fred.mathworks.com>...
> > > > By far your major time hog is this:
>
> > > > for n = 1:ncut
> > > > ? ? ?m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
> > > > ? ? ?m(n) = (m(n))/(N-n);
> > > > end
>
> > > > Much of the time is spent indexing into p and q. ?As such I see no way of dramatically speeding this up without creating large intermediate arrays (because your vectors are large to begin with). ?Now maybe someone else sees it, but I don't. ?I think your only choice is to go to MEX. ?Good luck.
>
> > > Thanks
>
> > Hi, James.
>
> > I was looking at your code and I just noticed the following lines:
>
> > %creates a vector with g random values
> > C(1:g) = (rand(1)*pi);
>
> > Are you aware that all of the entries of C are identical?  Is that
> > your intention?  If so, you could eliminate needless indexing into C
> > in the remainder of your code by replacing C with a scalar. If not,
> > perhaps you meant to write
>
> > C = rand(1,g) * pi;
>
> > --Peter
>
> Thanks, I wasn't aware that they were all the same

Hi, James.

I added the line

cj = sqrt(-1);

at the start of your function, and changed the innermost and next
outer enclosing loop as follows:

    Esq = E^2; % remove loop invariant calculation from below loop
    %calculates the p's and q's
    for f = 1:g
        q = exp(C(f) * cj * (1:N));
        p = cumsum(B .* q);
        %calculates the m(n)s for the range of c, m(n) is a vector
containing values
        %of the meansquare displacement
        for n = 1:ncut
            m(n) = norm(p((n+1):N) - p(1:(N-n)))^2 / (N-n);
        end
        %calculates the oscillatory terms
        Vosc = Esq*((1-real(q(1:ncut)))/(1-real(q(1))));
        %Works out the new meansquare displacement with the
oscillating
        %term taken out
        D = m - Vosc;
        %creates a vector containing the values of k for each value of
c
        K(f) = corr(X,D);
    end


The sines and cosines are evaluated in a single call to the
exponential function with imaginary argument and stored in q. p now
contains the old p values in its real part and the old q values in its
imaginary part. This saves on redundant indexing and also allows the
cosine values to be reused when computing Vosc without recomputing
them each time in the outer loop. I ran your original version and this
modified version using the following call:

tic; the021testgraphvectorized(3.5,4,0.001,0.1,1250,990,100); toc

The original code required 53.8 seconds and my modified version
required 43.5 seconds, a time savings of about 19%. Not bad for a
small amount of effort. If your code takes days to run this could add
up to many hours of savings.


--Peter

Peter <petersamsimon2@hotmail.com> wrote in message <11be7337-b5a7-4449-89c8-b2545c852be9@y6g2000prf.googlegroups.com>...

> Hi, James.
>
> I added the line
>
> cj = sqrt(-1);
>
> at the start of your function, and changed the innermost and next
> outer enclosing loop as follows:
>
> Esq = E^2; % remove loop invariant calculation from below loop
> %calculates the p's and q's
> for f = 1:g
> q = exp(C(f) * cj * (1:N));
> p = cumsum(B .* q);
> %calculates the m(n)s for the range of c, m(n) is a vector
> containing values
> %of the meansquare displacement
> for n = 1:ncut
> m(n) = norm(p((n+1):N) - p(1:(N-n)))^2 / (N-n);
> end
> %calculates the oscillatory terms
> Vosc = Esq*((1-real(q(1:ncut)))/(1-real(q(1))));
> %Works out the new meansquare displacement with the
> oscillating
> %term taken out
> D = m - Vosc;
> %creates a vector containing the values of k for each value of
> c
> K(f) = corr(X,D);
> end
>
>
> The sines and cosines are evaluated in a single call to the
> exponential function with imaginary argument and stored in q. p now
> contains the old p values in its real part and the old q values in its
> imaginary part. This saves on redundant indexing and also allows the
> cosine values to be reused when computing Vosc without recomputing
> them each time in the outer loop. I ran your original version and this
> modified version using the following call:
>
> tic; the021testgraphvectorized(3.5,4,0.001,0.1,1250,990,100); toc
>
> The original code required 53.8 seconds and my modified version
> required 43.5 seconds, a time savings of about 19%. Not bad for a
> small amount of effort. If your code takes days to run this could add
> up to many hours of savings.
>
>
> --Peter
Thank you for your help again, I was considering using mex files to speed up the program, though I've never used them before, do you know if they're significantly faster? If I can't make my program run multiple times quicker than it currently is, I think I'm going to have to write it in C instead.

"James Wright" <jameswright1001@yahoo.co.uk> wrote in message
> Thank you for your help again, I was considering using mex files to speed up the program, though I've never used them before, do you know if they're significantly faster? If I can't make my program run multiple times quicker than it currently is, I think I'm going to have to write it in C instead.

MEX can be many times faster, depending on the problem. Check this out for example:

http://www.mathworks.com/matlabcentral/fileexchange/21376

This has a MATLAB solution and a MEX-File of the same solution, so you can compare the speeds for different sizes of the problem. Good luck.

On Jun 19, 3:44 pm, "Matt Fig" <spama...@yahoo.com> wrote:
> "James Wright" <jameswright1...@yahoo.co.uk> wrote in message
> > Thank you for your help again, I was considering using mex files to speed up the program, though I've never used them before, do you know if they're significantly faster? If I can't make my program run multiple times quicker than it currently is, I think I'm going to have to write it in C instead.
>
> MEX can be many times faster, depending on the problem.  Check this out for example:
>
> http://www.mathworks.com/matlabcentral/fileexchange/21376
>
> This has a MATLAB solution and a MEX-File of the same solution, so you can compare the speeds for different sizes of the problem.  Good luck.

One other thing you might want to consider is to simply bypass the MEX
route and write this as a standalone Fortran or C (or C++ or whatever
compiler you're most comfortable with) program. The program could
take its inputs from the command line and write its outputs to a
file. The overhead of the file I/O would be negligible compared to
the computation time. In particular, Fortran is very good at fast
numerical computation, and Fortran 95 syntax is very similar to Matlab
syntax, so converting this purely numerical computation to Fortran
would be very straightforward. If you wanted to, you could invoke the
compiled executable from within Matlab (using the system command) and
then gather up the ouputs into your Matlab session for further
processing/plotting. This method is easier than MEX coding, and for
your application I don't think the inefficiency associated with using
files to communicate with the Matlab session is an important
consideration.

--Peter

"Matt Fig" <spamanon@yahoo.com> wrote in message <h1bc89$83t$1@fred.mathworks.com>...
> By far your major time hog is this:
>
> for n = 1:ncut
> m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
> m(n) = (m(n))/(N-n);
> end
>

Regardless it will be implemented in Matlab or Mex, I believe the above loop can be rewitten as few cumulative sums. Just developing the square expression and try to see what to add from one iteration to the next. This will reduce the complexity from O(n^2) to O(n).

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h1i25t$6i$1@fred.mathworks.com>...
> "Matt Fig" <spamanon@yahoo.com> wrote in message <h1bc89$83t$1@fred.mathworks.com>...
> > By far your major time hog is this:
> >
> > for n = 1:ncut
> > m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
> > m(n) = (m(n))/(N-n);
> > end
> >
>
> Regardless it will be implemented in Matlab or Mex, I believe the above loop can be rewitten as few cumulative sums. Just developing the square expression and try to see what to add from one iteration to the next. This will reduce the complexity from O(n^2) to O(n).

 Sorry, it must need an additional autocorrelation calculation of p and also for q. That's (n*log(n)).

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h1imem$b5f$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h1i25t$6i$1@fred.mathworks.com>...
> > "Matt Fig" <spamanon@yahoo.com> wrote in message <h1bc89$83t$1@fred.mathworks.com>...
> > > By far your major time hog is this:
> > >
> > > for n = 1:ncut
> > > m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
> > > m(n) = (m(n))/(N-n);
> > > end
> > >
> >
> > Regardless it will be implemented in Matlab or Mex, I believe the above loop can be rewitten as few cumulative sums. Just developing the square expression and try to see what to add from one iteration to the next. This will reduce the complexity from O(n^2) to O(n).

Here is the practice:

% Data
N=10;
ncut=7;

p=10*rand(1,N);
q=10*rand(1,N);
m=zeros(1,N);

% Engine for-loop
for n = 1:ncut
    m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
    m(n) = (m(n))/(N-n);
end

% Vectorized Engine
mm=zeros(1,N);

n=1:ncut;
cp=conv(p,p(end:-1:1),'full'); % could do better
p2 = p.^2;
cp2 = cumsum(p2);
cp2f = cumsum(p2(N+1-n));

cq=conv(q,q(end:-1:1),'full'); % could do better
q2 = q.^2;
cq2 = cumsum(q2);
cq2f = cumsum(q2(N+1-n));

mm(n) = 2*cp2(end) - (cp2(n) + cp2f(n) + 2*cp(N-n)) + ...
        2*cq2(end) - (cq2(n) + cq2f(n) + 2*cq(N-n)); % could simplify a bit
mm(n) = mm(n) ./ (N-(n));

% Compare
m
mm

% Bruno

Matlab CONV does not seem to use FFT (contrary to the doc, since it is very slow). I use my own FFT-based convolution; and here is a timing.

function test

% Data
N=5e4;
ncut=1e4;

p=10*rand(1,N);
q=10*rand(1,N);
m=zeros(1,N);

tic
% Engine for-loop
for n = 1:ncut
    m(n) = sum((p((n+1):N)-p(1:(N-n))).^2+(q((n+1):N)-q(1:(N-n))).^2);
    m(n) = (m(n))/(N-n);
end
toc % Elapsed time is 22.610620 seconds.
if N<20
    disp(m);
end
clear m

% Vectorized Engine
tic
mm=zeros(1,N);

n=1:ncut;
cp=MultiConv1(p,p(end:-1:1),2);
p2 = p.^2;
cp2 = cumsum(p2);
cp2f = cumsum(p2(N+1-n));

cq=MultiConv1(q,q(end:-1:1),2);
q2 = q.^2;
cq2 = cumsum(q2);
cq2f = cumsum(q2(N+1-n));

mm(n) = 2*cp2(end) - (cp2(n) + cp2f(n) + 2*cp(N-n)) + ...
        2*cq2(end) - (cq2(n) + cq2f(n) + 2*cq(N-n));
mm(n) = mm(n) ./ (N-(n));
toc % Elapsed time is 0.089810 seconds.

if N<20
    disp(mm);
end

end

function C = MultiConv1(A, B, dim)
%
% Use for Polynomial multiplication
%
if nargin<3 || isempty(dim)
    dim=1;
end
m = size(A,dim);
n = size(B,dim);
l = 2^nextpow2(m+n-1);
C = ifft(fft(A,l,dim).* fft(B,l,dim),l,dim);
subs(1:ndims(C)) = {':'};
subs{dim} = 1:m+n-1;
C = real(C(subs{:})); % remove real for complex polynomial

end % MultiConv1

Thanks for all the help, I think I've finally got something that will be quick enough (still a few days but better than the 125 it would have taken before) and accurate enough.
Couldn't have got anywhere near as fast without all of you.
Thanks!

James

"James Wright" <jameswright1001@yahoo.co.uk> wrote in message <h1ngnh$68o$1@fred.mathworks.com>...
> Thanks for all the help, I think I've finally got something that will be quick enough (still a few days but better than the 125 it would have taken before) and accurate enough.
> Couldn't have got anywhere near as fast without all of you.
> Thanks!
>
> James
Maybe i'm just saying something u already know, have u ever tried to use the profiler?
It gives a summary on the relative weights that each piece of code has on the total amount of time needed to execute it.

profile on
% type ur code/function
profile off
profile viewer

for an accurate description look for "profiling for improving performance" in the matlab general help.