# Find closest point, subject to linear inequality constraints

Given a point $p\in \mathcal{R}^2$, I want to compute the closest point $x \in \mathcal{R}^2$, subject to linear inequality constraints $Ax \leq b$. That is,

$$\begin{array}{ll} \text{minimize} & \|x-p\|_2\\ \text{subject to} & A x \leq b\end{array}$$

I believe that this can be done with an off-the-shelf quadratic programming solver, but I'm wondering if there is a more efficient algorithm specialized to two variables ($x \in \mathcal{R}^2$) and Euclidean distance ($\min \|x - p\|$).

I think that this question is a bit different than this one because I don't have a direct representation of the feasible region as a list of vertices.

• How many rows does $A$ have? Jul 2, 2016 at 11:51
• @Rodrigo $A$ can have an arbitrary number of rows. Jul 2, 2016 at 14:39
• @BmoreDaniel Have you heard of cvxgen.com? Solving quadratic programs really fast is what it does. Jul 3, 2016 at 0:57
• @RodrigodeAzevedo Thanks for the pointer. I did not know about cvxgen. I've requested a license, so I'll post here once I find out what code it produces for this problem. Jul 3, 2016 at 2:01
• There's nothing particularly unique about the algorithm in cvxgen. What makes it fast is that it unrolls the loops in the various linear algebra routines. Jul 4, 2016 at 3:57

This can be done in $O(N log(N))$ time, where $N$ is the number of rows in $A$.

The matrix inequality $Ax \leq b$ is equivalent to half-plane inequalities: $(a_i, x) \leq b_i, i = 1, \ldots, N$. We can obtain the vertices of polygon that is the intersection of these half-planes by the following procedure:

1. Transform lines $\{(a_i, x) = b_i\} \mapsto d_i$ using duality transformation: $\delta(\{a x + b y = 1\}) = (a, b)$.

2. Construct the convex hull of $d_i$ using one of known algorithms.

3. Map the vertices and sides of the obtained convex hull back to the primal space: $\delta((a, b)) = \{a x + b y = 1\}$, $\delta(\{a x + b y = 1\}) = (a, b)$ and obtain the polygon $C_1 C_2 \ldots C_N$ which is this the intersection of initial half-planes.

4. Now you know all the vertices of the polygon and can apply the above cited methods to your problem.

I'm not quite sure that this is the best solution, though.

If your task is repeatable, i.e. you have to repeat the problem for different $p$'s many times, you can construct the dynamic convex hull and find the distance in $O(log(N))$ for each query.

• I didn't know what the duality transformation was, so for anyone else in that position, I found the description within this article to be useful. Jul 8, 2016 at 20:53

In 2D the closest point problem reduces to two cases. Either it will be on a vertex or in one of the edges of the polyhedron described by $$Ax \leq b$$. And, of course, min distance will be $$0$$ if $$q$$ is inside this polyhedron. Converting to a vertex representation seems adequate in 2D, it will allow you to use the algorithms from computational geometry.

This problem can be set up as a minimization problem using the calculus of variations and then solved as a linear equation. Basically, you minimize the equation

$$J = \frac{1}{2}\sum_{i=1}^{N} (x_i - p_{i})\cdot(x_i - p_{i}) + \vec{\lambda}^T (A\vec{x} - \vec{b})$$

where $$\vec{x}$$ is the original vector of length $$N$$ and $$\vec{\lambda}$$ is a vector of unknown constants with a length equal to the number of equations in the set $$A\vec{x} = \vec{b}$$ (we will call this length $$M$$).

Once you have written this equation, take the partial derivatives of J with respect to the vectors $$\vec{x}$$ and $$\vec{\lambda}$$, set them equal to $$0$$

$$\begin{bmatrix} \frac{\partial J}{\partial \vec{x}} \\ \frac{\partial J}{\partial \vec{\lambda}} \end{bmatrix} = 0$$

and then solve for $$\vec{x}$$ and $$\vec{\lambda}$$. This should give you a linear equation of the form $$C\vec{y}=\vec{d}$$ where the values for $$C$$, $$\vec{y}$$, and $$\vec{d}$$ are respectively:

$$\begin{bmatrix} {I}_{NxN} & A^T \\ A & {0}_{MxM} \end{bmatrix} \begin{bmatrix} \vec{x} \\ \vec{\lambda} \end{bmatrix} = \begin{bmatrix} \vec{p} \\ \vec{b} \end{bmatrix}$$

$$N$$ = number of variables in $$\vec{x}$$
$$M$$ = number of constraint equations in $$A\vec{x}=\vec{b}$$

$$\vec{y} = C^{-1} \vec{d}$$ , and the first $$N$$ values of $$\vec{y}$$ will be your answer for $$x_i$$ (the remaining $$M$$ values will be the values of $$\vec{\lambda}$$, and are not actually necessary for your solution).

• Please tell me if there are any mistakes in my edit. By the way could you expand on "i dont usually answer these questions"? How come you do now? :) Mar 18, 2021 at 22:26
• I simply rarely see a question on these forums to which I have the answer and nobody else does. But I worked out this one as I was looking it up, and since nobody else had my solution listed... Thank you for helping me out with the edits; your work gave me the clues I needed to figure out how to write the rest of the equations properly (I don't usually use LaTeX either :-) Mar 19, 2021 at 17:40
• Yes! Latex is quite inutitive after like a half an hour! Mar 19, 2021 at 17:54