can this loop be vectorized?
3 views (last 30 days)
Show older comments
Here's a simplified example of some data I have. I'd like to eliminate the loop because my x and y arrays in reality have tens of thousands of elements, not just three.
X = repmat(0:5:1000,201,1);
Y = repmat((1000:-5:0)',1,201);
x = [50 326 800];
y = [900 429 600];
H = NaN(201,201,length(x));
for k = 1:length(x)
H(:,:,k) = hypot(X-x(k),Y-y(k));
end
0 Comments
Accepted Answer
Matt J
on 30 May 2014
Edited: Matt J
on 30 May 2014
Doing bsxfun(@minus,X,...) will be inefficient since X contains a lot of replicated data. So, I propose,
X0=(0:5:1000).'; sX=length(X0);
Y0=(1000:-5:0).'; sY=length(Y0);
[I,J]=ndgrid(1:sY,1:sX);
Xc=bsxfun(@minus,X0,x(:).').^2;
Yc=bsxfun(@minus,Y0,y(:).').^2;
H=reshape( sqrt( Yc(I,:)+Xc(J,:)) ,sY,sX,[]);
5 Comments
More Answers (3)
José-Luis
on 30 May 2014
Edited: José-Luis
on 30 May 2014
bsxfun() is your friend:
X = repmat(0:5:1000,201,1);
Y = repmat((1000:-5:0)',1,201);
x(1,1,1:3) = [50 326 800];
y(1,1,1:3) = [900 429 600];
tic
H = NaN(201,201,length(x));
for k = 1:length(x)
H(:,:,k) = hypot(X-x(k),Y-y(k));
end
toc
tic
H_new = hypot(bsxfun(@minus,X,x),bsxfun(@minus,Y,y));
toc
Vectorized code is not necessarily faster. If performance is really an issue, you could try a mex file.
EDIT
Another option:
tic
H_nuevo = hypot(repmat(X,[1 1 3]) - repmat(x,[201 201 1]) , ...
repmat(Y,[1 1 3]) - repmat(y,[201 201 1]) );
toc
An idea to increase the efficiency of your code is to try to compare squared distances instead, since the sqrt() is relatively expensive to perform.
6 Comments
José-Luis
on 2 Jun 2014
Thanks Sean. Good to know. I find it hard to know a-priori whether Matlab is going to perform decently or will be mired in overheads. Any good rule of thumb?
Sean de Wolski
on 2 Jun 2014
What we say at seminars is:
- You'll always see a speedup for signal processing, communications, or other applications where there is a persistent state that needs to be updated
- Anything with fixed-point
You won't see a speedup for:
- Linear Algebra, ffts, or functions that use IPP.
But, your milage will vary :)
Sean de Wolski
on 30 May 2014
I'm seeing a slight (1.25ish times) speedup with the bsxfun approach
function speedyhypot()
X = repmat(0:5:1000,201,1);
Y = repmat((1000:-5:0)',1,201);
x = [50 326 800];
y = [900 429 600];
t1 = timeit(@FirstWay);
t2 = timeit(@SecondWay);
disp([t1 t2])
function FirstWay
H = NaN(201,201,length(x));
for k = 1:length(x)
H(:,:,k) = hypot(X-x(k),Y-y(k));
end
end
function SecondWay
H2 = hypot(bsxfun(@minus,X,reshape(x,1,1,3)),bsxfun(@minus,Y,reshape(y,1,1,3)));
end
end
Jan
on 31 May 2014
Edited: Jan
on 31 May 2014
hypot is smart with avoiding overflows. But for limited values sqrt(a^2+b^2) is much faster:
for k = 1:length(x)
H(:,:,k) = sqrt((X-x(k)).^2 + (Y-y(k)).^2);
end
I get a speedup of factor 3, much more than the bsxfun tricks.
The power operator is expensive, so it is worth to try to reduce the squaring by using the identity: (a-b)^2 == a^2 - 2*a*b + b^2:
X2 = X .^ 2;
Y2 = Y .^ 2;
for k = 1:length(x)
H(:,:,k) = sqrt(X2 - X*(2*x(k)) + x(k)^2 + Y2 - Y*(2*y(k)) + y(k)^2);
end
2*X*x(k) is less efficient than X*(2*x(k)), because in the first case you have two multiplications of a scalar with a matrix, but in the second one one matrix operation only.
But this is slower than the first method. It seems like Matlab smartly performs a multiplikation for .^2 and not a power operation.
5 Comments
See Also
Categories
Find more on Logical in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!