Thread Subject: avoiding loops to built a z matrix of the most efficient form

Hello us, and everyone.
us showed me how develop a matrix z
avoiding loops....but....I felt really surprised how the algorithm with loops
still run faster than the algorithm using a arrayfunction (without lopps).

Why?

You can see it in the next results:

%piece of code

clear all
clc
tic

% Way 1: Algorithm using a function arrayfun to avoid loops.


nrv=2 % random variable

nrp=3 % number of 2D-points points.

v1=rand(nrp,2); % 2D-points
v2=rand(nrv,2); % 2D-parameters



     [ix2,ix1]=meshgrid(1:nrv,1:nrp);
     r=arrayfun(@(x,y) v1(x,:)*v2(y,:)',ix1,ix2,'uni',false);
     m=sum(cell2mat(r)) % <- matrix
     
     toc

% way 2: Using my algorithm with loops.

tic

z1=zeros(1,nrp);
z1=zeros(1,nrv);


for j=1:nrv


for i=1:nrp
z1(i)=v1(i,:)*v2(j,:)';
end

z(j)=sum(z1)'

end

toc


Results:

way 1: without loops, using the arrayfun suggested by us:

m =

    1.7735 1.1724

Elapsed time is 0.041280 seconds.

way 2: With loops (my code)

z =

    1.7735 1.1724

Elapsed time is 0.000276 seconds.

You can see how using loops the algorithm run faster than the other without loops sugested by us.

Anyone can help me how calculate the z matrix avoiding loops in other alternative way?

It is a nightmare!!
Thanks,

Jose.

Loops have not been slow in Matlab for several releases now. OFTEN, running a loop is going to be faster than vectorizing the code. Not always, but it should definitely not be a surprise to find this is the case in any given situation.

Now, the real question. Why are you so determined to solve your problem without loops? Are you looking for a performance boost? Do you just not like loops? Is this a requirement for a homework assignment?

"Jose " <jose.l.vega@gmail.com> wrote in message <gqdj93$qa4$1@fred.mathworks.com>...
> Hello us, and everyone.
> us showed me how develop a matrix z
> avoiding loops....but....I felt really surprised how the algorithm with loops
> still run faster than the algorithm using a arrayfunction (without lopps).
> .......

  Jose, have you tried out the solution I sent you in the thread "Built a matrix z(3x2) avoiding loops" at

 http://www.mathworks.com/matlabcentral/newsreader/view_thread/247361

which appears to be related to the question in this thread? If so, how well did it perform timewise as against a for-loop method?

Roger Stafford

"Jose " <jose.l.vega@gmail.com> wrote in message <gqdj93$qa4$1@fred.mathworks.com>...

>
> You can see how using loops the algorithm run faster than the other without loops sugested by us.
>

I can't see the point to optimize/benchmark for such a small-size problem, unless if you repeat it many time, but we haven't told that.

The overhead of calling functions are expensive, and not computation. In this case, just stick with the for-loop.

Otherwise, for big size here is one way:


%piece of code
tic
nrv=1000% random variable
nrp=1001 % number of 2D-points points.

v1=rand(nrp,2); % 2D-points
v2=rand(nrv,2); % 2D-parameters

% way 2: Using my algorithm with loops.

tic
z1=zeros(1,nrp);
zloop=zeros(1,nrv);

for j=1:nrv
    for i=1:nrp
        z1(i)=v1(i,:)*v2(j,:)';
    end
    zloop(j)=sum(z1)';
end
toc % 1.235350 seconds.


% way three

tic
z=sum(reshape(sum(bsxfun(@times,...
      reshape(v1,[nrp 1 2]),...
      reshape(v2,[1 nrv 2])),3),...
      [nrp nrv]),1);
toc % 0.029425 seconds.

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gqdn0i$a5r$1@fred.mathworks.com>...
> "Jose " <jose.l.vega@gmail.com> wrote in message <gqdj93$qa4$1@fred.mathworks.com>...
>
> >
> > You can see how using loops the algorithm run faster than the other without loops sugested by us.
> >
>
> I can't see the point to optimize/benchmark for such a small-size problem, unless if you repeat it many time, but we haven't told that.
>
> The overhead of calling functions are expensive, and not computation. In this case, just stick with the for-loop.
>
> Otherwise, for big size here is one way:
>
>
> %piece of code
> tic
> nrv=1000% random variable
> nrp=1001 % number of 2D-points points.
>
> v1=rand(nrp,2); % 2D-points
> v2=rand(nrv,2); % 2D-parameters
>
> % way 2: Using my algorithm with loops.
>
> tic
> z1=zeros(1,nrp);
> zloop=zeros(1,nrv);
>
> for j=1:nrv
> for i=1:nrp
> z1(i)=v1(i,:)*v2(j,:)';
> end
> zloop(j)=sum(z1)';
> end
> toc % 1.235350 seconds.
>
>
> % way three
>
> tic
> z=sum(reshape(sum(bsxfun(@times,...
> reshape(v1,[nrp 1 2]),...
> reshape(v2,[1 nrv 2])),3),...
> [nrp nrv]),1);
> toc % 0.029425 seconds.
>
> % Bruno


Thanks very much guys : Roger, us, Matt and Bruno.


Here we have the answer:

This is the different codes, 1 ,2 and 3:

clear all
clc

1. suggested by us.

tic

nrv=1000; % random variable

nrp=1000; % number of 2D-points x(nrp,:)

v1=rand(nrp,2);
v2=rand(nrv,2);

inversa=[1 2; 3 4];

     [ix2,ix1]=meshgrid(1:nrv,1:nrp);
     
     r=arrayfun(@(x,y) (v1(x,:)-v2(y,:))*inversa*(v1(x,:)-v2(y,:))',ix1,ix2,'uni',false);
     %[r{1,:}]; % <- first row
     %[r{nrp,:}]; % <- last row
     m=sum(cell2mat(r)); % <- matrix
     
     toc

     
     
     
2. Suggested by jose. (The original one with loops).

tic

z1=zeros(1,nrp);
z=zeros(1,nrv);


for j=1:nrv


for i=1:nrp
z1(i)=(v1(i,:)-v2(j,:))*inversa*(v1(i,:)-v2(j,:))';
end

z(j)=sum(z1)';

end

toc


3. Suggested by Roger.



tic

 z = repmat(v1,nrv,1)-reshape(repmat(reshape(v2,1,[]),nrp,1),[],2);
 z = sum(reshape(sum(((z*inversa).*z),2),nrp,nrv));


toc

4. Suggested by Bruno, pending to test for unknowing the right function:
(v1(i,:)-v2(j,:))*inversa*(v1(i,:)-v2(j,:))'

% tic
% z=sum(reshape(sum(bsxfun(@times,...
% reshape(v1,[nrp 1 2]),...
% reshape(v2,[1 nrv 2])),3),...
% [nrp nrv]),1)
%
% toc

Result:

us: Elapsed time is 17.815049 seconds.
Jose: Elapsed time is 5.732546 seconds.
Roger Elapsed time is 0.273193 seconds.

I did not check the Bruno code, because He used the function v1*v2' and
i am interested in (v1-v2)*inversa*(v1-v2)'.
It would be interested if Bruno can reply with the correct function I am interested to to check it.

Then, the conclusion up to the moment is that using "rematting" and "reshapping" is the faster way to do it as Roger suggest when we consider big number of data v1 and parameters v2.

Thank you guys for your help.

Jose.

"Jose " <jose.l.vega@gmail.com> wrote in message <gqdr01$jp4$1@fred.mathworks.com>...
> ......
> Then, the conclusion up to the moment is that using "rematting" and "reshapping" is the faster way to do it as Roger suggest when we consider big number of data v1 and parameters v2.
> ......

  I don't think you have made a fair comparison here, Jose. The question asked here in this thread and answered by Urs and Bruno is different, though related, from the one I have answered in the other thread "Built a matrix z(3x2) avoiding loops". The latter thread involved the 'inversa' matrix and this thread didn't. You are comparing "apples and oranges" here.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gqdsgp$aic$1@fred.mathworks.com>...
> "Jose " <jose.l.vega@gmail.com> wrote in message <gqdr01$jp4$1@fred.mathworks.com>...
> > ......
> > Then, the conclusion up to the moment is that using "rematting" and "reshapping" is the faster way to do it as Roger suggest when we consider big number of data v1 and parameters v2.
> > ......
>
> I don't think you have made a fair comparison here, Jose. The question asked here in this thread and answered by Urs and Bruno is different, though related, from the one I have answered in the other thread "Built a matrix z(3x2) avoiding loops". The latter thread involved the 'inversa' matrix and this thread didn't. You are comparing "apples and oranges" here.
>
> Roger Stafford

Ok roger, "apples and oranges" in different thread, but "apples with apples" in
the test I showed above with the different codes for this particular example.

"Jose " <jose.l.vega@gmail.com> wrote in message <gqdr01$jp4$1@fred.mathworks.com>...

>
> I did not check the Bruno code, because He used the function v1*v2' and
> i am interested in (v1-v2)*inversa*(v1-v2)'.

Please for the next time make an minimum effort to post what you really wanted.

vdiff=reshape(bsxfun(@minus,...
              reshape(v1,[nrp 1 2]),...
              reshape(v2,[1 nrv 2])), ...
             [],2);
z = sum(reshape(sum((vdiff*inversa).*vdiff,2), ...
         [nrp nrv]),1);

% Bruno, who feels wasting his time

"Jose " <jose.l.vega@gmail.com> wrote in message <gqdtvl$n5j$1@fred.mathworks.com>...

>
> Ok roger, "apples and oranges" in different thread, but "apples with apples" in
> the test I showed above with the different codes for this particular example.

Common! I was told at first to grow oranges before knowing apple pies that will be made instead, and all that are in the same thread.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gqduat$hug$1@fred.mathworks.com>...
> "Jose " <jose.l.vega@gmail.com> wrote in message <gqdr01$jp4$1@fred.mathworks.com>...
>
> >
> > I did not check the Bruno code, because He used the function v1*v2' and
> > i am interested in (v1-v2)*inversa*(v1-v2)'.
>
> Please for the next time make an minimum effort to post what you really wanted.
>
> vdiff=reshape(bsxfun(@minus,...
> reshape(v1,[nrp 1 2]),...
> reshape(v2,[1 nrv 2])), ...
> [],2);
> z = sum(reshape(sum((vdiff*inversa).*vdiff,2), ...
> [nrp nrv]),1);
>
> % Bruno, who feels wasting his time


My apologize Bruno I made you waste your time,
sorry about that.

%piece of code.

nrv=3000; % random variable

nrp=5000; % number of 2D-points x(nrp,:)

v1=rand(nrp,2);
v2=rand(nrv,2);

inversa=[1 2; 3 4];

1. Roger

tic
 z = repmat(v1,nrv,1)-reshape(repmat(reshape(v2,1,[]),nrp,1),[],2);
 z = sum(reshape(sum(((z*inversa).*z),2),nrp,nrv));
toc


2. Bruno.

tic
vdiff=reshape(bsxfun(@minus,...
              reshape(v1,[nrp 1 2]),...
              reshape(v2,[1 nrv 2])), ...
             [],2);
z = sum(reshape(sum((vdiff*inversa).*vdiff,2), ...
         [nrp nrv]),1);
toc

Results:

Roger: Elapsed time is 3.702352 seconds.
Bruno: Elapsed time is 2.334867 seconds.

Thanks Bruno for your help.

Jose.

Yet another way:

bfun=@(x) sum((x*inversa).*x,2);
qfun=@(x,y) (x*(inversa+inversa.')*y.');
z=sum(bsxfun(@plus,bfun(v1),bfun(v2).')-qfun(v1,v2),1);

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gqe1nl$cra$1@fred.mathworks.com>...
> Yet another way:
>
> bfun=@(x) sum((x*inversa).*x,2);
> qfun=@(x,y) (x*(inversa+inversa.')*y.');
> z=sum(bsxfun(@plus,bfun(v1),bfun(v2).')-qfun(v1,v2),1);
>
> % Bruno



It is a very very clever code Bruno:

nrv=3000; % random variable

nrp=5000; % number of 2D-points x(nrp,:)

Elapsed time is 3993.921122 seconds. (With loops)
Elapsed time is 4.071653 seconds. (Roger)
Elapsed time is 2.480476 seconds. (You)
Elapsed time is 0.713359 seconds. (You: improved)

Really well done,

Thank you very much. I love this three!!! :)

Jose.

And what about this?

bfun=@(x) sum((x*inversa).*x,2);
z=(sum(bfun(v1))+nrp*bfun(v2) - ...
    (v2*((inversa+inversa.')*sum(v1,1).'))).';

Bruno

Just to show that the old For loop isn't as uncompetitive as it seems, try this out for giggles:


tic
z = zeros(1,nrv);
r = ones(nrp,1);
for jj = 1:nrv
    tmp = v2(jj,:);
    tmp2 = v1 - tmp(r,:);
    tmp3 = tmp2*inversa;
    z(jj) = sum(tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2));
end
toc

As is often the case, I noticed an improvement that could be made after I posted. On my machine this is as fast as any other. Though definitely not as cool as Bruno's final offer.


tic
z = zeros(1,nrv);
for jj = 1:nrv
    tmp2 = [v1(:,1)-v2(jj,1) v1(:,2)-v2(jj,2)];
    tmp3 = tmp2*inversa;
    zm(jj) = sum(tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2));
end
toc

"Matt Fig" <spamanon@yahoo.com> wrote in message <gqe7j5$evd$1@fred.mathworks.com>...
> As is often the case, I noticed an improvement that could be made after I posted. On my machine this is as fast as any other. Though definitely not as cool as Bruno's final offer.
>
>
> tic
> z = zeros(1,nrv);
> for jj = 1:nrv
> tmp2 = [v1(:,1)-v2(jj,1) v1(:,2)-v2(jj,2)];
> tmp3 = tmp2*inversa;
> zm(jj) = sum(tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2));
> end
> toc



Notice I didn't pre-allocate correctly (embarrassed z ~=zm)! With the correct pre-allocation, this is fastest of all. Good old For loops!

"Matt Fig" <spamanon@yahoo.com> wrote in message <gqe7ud$9vm$1@fred.mathworks.com>...
> "Matt Fig" <spamanon@yahoo.com> wrote in message <gqe7j5$evd$1@fred.mathworks.com>...
> > As is often the case, I noticed an improvement that could be made after I posted. On my machine this is as fast as any other. Though definitely not as cool as Bruno's final offer.
> >
> >
> > tic
> > z = zeros(1,nrv);
> > for jj = 1:nrv
> > tmp2 = [v1(:,1)-v2(jj,1) v1(:,2)-v2(jj,2)];
> > tmp3 = tmp2*inversa;
> > zm(jj) = sum(tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2));
> > end
> > toc
>
>
>
> Notice I didn't pre-allocate correctly (embarrassed z ~=zm)! With the correct pre-allocation, this is fastest of all. Good old For loops!


Hello guys,
in summary up to the moment the most faster code to work with big 2D-data points v1 and 2D-parameters v2 is done for Bruno (without loops), number 4.

(Look at the results).


%piece of code
clear all
clc


nrv=3000; % number of 2D-parameters

nrp=5000; % number of 2D-data points

v1=rand(nrp,2);
v2=rand(nrv,2);

inversa=[1 2; 3 4];

     
1. Roger

tic

 z = repmat(v1,nrv,1)-reshape(repmat(reshape(v2,1,[]),nrp,1),[],2);
 z = sum(reshape(sum(((z*inversa).*z),2),nrp,nrv));


toc

2. Bruno

tic
vdiff=reshape(bsxfun(@minus,...
              reshape(v1,[nrp 1 2]),...
              reshape(v2,[1 nrv 2])), ...
             [],2);
z = sum(reshape(sum((vdiff*inversa).*vdiff,2), ...
         [nrp nrv]),1);

toc


3. Bruno (improved)

tic
bfun=@(x) sum((x*inversa).*x,2);
qfun=@(x,y) (x*(inversa+inversa.')*y.');
z=sum(bsxfun(@plus,bfun(v1),bfun(v2).')-qfun(v1,v2),1);
toc

4. Bruno (Much better)

tic
bfun=@(x) sum((x*inversa).*x,2);
z=(sum(bfun(v1))+nrp*bfun(v2) - ...
    (v2*((inversa+inversa.')*sum(v1,1).'))).';

toc


5. Matt

tic
z = zeros(1,nrv);
r = ones(nrp,1);
for jj = 1:nrv
    tmp = v2(jj,:);
    tmp2 = v1 - tmp(r,:);
    tmp3 = tmp2*inversa;
    z(jj) = sum(tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2));
end
toc


6. Matt (improved)

tic
z = zeros(1,nrv);
for jj = 1:nrv
    tmp2 = [v1(:,1)-v2(jj,1) v1(:,2)-v2(jj,2)];
    tmp3 = tmp2*inversa;
    zm(jj) = sum(tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2));
end
toc


7. Matt (much better)

tic
z = zeros(1,nrv);
for jj = 1:nrv
tmp2 = [v1(:,1)-v2(jj,1) v1(:,2)-v2(jj,2)];
tmp3 = tmp2*inversa;
zm(jj) = sum(tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2));
end
toc

Results:

without loops:
1.Roger : Elapsed time is 4.267860 seconds.
2.Bruno :Elapsed time is 2.731296 seconds.
3.Bruno (improved): Elapsed time is 0.818409 seconds.
4.Bruno (much better) Elapsed time is 0.015496 seconds.

with lopps:

5.Matt :Elapsed time is 1.031870 seconds.
6.Matt (improve) Elapsed time is 0.462994 seconds.
7.Matt (much better)Elapsed time is 0.449142 seconds.

The code number 4 made for Bruno (without loops) is the most efficient code to built the z matrix with big 2D-data points v1 and 2D-parameters v2.

Althought the code (with loops) is not so bad result as I showed with my inefficient code called Jose, therefore Matt is right.

Jose.

> 7. Matt (much better)
>
> tic
> z = zeros(1,nrv);
> for jj = 1:nrv
> tmp2 = [v1(:,1)-v2(jj,1) v1(:,2)-v2(jj,2)];
> tmp3 = tmp2*inversa;
> zm(jj) = sum(tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2));
> end
> toc
>
> Results:
>
> without loops:
> 1.Roger : Elapsed time is 4.267860 seconds.
> 2.Bruno :Elapsed time is 2.731296 seconds.
> 3.Bruno (improved): Elapsed time is 0.818409 seconds.
> 4.Bruno (much better) Elapsed time is 0.015496 seconds.
>
> with lopps:
>
> 5.Matt :Elapsed time is 1.031870 seconds.
> 6.Matt (improve) Elapsed time is 0.462994 seconds.
> 7.Matt (much better)Elapsed time is 0.449142 seconds.
>






Jose, change the pre-allocation in the last loop code to zm instead of z, then try it. On my machine this is faster than any other.

Hats off to Bruno! I didn't see the shorter code he posted earlier. I don't think I can beat that one with a For loop. Nice ;)

"Matt Fig" <spamanon@yahoo.com> wrote in message <gqehad$k7c$1@fred.mathworks.com>...
>
>
> > 7. Matt (much better)
> >
> > tic
> > z = zeros(1,nrv);
> > for jj = 1:nrv
> > tmp2 = [v1(:,1)-v2(jj,1) v1(:,2)-v2(jj,2)];
> > tmp3 = tmp2*inversa;
> > zm(jj) = sum(tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2));
> > end
> > toc
> >
> > Results:
> >
> > without loops:
> > 1.Roger : Elapsed time is 4.267860 seconds.
> > 2.Bruno :Elapsed time is 2.731296 seconds.
> > 3.Bruno (improved): Elapsed time is 0.818409 seconds.
> > 4.Bruno (much better) Elapsed time is 0.015496 seconds.
> >
> > with lopps:
> >
> > 5.Matt :Elapsed time is 1.031870 seconds.
> > 6.Matt (improve) Elapsed time is 0.462994 seconds.
> > 7.Matt (much better)Elapsed time is 0.449142 seconds.
> >
>
>
>
>
>
>
> Jose, change the pre-allocation in the last loop code to zm instead of z, then try it. On my machine this is faster than any other.


Hello Matt, in my computer it does not change much:

7'. Matt code (With the right pre-allocation, zm instead of z):

tic
zm = zeros(1,nrv);
for jj = 1:nrv
tmp2 = [v1(:,1)-v2(jj,1) v1(:,2)-v2(jj,2)];
tmp3 = tmp2*inversa;
zm(jj) = sum(tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2));
end
toc

1.Elapsed time is 4.217884 seconds.
2.Elapsed time is 2.761333 seconds.
3.Elapsed time is 0.746436 seconds.
4.Elapsed time is 0.017458 seconds. (Bruno code)
5.Elapsed time is 1.037804 seconds.
6.Elapsed time is 0.477522 seconds.
7.Elapsed time is 0.455439 seconds.

7'. Matt code with the right pre-allocation (zm instead of z)
Elapsed time is 0.457492 seconds.

"Jose"
> us showed me how develop a matrix z
> avoiding loops....but....I felt really surprised how the algorithm with loops
> still run faster than the algorithm using a arrayfunction (without lopps).
> You can see it in the next results:
> nrv=2 % random variable
> nrp=3 % number of 2D-points points.
> v1=rand(nrp,2); % 2D-points
> v2=rand(nrv,2); % 2D-parameters
> [ix2,ix1]=meshgrid(1:nrv,1:nrp);
> r=arrayfun(@(x,y) v1(x,:)*v2(y,:)',ix1,ix2,'uni',false);
> m=sum(cell2mat(r)) % <- matrix

well, well, well...
it's a bit late to drop in, however
- you refer to an OP where you asked for a FULL BLOWN MATRIX -> R

http://www.mathworks.com/matlabcentral/newsreader/view_thread/247466

- NOW you really only want the sum of each col of R
     m=sum(cell2mat(r));
- this is an entirely different story and - believe me:
- no CSSMer in his/her right mind would use the above approach to achieve this...
- ...as is clearly demonstrated by senior CSSMers roger, matt (fig), and bruno
  by various, nice computational engines employing a col-based algorithm, which
  may use fast stock functions (eg, BSXFUN) and is less memory taxing, of course
- in summary:
  don't use a sledgehammer to fold a napkin...

http://www.napkinfoldingguide.com

just a thought...
us

I must agree with us here. If there is one thing to be learn from the thread is: Do not describe only half of the problem; All aspects of the problem are their importance when optimizing. If some characteristics are hidden, no effective optimization could be achieved.

Also it seems to me many complete/fast solutions have been proposed in other related threads (us & Roger), but they were overlooked for some reason. That's quite frustrating for people who spent their times to help, to finally see the same thread post over and over again.

Bruno



 

"Jose " <jose.l.vega@gmail.com> wrote in message <gqeh2t$aqc$1@fred.mathworks.com>...

>
> Althought the code (with loops) is not so bad result as I showed with my inefficient code called Jose, therefore Matt is right.
>

Under 100 points, for-loop might be the best choice. Things get worse at 1000. At 1e6 points, it even not possible to run it at reasonable time. But I'm sure Matt has not give his last word, and he might be able to modify his for-loop to overcome this.

But that's not so important. The point is: do not bother to vectorize for small problem. Some time readability and flexibility of for-loop is to be preferred.

You would not be happy to come back a year later and wonder how those obscure lines of highly optimized code work.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gqf9q6$gad$1@fred.mathworks.com>...
> "Jose " <jose.l.vega@gmail.com> wrote in message <gqeh2t$aqc$1@fred.mathworks.com>...
>
> >
> > Althought the code (with loops) is not so bad result as I showed with my inefficient code called Jose, therefore Matt is right.
> >
>
> Under 100 points, for-loop might be the best choice. Things get worse at 1000. At 1e6 points, it even not possible to run it at reasonable time. But I'm sure Matt has not give his last word, and he might be able to modify his for-loop to overcome this.
>
> But that's not so important. The point is: do not bother to vectorize for small problem. Some time readability and flexibility of for-loop is to be preferred.
>
> You would not be happy to come back a year later and wonder how those obscure lines of highly optimized code work.
>
> Bruno

Thank you guys for spent part of your time to solve my problem.
I learnt a lot about loops, stock functions and threads.
Next time I will ask you a specific question waiting for answers, and definitelly
will not generate different threads.
Thank Roger, Us, Matt and Bruno for your patience and your time.
Jose.

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gqf9q6$gad$1@fred.mathworks.com>...
> "Jose " <jose.l.vega@gmail.com> wrote in message <gqeh2t$aqc$1@fred.mathworks.com>...
>
> >
> > Althought the code (with loops) is not so bad result as I showed with my inefficient code called Jose, therefore Matt is right.
> >
>
> Under 100 points, for-loop might be the best choice. Things get worse at 1000. At 1e6 points, it even not possible to run it at reasonable time. But I'm sure Matt has not give his last word, and he might be able to modify his for-loop to overcome this.
>
> But that's not so important. The point is: do not bother to vectorize for small problem. Some time readability and flexibility of for-loop is to be preferred.
>
> You would not be happy to come back a year later and wonder how those obscure lines of highly optimized code work.
>
> Bruno


Hello everyone again, I am interested to solve my final problem using the most efficient algorithm, and I would like to use the stock function used by Bruno.

Problem:
In the next piece of code you can see the most efficient code in this thread with a for-loop suggested by Matt (Message 16) where I added a new vectors pm and lpm:


%piece of code

clear all
clc


nrv=3000; % random variable

nrp=5000; % number of 2D-points

v1=rand(nrp,2);
v2=rand(nrv,2);

inversa=[1 2; 3 4];

determ=2;

d=2;


%Matt code of Message 16 (with lopps) slightly modified:


tic
zm = zeros(1,nrv);
for jj = 1:nrv
    tmp2 = [v1(:,1)-v2(jj,1) v1(:,2)-v2(jj,2)];
    tmp3 = tmp2*inversa;
    zm = (tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2))'
% New vectors
    pm=(1/((2*pi)^(0.5*d)*(determ)^(0.5))*exp(-0.5*zm));
    lpm=log(1/((2*pi)^(0.5*d)*(determ)^(0.5)))-0.5*zm;
    lpms(jj)=sum(lpm);
end

toc

Question:

I am interested to keep the (lpm) vector, not only the sum of this vector, (lpms).
With the above modified efficient algorithm of Matt, I can obtain it very well,
but my question is:

Is possible to get it (lpm and lpms), using the stock function used by Bruno.

I would like to know if it is possible to modified the bruno code
write down in Message 13 to obtain the vector lpm and its sum (lpms), because are
two different things (thanks us for your explanation):

Bruno code from Message 13:

tic
 bfun=@(x) sum((x*inversa).*x,2);
 z=(sum(bfun(v1))+nrp*bfun(vt) - ...
     (vt*((inversa+inversa.')*sum(v1,1).'))).'
 
 toc

Thanks in advance

Jose.

"Jose " <jose.l.vega@gmail.com> wrote in message <gqg8go$o5j$1@fred.mathworks.com>...

>
> I would like to know if it is possible to modified the bruno code
> write down in Message 13 to obtain the vector lpm and its sum (lpms), because are
> two different things (thanks us for your explanation):
>

us was right. It will not be straightforward to adapt the code.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gqgg4e$98j$1@fred.mathworks.com>...
> "Jose " <jose.l.vega@gmail.com> wrote in message <gqg8go$o5j$1@fred.mathworks.com>...
>
> >
> > I would like to know if it is possible to modified the bruno code
> > write down in Message 13 to obtain the vector lpm and its sum (lpms), because are
> > two different things (thanks us for your explanation):
> >
>
> us was right. It will not be straightforward to adapt the code.
>
> Bruno

Thanks Bruno for your quickly answer!!

I was checking Roger's code as well to try to do it faster than Matt's code , and
the conclusion for my particular problem is that
one for-loop over (nrv) is faster than "repmatting" and "reshapping" avoiding the for-loops, even with big numbers.

cheers,

Jose.


Jose.

Why is pm being calculated when it is not being used? Why are you pre-allocating zm but not the array that is actually growing in the loop (lpms)?

"Matt Fig" <spamanon@yahoo.com> wrote in message <gqgj7p$dno$1@fred.mathworks.com>...
> Why is pm being calculated when it is not being used? Why are you pre-allocating zm but not the array that is actually growing in the loop (lpms)?

Good questions Matt.

I need pm, and you are right in the second question, I must pre-allocated lpms instead of zm.

Thanks,

Jose.

> I need pm
> Jose.

Since pm is being overwritten with each iteration, are you saying you only need the final pm (and lpm)? Or would you like to save all of the pm and lpm arrays calculated? If so they will have to be indexed too.

"Matt Fig" <spamanon@yahoo.com> wrote in message <gqgkiu$il3$1@fred.mathworks.com>...
> > I need pm
> > Jose.
>
> Since pm is being overwritten with each iteration, are you saying you only need the final pm (and lpm)? Or would you like to save all of the pm and lpm arrays calculated? If so they will have to be indexed too.

Hello Matt,
I am only interested in pm and lpm vectors for every iteration jj.
 I need these vectors to calculate a new vector (pm*lpm* (other vectors)).
After that, I do not need pm and lpm anymore and it will be overwriten in the
next iteration jj of the for-loop of your code.
For this reason they do not need to be indexed.

My code using your "for-loop code" run now really most faster than before.

Thanks for your interest Matt.

Cheers,
Jose.

"Jose " <jose.l.vega@gmail.com> wrote in message <gqgm9o$jif$1@fred.mathworks.com>...

>
> Hello Matt,
> I am only interested in pm and lpm vectors for every iteration jj.
> I need these vectors to calculate a new vector (pm*lpm* (other vectors)).
> After that, I do not need pm and lpm anymore and it will be overwriten in the
> next iteration jj of the for-loop of your code.

But you said you "I am interested to keep the (lpm) vector", then you erase them at every iteration???

PLEASE BE CLEAR WHAT YOU WANT TO PRESERVE AND WHAT YOU DON'T AT THE END.... OTHERWISE WE CAN'T OPTIMIZE THE CODE.

Bruno

"Bruno Luong"
> PLEASE BE CLEAR WHAT YOU WANT TO PRESERVE AND WHAT YOU DON'T AT THE END.... OTHERWISE WE CAN'T OPTIMIZE THE CODE...

(senior CSSMers now will certainly understand why the scolding tag was used in another thread using the same pattern...)

:-)
us

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gqgmkh$cj5$1@fred.mathworks.com>...
> "Jose " <jose.l.vega@gmail.com> wrote in message <gqgm9o$jif$1@fred.mathworks.com>...
>
> >
> > Hello Matt,
> > I am only interested in pm and lpm vectors for every iteration jj.
> > I need these vectors to calculate a new vector (pm*lpm* (other vectors)).
> > After that, I do not need pm and lpm anymore and it will be overwriten in the
> > next iteration jj of the for-loop of your code.
>
> But you said you "I am interested to keep the (lpm) vector", then you erase them at every iteration???
>
> PLEASE BE CLEAR WHAT YOU WANT TO PRESERVE AND WHAT YOU DON'T AT THE END.... OTHERWISE WE CAN'T OPTIMIZE THE CODE.
>
> Bruno


Bruno, yes you are right, but sometime is difficutl to explain it.
Keep lpm...I mean, I need this vector and it can't be obtained in your code,
because again I confused you asking for the sum of the vector, and not the vector,
two different things, well clarified for "us".

Now, With the Matt code I obtain "lpm" explicitelly
and I can use it in every iteration jj, as well as the sum of this vector called "lpms" growing in the only one "for-loop".

Sorry, I am not English, and "keep lpm" was the wrong expression to explain it.

cheers,
Jose.

Well, you can shave off about 10% of the run-time by using the following:



nrv = 3000; % random variable
nrp = 5000; % number of 2D-points
v1 = rand(nrp,2);
v2 = rand(nrv,2);
ia = [1 2; 3 4];
determ = 2;
d = 2;
% There is no reason to keep calculating these guys in the loop!
fact1 = 1/((2*pi)^(0.5*d)*(determ)^(0.5));
fact2 = log(fact1);

tic
lpms = zeros(1,nrv);
for jj = 1:nrv
    t1 = v1(:,1) - v2(jj,1);
    t2 = v1(:,2) - v2(jj,2);
    zm = -0.5*(t1.^2*ia(1) + (ia(2)+ia(3))*t1.*t2 + ia(4)*t2.^2);
    pm = fact1 * exp(zm);
    lpm = fact2 + zm;
    lpms(jj) = sum(lpm);
end
toc


The profiler says that the exp function takes 70% of run-time. This being the case, I would be surprised if there was a method that was WAY faster out there. Of course Bruno (and others!) has surprised me before, but there have been several threads in the past dealing with trying to speed up exp. This being the case, short of a faster approximating scheme to exp, I think this may be close to bottomed out.

Good luck.

"Matt Fig" <spamanon@yahoo.com> wrote in message <gqhk3l$4eb$1@fred.mathworks.com>...
> Well, you can shave off about 10% of the run-time by using the following:
>
>
>
> nrv = 3000; % random variable
> nrp = 5000; % number of 2D-points
> v1 = rand(nrp,2);
> v2 = rand(nrv,2);
> ia = [1 2; 3 4];
> determ = 2;
> d = 2;
> % There is no reason to keep calculating these guys in the loop!
> fact1 = 1/((2*pi)^(0.5*d)*(determ)^(0.5));
> fact2 = log(fact1);
>
> tic
> lpms = zeros(1,nrv);
> for jj = 1:nrv
> t1 = v1(:,1) - v2(jj,1);
> t2 = v1(:,2) - v2(jj,2);
> zm = -0.5*(t1.^2*ia(1) + (ia(2)+ia(3))*t1.*t2 + ia(4)*t2.^2);
> pm = fact1 * exp(zm);
> lpm = fact2 + zm;
> lpms(jj) = sum(lpm);
> end
> toc
>
>
> The profiler says that the exp function takes 70% of run-time. This being the case, I would be surprised if there was a method that was WAY faster out there. Of course Bruno (and others!) has surprised me before, but there have been several threads in the past dealing with trying to speed up exp. This being the case, short of a faster approximating scheme to exp, I think this may be close to bottomed out.
>
> Good luck.


Thanks Matt for your help,
You are right I have compared the Roger Code, Your code, and Your
last code Modified, and I can shave off about 10% of the run-time
from your first one.


% code


nrv=2000; % number of 2D-variable

nrp=3000; % number of 2D-points

v1=rand(nrp,2);
v2=rand(nrv,2);

ia=[1 2; 3 4];

determ=2;

d=2;


fact1 = 1/((2*pi)^(0.5*d)*(determ)^(0.5));
fact2 = log(fact1);

1. Roger

tic
  zx = repmat(v1,nrv,1)-reshape(repmat(reshape(v2,1,[]),nrp,1),[],2);
  y = (reshape(sum(((zx*ia).*zx),2),nrp,nrv));
  p1=fact1*exp(-0.5*y);
  lp1=fact2-0.5*y;
  zx2=sum(lp1);
 toc

1.Elapsed time is 2.093532 seconds.



2. Matt Code.

tic
lpms = zeros(1,nrv);
for jj = 1:nrv
    tmp2 = [v1(:,1)-v2(jj,1) v1(:,2)-v2(jj,2)];
    tmp3 = tmp2*ia;
    zm = (tmp3(:,1).*tmp2(:,1) + tmp3(:,2).*tmp2(:,2))';
    pm=fact1*exp(-0.5*zm);
    lpm=fact2-0.5*zm;
    lpms(jj)=sum(lpm);
end

toc


2.Elapsed time is 0.641402 seconds.


3. Matt code (Modified) shaving of 10 % of run-time.


tic
lpms = zeros(1,nrv);
for jj = 1:nrv
    t1 = v1(:,1) - v2(jj,1);
    t2 = v1(:,2) - v2(jj,2);
    zm = -0.5*(t1.^2*ia(1) + (ia(2)+ia(3))*t1.*t2 + ia(4)*t2.^2);
    pm = fact1 * exp(zm);
    lpm = fact2 + zm;
    lpms(jj) = sum(lpm);
end
toc

3.Elapsed time is 0.540709 seconds.

It is a good point Matt for working out over the exponential function.
Thanks again.

Many thanks guys for your help.
 I learnt a lot from you in this thread,
my apologize again for confusing you with my english.

Cheers,

Jose.