Your a
matrix is a 1D vector and is incompatible with the nested loop, which computes distance in 2D space from each point to each other point. So the following answer applies to the problem of finding all pairwise distances in a N-by-D
matrix, as your loop does for the case of D=2
.
Option 1 - pdist
I think you are looking for pdist
with the 'euclidean'
distance option.
a = randn(10, 2); %// 2D, 10 samples
D = pdist(a,'euclidean'); %// euclidean distance
Follow that by squareform
to get the square matrix with zero on the diagonal as you want it:
distances = squareform(D);
Option 2 - bsxfun
If you don't have pdist
, which is in the Statistics Toolbox, you can do this easily with bsxfun
:
da = bsxfun(@minus,a,permute(a,[3 2 1]));
distances = squeeze(sqrt(sum(da.^2,2)));
Option 3 - reformulated equation
You can also use an alternate form of Euclidean (2-norm) distance,
||A-B|| = sqrt ( ||A||^2 + ||B||^2 - 2*A.B )
Writing this in MATLAB for two data arrays u
and v
of size NxD
,
dot(u-v,u-v,2) == dot(u,u,2) + dot(v,v,2) - 2*dot(u,v,2) % useful identity
%// there are actually small differences from floating point precision, but...
abs(dot(u-v,u-v,2) - (dot(u,u,2) + dot(v,v,2) - 2*dot(u,v,2))) < 1e-15
With the reformulated equation, the solution becomes:
aa = a*a';
a2 = sum(a.*a,2); % diag(aa)
a2 = bsxfun(@plus,a2,a2');
distances = sqrt(a2 - 2*aa);
You might use this method if Option 2 eats up too much memory.
Timings
For a random data matrix of size 1e3-by-3 (N-by-D), here are timings for 100 runs (Core 2 Quad, 4GB DDR2, R2013a).
- Option 1 (
pdist
): 1.561150 sec (0.560947 sec in pdist
)
- Option 2 (
bsxfun
): 2.695059 sec
- Option 3 (
bsxfun
alt): 1.334880 sec
Findings: (i) Do computations with bsxfun
, use the alternate formula. (ii) the pdist
+squareform
option has comparable performance. (iii) The reason why squareform
takes twice as much time as pdist
is probably because pdist
only computes the triangular matrix since the distance matrix is symmetric. If you can do without the square matrix, then you can avoid squareform
and do your computations in about 40% of the time required to do it manually with bsxfun
(0.5609/1.3348).