How to vectorize a for loop?

6 次查看（过去 30 天）
Tony
Tony 2012-9-5
How can i make that code run faster?
data_file=importdata('Hz.txt');
xA=data_file(:,1);
yA=data_file(:,2);
zA=data_file(:,3);
[ndat,yyy]=size(xA); #ndat=540
aA=1./(xA+1); %Number of randoms Ran=10000;
X1=randi([-100,150],[1,Ran])/100;
s1=randi([-100,150],[1,Ran])/100;
s2=randi([-100,150],[1,Ran])/100;
sw=randi([50,130],[1,Ran]);
T=zeros([ndat,Ran]);
H=zeros([ndat,Ran]);
XsA=zeros([1,ndat]);
Chit=zeros([1,Ran]);
t=0:0.0000000001:0.04;
for j=1:Ran
x1=X1(j);
y1=s1(j);
y2=s2(j);
w=(sw(j));
x2=-x1-1/4*y1^2-3/2*y1*y2-1/4*y2^2;
As=x1*y1*y2*exp(-w*t)+x2;
for i=1:ndat
dA=abs(As-aA(i));
[v T(i,j)]=min(dA);
I=T(i,j);
H(i,j)=x1*y2+exp(-w*t(I))+exp(x2);
XsA(i)=(H(i,j)-yA(i)).^2 ./(zA(i).^2);
Chit(j)=sum(XsA);
end
end
end
  1 个评论
Jan
Jan 2012-9-6
编辑：Jan 2012-9-6
An optimization needs a look in the profiler (although this disables the important JIT!) and the possibility to run the code. Even advanced programmers cannot predict, how the JIT will handle modifications of the code. Therefore it would be helpful, if you provide test data. E.g. RAND with a fixed seed is a good and cheap method in the forum - if such data are valid for the tests.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2012-9-6
编辑：Jan 2012-9-6
At first follow the standard rule to avoid all repeated calculations by moving them befor the loop:
c1 = exp(-w * t);
c3 = 1 ./ (zA .^ 2);
for j=1:Ran
x1 = X1(j);
y1 = s1(j);
y2 = s2(j);
w = sw(j);
x2 = -x1 - 0.25 * y1^2 - 1.5 * y1 * y2 - 0.25 * y2^2;
As = x1 * y1 * y2 * c1 + x2;
c4 = x1 * y2 + exp(x2);
for i = 1:ndat
dA = abs(As - aA(i));
[v, I] = min(dA);
T(i,j) = I;
H(i,j) = c4 + c1(I);
XsA(i) = (H(i,j) - yA(i)) .^ 2 .* c3(i);
end
Chit(j) = sum(XsA); % I guess this should be here
end
From this point a vectorizeation would mean to replace "i" inside the loop by "1:ndat".
  1 个评论
Tony
Tony 2012-9-6
the Chit is ok. and also c4 and c3. But w,t are vectors so c1 is not working. Also i have the same code for another problem where inside the 1:ndat loop there is a quad command. i.e. H=quad(As,0,t(I))*aA(i).

请先登录，再进行评论。

更多回答（2 个）

Doug Hull
Doug Hull 2012-9-6
400 Million long vector t seems excessive to me.
What if you make it significantly shorter? I would lower your number of iterations on everything. See if that gives good results. Then start cranking up the number of iterations, and resolution. If it is not fine enough, then run it through the profiler to see where the bottlenecks are. My guess is that doing the above will get the results you need without code changes.
  1 个评论
Tony
Tony 2012-9-6
When i define a smaller vector t , i don't have accuracy in the dA. But if that cost, i think that will be ok if i can set 4M long t.

请先登录，再进行评论。

Robert Cumming
Robert Cumming 2012-9-6
you over write:
Chit(j)=sum(XsA);
in every i loop - are you sure you want to do that?
Profile the code and you will see which lines are the most expensive.
  2 个评论
Tony
Tony 2012-9-6
I have correct that.
Jan
Jan 2012-9-6
编辑：Jan 2012-9-6
Where have you corrected this? Is sum(XsA) calculated inside or outside the for i loop?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by