Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
vectorization question

Subject: vectorization question

From: Peter Schreiber

Date: 13 Jun, 2011 09:40:05

Message: 1 of 4

Hello,
Is there any chance to get rid of this for loop by vectorization?
Is it possible to do this by utilizing three-dim. arrays ?

Thanks in advance,
Peter

for counter = 1:max_count
A=(P5(:,counter)-P4(:,counter))/norm(P5(:,counter)-P4(:,counter));
As=repmat(A',n_size_of ,1) - repmat(2*sum(N.*(repmat(A',n_size_of ,1)),2),1,3).*N;
end


where
                                             
max_count = 10

size(A)
ans = 3 1


size(As)
ans = 5000 3


size(P5)
ans =
           3 1000

size(P4)
ans =
           3 1000

size(N)
ans =
        5000 3

n_size_of = 5000

Subject: vectorization question

From: Matt J

Date: 13 Jun, 2011 13:38:04

Message: 2 of 4

"Peter Schreiber" <schreiber.peter15@gmail.com> wrote in message <it4ltl$n6p$1@newscl01ah.mathworks.com>...
> Hello,
> Is there any chance to get rid of this for loop by vectorization?
> Is it possible to do this by utilizing three-dim. arrays ?
>
> Thanks in advance,
> Peter
>
> for counter = 1:max_count
> A=(P5(:,counter)-P4(:,counter))/norm(P5(:,counter)-P4(:,counter));
> As=repmat(A',n_size_of ,1) - repmat(2*sum(N.*(repmat(A',n_size_of ,1)),2),1,3).*N;
> end
>

You have a typo somewhere. Every iteration of the loop overwrites the previous one, so there is little here that makes sense.

Also, if max_count is really so small, there seems to be little reason to get rid of the loop. Your main overhead comes from the use of all the REPMATs in the computation of As. You can get rid of those as follows

 v= N*(2*A);
 As=bsxfun(@times,v,N);
 As=bsxfun(@minus,A',As);

In fact though, it will probably be more efficient in this case to use an additional for-loop instead of BSXFUN since you only have to loop from 1 to 3:

v= N*(2*A);
As=zeros(size(N));
for i=1:3
  As(:,i)=A(i)- N(:,i)*v(i);
end

Subject: vectorization question

From: Matt J

Date: 13 Jun, 2011 13:55:04

Message: 3 of 4

"Matt J" wrote in message <it53rs$q5$1@newscl01ah.mathworks.com>...
>
 In fact though, it will probably be more efficient in this case to use an additional for-loop instead of BSXFUN since you only have to loop from 1 to 3:
>
> v= N*(2*A);
> As=zeros(size(N));
> for i=1:3
> As(:,i)=A(i)- N(:,i)*v(i);
> end
==================

Sorry, that should be

v= N*(2*A);
As=N;
for i=1:3
  As(:,i)=A(i)- N(:,i).*v;
end
 

Subject: vectorization question

From: Roger Stafford

Date: 13 Jun, 2011 16:18:02

Message: 4 of 4

"Peter Schreiber" <schreiber.peter15@gmail.com> wrote in message <it4ltl$n6p$1@newscl01ah.mathworks.com>...
> for counter = 1:max_count
> A=(P5(:,counter)-P4(:,counter))/norm(P5(:,counter)-P4(:,counter));
> As=repmat(A',n_size_of ,1) - repmat(2*sum(N.*(repmat(A',n_size_of ,1)),2),1,3).*N;
> end
- - - - - - - -
  My guess is you are doing ray tracing off reflective surfaces and N is a set of 5000 unit normals to surface points. I guess further that you would actually like to have a max_count of 1000 instead of 10 since that is the size of P4 and P5 rows.

  Your computation can be vectorized to produce a 3 x 5000 x 1000 3D array of (ray?) vectors using 'bsxfun' (similar to Matt's code) provided you have enough memory to store it. Otherwise you would have to cut down the 1000 to lower values, as you have done. Also I am not sure this would be any faster than something involving for-loops, but you can test that yourself.

 A = P5-P4; % 3 x 1000
 A = bsxfun(@rdivide,A,sqrt(sum(A.^2,1))); % 3 x 1000)
 As = bsxfun(@times,reshape(N*(2*A),1,5000,1000),N.'); % 3 x 5000 x 1000
 As = bsxfun(@minus,reshape(A,3,1,1000),As) % 3 x 5000 x 1000

Roger Stafford

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us