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. 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. 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.