Plane constraints in R3

Question

I have multiple plane constraints in $\mathbb{R}^3$ of the form:

$$n_i \cdot x \ge \delta_i$$

Where $n_i$ is the $i$th plane normal (in form (x, y, z)), $x$ is a point in space, and $\delta_i$ is the plane constant (distance of plane from origin).

I want to find $x$ such that the above is satisfied and the distance from $x$ to some other point $p$ is minimized. I believe this is a convex quadratic programming problem (the union of half spaces is always convex, IIRC). Normally I'd try to find an off-the-shelf library that could solve generalized quadratic programming problems, but I need a solution I can hand over to a client, and most of the solvers I know are either license restrictive (GPL or proprietary), gigantic libraries or both (CGAL).

I just need to solve this very specific problem form, and I'm willing to be relatively obtuse about how I do it. Does anyone know of any resources or way to approach this? I'm willing to take an off-the-shelf solver if it's approximately a single file and MIT-license equivalent, and I'm willing to write my own based on an algorithm description.

Answer 1

It's a quadratic program with linear inequality constraints, which are efficiently solved by the active set method. I pose implementing this method as a homework in my optimization courses and it doesn't take my students more than at most a couple of hours, assuming you have a solver for the linear systems you need to solve in each iteration.

I would suggest to just follow the description in Nocedal-Wright (Nocedal and Wright: "Numerical Optimization"). The algorithm is sufficiently simple that you will probably not spend much more time to implement it than to get used to the interface of some other library. Even better, if one of the other answers and comments turns up a free QP solver, then wrapping the active set method around it is pretty straight-forward.

Answer 2

If I'm reading your description correctly, you want to solve:

\begin{align} & \min_{x} \|x - p\|^{2}_{2} \\ \textrm{s.t.} & n_{i} \cdot x \geq \delta_{i}, i = 1,\ldots, N, \end{align}

assuming you have $N$ hyperplanes. As you point out, this program is a convex quadratic program (QP).

If you know Python, you could try using the BSD-licensed constrained QP (or nonlinear programming) solvers in SciPy and OpenOpt.

Since this problem is a convex QP, it can also be reformulated as a semidefinite program (SDP). DSDP has a BSD-like license, so you could try using that solver as well (but read the license to make sure; IANAL).

Regrettably, optimization software lends itself to making money, so any good software is either commercial, or a "free" license comes with many restrictions (GPL, academic use, non-commercial use, personal use, etc.)

I would not be surprised if there were a fast algorithm that exploited your particular problem structure.

Answer 3

Not a great solution, but what I'm doing right now:

Since the problem is in $R^3$, there can be at most 3 active plane constraints (inequalities turned equalities) at the correct solution. Any more, and at least one of them would be redundant. So I iterate through all possible combinations of 3, 2, 1, and 0 planes from among the $n$ inequalities ($\tbinom{n}{3}$ + $\tbinom{n}{2}$ + $\tbinom{n}{1}$ + 1) and find the one that produces a point that satisfies all other inequality constraints and is closest to the target point (or actually, I transform everything so the target point is the origin, and find the candidate point with the smallest norm that satisfies all constraints).

Given up to three plane constraints (the number of plane constraints being $d$), I stack them row by row to produce a $d$x$3$ matrix $A$ and similarly to produce a $d$x$1$ vector $\delta$. I construct the system $A x = \delta$, and solve it in a least squares way using normal equations. eg: $x = (A^TA)^+A^T\delta$. I find the pseudoinverse of $A^TA$ using a spectral decomposition for 3x3 symmetric matrices based on quaternions (see clicky). That is, $(A^TA) = Q D Q^T$, for an orthogonal matrix $Q$ and a diagonal matrix $D$, so $(A^TA)^+ = Q D^+ Q^T$. $D^+_{i,i} = 1/D_{i,i}$ if $D_{i,i} > 0$, and $0$ otherwise.

To check if a point satisfies all the plane inequalities I just check it against all planes one at a time.

...

This produces an $O(n^4)$ algorithm: $O(n^3)$ combinations of 3 planes producing $O(n^3)$ points that need to be checked. Each check needs to be against all $O(n)$ points.

It's not particularly clever, and for large $n$ it's going to be too slow, but it was relatively simple to code and test and was able to leverage existing quaternion and matrix3 code I had.

Is there anything more clever I could be doing that doesn't require more sophisticated mathematical tools? I don't have a general $n$x$n$ matrix library to leverage, just quaternions, 3x3 matrices, 4x4 matrices, and 3d/4d vectors.