# MATLAB Puzzler: Finding the two closest points 23

Posted by **Doug Hull**,

```
n = 10;
x = rand(1,n);
y = rand(1,n);
```

Brett and I took very different paths to solving this problem, he knows more about the statistics toolbox than I do and found a much shorter solution because there is a function that basically does this. As a benchmark. sticking with base MATLAB, I solved this in seven lines of code. Post your solutions in the comments area along with how long it took you to solve.

## 23 CommentsOldest to Newest

**1**of 23

```
n = 10;
x = rand(1,n);
y = rand(1,n);
min_dist = Inf;
min_p1 = 0;
min_p2 = 0;
for p1 = 1:n
for p2 = 2:n
d = (x(p1)-x(p2))^2+(y(p1)-y(p2))^2;
if d < min_dist && p1 ~= p2
min_p1 = p1;
min_p2 = p2;
min_dist = d;
end
end
end
plot(x,y,'b.',[x(min_p1),x(min_p2)],[y(min_p1),y(min_p2)],'ro');
```

I like your solution for putting more of Matlab's unique features to use. For large N, I have a suspicion that the matrix based approach will become a bit of a memory issue, but I it'll be much faster for smaller N, taking advantage of the fast matrix/vector operations in Matlab.
The code above is the trivial solution of forming all combinations of points in sequence, measuring the distance and keeping track of the minimum (eliminating the case where both chosen points are the same point). **2**of 23

**3**of 23

**4**of 23

```
min_dist = Inf;
for p1 = 1:(n-1)
for p2 = (p1+1):n
d = (x(p1)-x(p2))^2+(y(p1)-y(p2))^2;
if d < min_dist
min_p1 = p1;
min_p2 = p2;
min_dist = d;
end
end
end
```

@ Rob
Good thinking to use the built in combination generator.
@Doug
Keep putting up these problems, they're fun and short enough that I can still get some real work done. **5**of 23

```
dist = squareform(pdist([x',y']));
```

WITHOUT the toolbox, check this out:
```
[xs,ys] = meshgrid(x,y);
dists = sqrt((xs-xs').^2 + (ys-ys').^2);
```

Cheers!
Brett **6**of 23

```
[ii,jj] = find(dist == min(dist(dist~=0)));
```

(or DISTS for the second case.)
Brett again **7**of 23

`[ix,iy] = find(dist==min(dist));`

This is an operation of O(2n), first find the min(O(n)) and then find the index if the min (O(n)). For large N, the following might be faster (O(n)), although of course it requires 2 lines due to matlab syntax.
```
[minVal, minIndex] = min(dist);
[ix,iy] = ind2sub(size(dist),minIndex);
```

--DA **8**of 23

```
%% puzzler
n = 10;
x = rand(1,n);
y = rand(1,n);
%% solution
[A,B] = meshgrid(2:n,1:n-1); %pairwise indices - all combinations
D = (x(A)-x(B)).^2 + (y(A)-y(B)).^2; %distance squared between each pair
D(tril(ones(n-1),-1)==1) = inf; %all values below diagonal => inf
[row,col] = find(D==min(D(:))); %locate smallest distance
closest = [A(row,col) B(row,col)]; %indices of closest 2 points
%% display results
plot(x,y,'b.',x(closest),y(closest),'ro')
axis equal
```

Then I looked at other solutions. Very interesting ideas. I learned a few things this morning. :)
Cheers
Bob **9**of 23

```
% Given:
n = 10;
x = rand(1,n);
y = rand(1,n);
% Solution:
[a,b] = meshgrid(1:n,1:n);
dmat_sq = (x(a) - x(b)).^2 + (y(a) - y(b)).^2;
dmat_sq(~dmat_sq) = Inf;
[I,J] = find(dmat_sq == min(dmat_sq(:)));
plot(x,y,'.',x(I),y(I),'o')
```

Five lines. Five minutes. After I wrote it, I looked at Doug's solution, and I think I took the long way around with my use of the meshgrid function (generating indices rather than x,y matrices). Also, mine would fail if two of the x,y points were identical... **10**of 23

```
% Find closest pair of points
% starter code
n=10;
x=rand(n,1);
y=rand(n,1);
% go to complex vector representation
z = x + i*y;
% find distances between all pairs
a = kron(z, ones(1, length(z)));
b = kron(z.', ones(size(z)));
zdist = abs(a-b);
% only need to look in lower triangle for symmetric pair-wise dists
mindist = min(zdist(tril(zdist)~=0));
%row and col; will contain indexes for the pair, but only need to keep one
[row, col] = find(zdist==mindist);
%plot
figure(1);
scatter(real(z), imag(z), '.b');
hold on;
scatter(real(z(row)), imag(z(row)), 'or');
hold off;
```

% I always forget to use the meshgrid function unless I'm doing something 3-d, but some people might want to know about kron products anyway.
% Tril (and triu) can be convenient filters.
% And for 2-d, complex variables work ok, but going into 3+ dimensions, then I imagine HYPOT or it's higher-dim cousins would keep working. **11**of 23

`n = 10; x = rand(1,n); y = rand(1,n); % Find the closest points [xx,yy]=meshgrid(x,y); %(1) zz = triu((xx-xx.').^2 + (yy-yy.').^2); %(2) zz(zz==0)=inf; %(3) [d,i]=min(zz(:)); %(4) should use ind2sub, but requires 2 outputs disp(['The closes points are: ' num2str([mod(i,10)+(mod(i,10)==0), ceil(i/10)])]) %(5) %% Plotting code ij = [mod(i,10)+(mod(i,10)==0), ceil(i/10)]; plot(x,y,'*b', x(ij),y(ij),'or-') text(x(ij(1)),y(ij(1))-.05,num2str(d)) axis equal`

I think it takes a full 5 lines to get the solution. I took about 15 minutes and had to add the graphing code to make sure I was getting it right. I definitely wanted a shortcut syntax to get the index of the minimum value without having to get both outputs... Thanks for the distraction DOUG!

**12**of 23

```
dd=(x(ones(n,1),:)-x(ones(n,1),:)').^2+(y(ones(n,1),:)-y(ones(n,1),:)').^2;
dd(1+(0:(n-1))*(n+1))=NaN;
[vmin,imin]=min(dd(:)); % uniqueness
[xi,yi]=ind2sub([n,n],imin);
```

two comments
1) setting the diagonal elements directly - more efficient memory-wise. for readability may use
sub2ind([n,n],1:n,1:n)
instead of
1+(0:(n-1))*(n+1)
2) by using the min function with two outputs, get a unique answer (instead of at least two when using find) **13**of 23

```
n=10;
x=rand(n,1);
y=rand(n,1);
```

This version uses complex numbers to simplify the expressions. I think it's very readable, too (except for the odd "nan" thingy which tells "min" to ignore the diagonal). I came up with it in about 5 minutes, without reading any of the above first:
```
[X,Y] = meshgrid (x + 1i*y) ;
D = abs (X-Y) + diag (nan (n,1)) ;
[d k] = min (D (:)) ;
i = mod (k-1,n) + 1 ;
j = fix ((k-1)/n) + 1 ;
```

**14**of 23

```
code goes here.
```

</pre></code>
[edit by Doug to have both <pre> and <code> tags in suggestion. Also edited old comment to do this. **15**of 23

**16**of 23

```
dist=sqrt(bsxfun(@minus,x,x').^2+bsxfun(@minus,y,y').^2);
dist(dist==0)=10; % Cancel the main diag
[i,j]=find(dist==min(dist(:)))
plot(x,y,'.',x(i),y(i),'o')
```

**17**of 23

**18**of 23

**19**of 23

```
n = 10;
x = rand(1,n);
y = rand(1,n);
allDist = pdist([x;y]');
[minDist,minLoc] = min(allDist);
[loc1,loc2] = find(tril(ones(n),-1));
minLoc1 = loc1(minLoc);
minLoc2 = loc2(minLoc);
```

Alternatively, we can use bsxfun to build up the distance map.
```
z = x + i*y;
distMap = tril(abs(bsxfun(@minus,z,z.')));
distMap(~distMap) = Inf;
[minCol,rowLoc] = min(distMap);
[minDist,colLoc] = min(minCol);
rowLoc = rowLoc(colLoc);
```

**20**of 23

`[i,j]=find(squareform(pdist([x; y]'))==min(pdist([x; y]')),1,'first')`

**21**of 23

```
d = pdist([x;y]');
[i,j] = find(squareform(d)==min(d),1,'first');
```

The length of code can be checked as follows:
`t = mtree('[i,j]=find(squareform(pdist([x; y]''))==min(pdist([x; y]'')),1,''first'')')`

which gives 33 nodes.
`t = mtree('d=pdist([x; y]'');[i,j]=find(squareform(d)==min(d),1,''first'')')`

only 28 nodes.
hth
Yi **22**of 23

```
flagIsBig = (variable > 1000);
if flagIsBig
%Do big code
else %not flagIsBig
%Do not big code
end
```

This takes more characters, but is very readable. I think only community voting can clear this up. What do you think?
-Doug **23**of 23

## Recent Comments