# Closest interior point on integer grid to a vertex of a convex polyhedron

I have a 3 dimensional convex polyhedron whose vertex coordinates are rational. For one of these vertices, I would like to find the nearest integer grid point (under the Euclidean metric) that is contained within the polyhedron. The polyhedron can be specified by a collection of inequalities, or as a collection of vertices with appropriate combinatorial information.

I assume I'm looking for some sort of integer program, but I can't quite match up what I want with what I'm reading. How do you solve this problem, and is there efficient code (preferably C/C++) available to do it?

Thanks!

In the comments to Johan's post I said it seems a shame to throw a full MIQP solver at this. For a general $n$-dimensional polyhedron, I'd certainly hold to that. But since this is a 3-dimensional problem, it might be competitive to do an intelligent exhaustive search. I suppose it depends on the application.

First, suppose we have constructed a generator that enumerates all of the lattice points in (partial) order of increasing Manhattan distance from the target vertex $x_k$. That is, if $k>l$, then $\|x_k-z\|_\infty \geq \|x_l-z\|_\infty$. Mahnattan distance is easier for enumeration purposes; we'll handle the need to pick the best Euclidean distance in the filtering algorithm below.

Now for the filter: let $d_{\text{max}}=\max_i\|v_i-z\|_\infty$, where $v_i$ are the vertices of the polygon, and set $d_{\text{best}}=+\infty$. Then, for $k=1,2,3,...$:

1. If $\|x_k-z\|_\infty > d_{\text{max}}$, terminate.
2. If $x_k$ falls outside of the polygon, repeat steps 1-5 with $k\leftarrow k+1$.
3. If $\|x_k-z\|_\infty > d_{\text{best}}$, terminate.
4. If $\|x_k-z\|_2 < d_{\text{best}}$, set $x_{\text{best}}\leftarrow x_k$ and $d_{\text{best}}=\|x_k-z\|_2$.
5. Repeat steps 1-5 with $k\leftarrow k+1$.

This loop will terminate in step 1 if there are no lattice points inside the polygon, and in either step 1 or step 3 when no points closer than $d_{\text{best}}$ remain. If $d_{\text{best}}=+\infty$ on termination, no lattice point was found within the polyhedron.

Understanding step 3 is important: since $\|x-z\|_2 \geq \|x-z\|_\infty$, this guarantees that no subsequent points from the generator will have a lower Euclidean distance.

So construct a generator and we're all set, right? Well, guaranteeing the partial ordering $k\geq l ~\rightarrow~ \|x_k-z\|_\infty \geq \|x_l-z\|_\infty$ might be difficult. In particular, it will likely require an internal buffer for sorting, which would be worth avoiding. So let's go simpler: let $\tilde{z}$ be the lattice point obtained by rounding the coordinates of $z$, and consider a a generator that guarantees $\|x_k-\tilde{z}\|_\infty \geq \|x_l-\tilde{z}\|_\infty$ for $k\geq l$. This is much simpler: for a 3-dimensional lattice, this can be programmed with 3 or 4 nested for loops.

The challenge now is that you need to relax Step 3 a bit, because the points coming out of the generator no longer satisfy the original partial ordering. But $$\|x_k-z\|_\infty = \|x_k-\tilde{z}+(z-\tilde{z})\|_\infty \leq \|x_k-\tilde{z}\|_\infty + \|z-\tilde{z} \|_\infty \leq \|x_k-\tilde{z}\|_\infty + 0.5$$ So I believe that relaxing the test in Step 3 to $\|x_k-\tilde{z}\|_\infty + 0.5 > d_{\text{best}}$ will restore correctness. I'd change $d_{\text{max}}=\max_i\|v_i-\tilde{z}\|_\infty$ and modify Step 1 to test $\|x_k-\tilde{z}\|_\infty$ as well.

You're going to be testing $O(\lceil d_{\text{best}}\rceil^3)$ points with this approach, with a worst case of $O(\lceil d_{\text{max}} \rceil^3)$. But depending upon your application, it might be better than using a general-purpose MIQP solver. And it will be easier than coding one of your own for the same purpose.

If you simply want a formulation (not necessarily a good way to solve the problem), you can state it as minimizing the distance $||x-z||$ where $z$ is an integer vector, satisfying $Az\leq b$ (inside polyhedron) and $x$ is a binary combination of the vertices $x = \sum_{i=1}^N \delta_i v_i$ where $\sum_{i=1}^N \delta_i = 1$, $\delta$ binary, and $v_i$ are the vertices of the polytope (if you only have the vertices of the polytope, you can replace the constraint $Az\leq b$ by $z = \sum_{i=1}^N \lambda_i v_i$ where $\sum_{i=1}^N \lambda_i = 1$, $\lambda \geq 0$.) Depending on the distance measure you use, both models generate mixed-integer linear or quadratic program.

The following is an implementation in the MATLAB Toolbox YALMIP (developed by me). It assumes you have a reasonably good MIQP solver installed (although the problem seems to be trivially simple for the solvers I tried)

% Data
V = randn(2,50)*2;
k = convhull(V');
V = V(:,k);
clf
plot(V(1,:),V(2,:),'*');
N = size(V,2);

delta  = binvar(N,1);
lambda = sdpvar(N,1);
z = intvar(2,1);
x = sdpvar(2,1);
Inside = [z == V*lambda, lambda >=0, sum(lambda)==1];
Corner = [x == V*delta, sum(delta)==1];
solvesdp([Inside,Corner], (x-z)'*(x-z))
hold on
grid
plot(double(x(1)),double(x(2)),'ob')
plot(double(z(1)),double(z(2)),'pk')

But as I said, I am pretty sure this can be solved more efficiently directly somehow

• Thanks! I forgot to specify that I'm working in C++ but I've edited the original post. I agree that it seems like using optimization techniques is breaking out heavy machinery, but for the time being I'd be happy just being able to compute correct answers. – geomnerd Mar 26 '13 at 7:49
• It seems a shame to have to throw a full-fledged MIQP solver at this, I agree; but I don't know of another way around it. – Michael Grant Mar 26 '13 at 13:15

Both linear and integer programming can be solved in fast polynomial time if the dimension is fixed. Here is an example reference for the integer case:

Eisenbrand, "Fast integer programming in fixed dimension", 2003.