I have a set of points $P$ and want to find the ellipse with the smallest area that covers at least a fraction $f$ of these points. How can I do this?
These questions ask the same thing, but folks have only provided approximate answers.
This question asks for the smallest ellipse that bounds all the points ($f=1$).
I don't want an approximate answer. I really do want the smallest ellipse.
You asked for the smallest ellipse. An ellipse so small that its smallestness needs to be italicized. Others have provided answers that identify smallish ellipses, but, as Miracle Max says, "Mostly small is slightly large." Let's find the actual smallest ellipse.
To do so, we'll use mathematical programming. In "standard" programming we write code that progressively builds the solutions we're looking for. In mathematical programming, we specify a mathematical description of what we want and let powerful tools use that description to find a suitable solution.
Though this isn't exhaustive, you can think of mathematical programming as a set of nested classes of problem types that trade-off between speed, problem size, and generality. Each class of problems contains the easier, more specific classes as a special case. Linear Programming (LP) is very easy. We know everything about it and can solve problems of millions of variables in short order. Unfortunately, as the name suggests, linear programming is limited to linear equalities and inequalities. Quadratic programming (QP) is harder and the size of the problems we can reasonably solve is smaller. Second-order cone programming (SOCP) is harder still. The most general class of programming for which we have reliable methods is semi-definite programming (SDP). SDPs can be very challenging to solve.
LP QP SOCP SDP EXP POW MIP CBC X X GLPK X GLPK_MI X X OSQP X X CPLEX X X X X NAG X X X ECOS X X X X GUROBI X X X X MOSEK X X X X X X X* CVXOPT X X X X SCS X X X X X X SCIP X X X X XPRESS X X X X SCIPY X (*) Except mixed-integer SDP.
(From CVXPY documentation)
As we use harder and more general techniques, the number of software packages available to solve our problems decreases and their cost increases.
Unfortunately, at first glance, your problem is an SDP.
Note that an ellipsoid can be defined by a positive semi-definite matrix $Q$ and an offset $\vec{b}$ via the equation $$ \lVert Q\vec{x}+\vec{b}\rVert_2\le1 $$
The volume of the ellipsoid defined by the matrix $Q$ above is given by $\det(Q^{-1})$. Therefore, we would like to build an optimization problem with the objective $$ \min \det(Q⁻¹) $$ Taking inverses in optimization problems is bad, but if we apply a logarithm we can make the inverse go away giving us the new optimization problem: $$ \min -\log \det(Q) $$ The value we find for $-\log \det(Q)$ will not be the same value we would get for $\det(Q^{-1})$, but if we minimize $-\log \det (Q)$, then we will have also minimized $\det(Q^{-1})$.
The matrix determinant above suggests that we're trying to find a PSD matrix $Q$, which is worrisome because this would make our problem SDP---one of the hardest class of problems to solve. Fortunately, through a few tricks we can convert this problem into a second-order cone program (SOCP).
To do so, we leverage the fact that we're working in 2D where a PSD matrix has the form
$$
Q=\begin{bmatrix}
x & z \\
z & y
\end{bmatrix}
$$
The determinant of this matrix is
$$
det(Q) = xy-z^2
$$
Therefore, minimizing $-\sqrt{\det(Q)}$ is equivalent to minimizing a variable $-u$ subject to the constraint
$$
xy - z^2 >= u^2
$$
If we manipulate this a bit we get
$$
\begin{aligned}
xy - z^2 &>= u^2 \\
xy &>= u^2 + z^2 \\
x &>= (u^2 + z^2)/y
\end{aligned}
$$
We do this manipulation because, in order to program the optimization problem we need to put it into a form understood by the domain-specific language we'll use to describe the problem. In this case, we'll use cvxpy (because I don't know how to do this in R) which has a function called quad_over_lin which matches the form of the left-hand side ($(u^2 + z^2)/y$) of the above equation.
The optimization problem we've developed so far is:
minimize -u
such that
x >= quad_over_lin((u, z), y)
We need to incorporate our data points into this.
To do so, we return to our definition of an ellipsoid $$ \lVert Q\vec{x}+\vec{b}\rVert_2\le1 $$ If we introduce this constraint for each point $\vec{p_i}$, but use the same $Q$ and $\vec{b}$ for each, then we'll be requiring that the same ellipsoid contain all the points and optimizing for the ellipsoid which does this with minimal volume.
This gives us the problem
minimize -u
such that
x ≥ quad_over_lin((u, z), y)
||Q pᵢ +b ||₂ ≤ 1 for all points pᵢ
But you didn't want all the points to be contained by the ellipse, only some of them.
So let's introduce a set of $N$ binary variables $o$ (for off). If $o_i=0$ then we'll require that the point $\vec{p_i}$ be included in the ellipse defined by $(Q,\vec{b})$; if $o_i=1$, then the point need not be included in the ellipse.
To represent this, we return to our ellipse constraints: $$ \lVert Q\vec{p_i}+\vec{b}\rVert_2\le1 $$ and modify them like so: $$ \lVert Q\vec{p_i}+\vec{b}\rVert_2\le1+1000o_i $$ Now, if $o_i=1$, the constraint is easy to satisfy regardless of the values of $Q$ and $\vec{b}$; however, if $o_i=0$, then the constraint matches the original requirement that $p_i$ be in the ellipse.
The value 1,000 chosen above is arbitrary and represents a "big M constraint". Choosing this value wisely can impact the speed of the solution, as explained here and here.
Now if we have $N$ points and we want some fraction $f$ of them to be included in the ellipse, this is the same as saying that we want: $$ \sum_i o_i\le \lceil(1-f)N\rceil $$
The SOCP class of problems is efficient to solve because the problems are convex, which means that any local minimum (or maximum) is also the global minimum or maximum.
Introducing the binary variable above destroyed the convexity of this problem and brought us into the world of mixed integer programming where there are no algorithms that can solve all classes of problems efficiently. Therefore, our expectation is that for large problem sizes it will be challenging to find exact answers, even if we can find suitable answers quickly.
In actuality, the situation is even more challenging: we're working with a mixed-integer second-order cone program.
Consulting the chart of solvers above indicates a limited set of solvers that can handle such problems: CPLEX, Gurobi, MOSEK, SCIP, and XPRESS. By process of elimination, we need to use XPRESS.
pip3 install xpress. In my testing, XPRESS also solved the problem much faster than our alternative, SCIP.A working implementation is below using XPRESS as a solver and cvxpy as a domain-specific language.
Exploring the entire space of binary variables would take $O(2^N)$ time, so we know that modest increases in $N$ can have, in the worst case, an exponentially detrimental effect on wall-time. The reason we're able to explore spaces of $N=100$ or so is that our MISOCP solver (XPRESS) is able to cleverly prove that certain combinations of points will never be optimal and, therefore, doesn't have to explore them.
The number of points included might have a U-shaped effect on wall-time. If we call the number of points we want included in the ellipse $k$, then there are $$ {n \choose x} = \frac{n!}{r!(n-r)!} $$ unique combinations to check. This equation has a parabolic shape
implying that the most challenging problems will be for intermediate numbers of points include in the ellipse.
Intuitively, this makes sense. If we have 100 points and want 99 included in the ellipse, then the problem is easily solved by looking at $O(N)$ solutions. If we have 100 points and want only 1 included, then the problem's solved in $O(1)$ time by simply choosing any point.
As a quick sanity check we generate a dispersed cloud of data with a large number of points near the middle. Any ellipse excluding a small enough number of points should include the dense region. This is implemented in the sanity_check() function below.
We can then run for various problem sizes:
The wall-times for these results running on a 12-core Intel(R) Xeon(R) E-2176M CPU @ 2.70GHz were
N Include Wall-time
100 80 32s
100 20 5325s
100 10 909s
50 20 62s
this roughly matches the intuition we developed with the theory, though it's surprising that (100,80) takes much less time than (100,20) despite the number of combinations to be considered being nominally equal. This indicates that if we want to understand the problem better, we'll need to think harder about the theory.
If we look at the diagnostic output of the solver, we see:
Node BestSoln BestBound Sols Active Depth Gap GInf Time
Will try to keep branch and bound tree memory usage below 23.0GB
E 1822 -3.986475 -1462.273140 2 1122 22 99.73% 0 2
E 1908 -4.012773 -1462.273140 3 1240 29 99.73% 0 2
E 2317 -4.022067 -1462.273140 4 1372 43 99.72% 0 2
E 3170 -4.025301 -1462.273140 5 1840 57 99.72% 0 3
E 3491 -4.119326 -1462.273140 6 2022 76 99.72% 0 3
* 4210 -4.303862 -1462.273140 7 2247 84 99.71% 0 3
E 4226 -4.486650 -1462.273140 8 2247 80 99.69% 0 3
E 344248 -4.502819 -134.526887 9 175810 65 96.65% 0 181
* 440097 -4.531829 -110.826023 10 210721 86 95.91% 0 212
* 516576 -4.707305 -101.282451 11 237464 69 95.35% 0 239
E 561831 -4.833001 -97.186140 12 246522 81 95.03% 0 255
E 640486 -4.886923 -91.351188 13 269620 64 94.65% 0 287
E 817133 -4.957405 -80.818377 14 331871 68 93.87% 0 361
* 968180 -5.315599 -74.068716 15 381383 61 92.82% 0 424
E2708500 -5.427174 -12.636859 16 716078 88 57.05% 0 1236
*5250416 -5.819478 -8.857125 17 1049483 81 34.30% 0 2116
*8767053 -5.880111 -7.351761 18 1001197 98 20.02% 0 3241
*10587923 -5.920791 -6.935110 19 931248 81 14.63% 0 3850
*11736596 -5.974790 -6.726221 20 841014 83 11.17% 0 4270
The gap indicates the difference between the best solution (BestSoln) found so far and the bound we have on the best solution (BestBound). If the gap goes to 0% then we've proved that the BestSoln we've found is the best solution there is. The rightmost column indicates running time. After additional computation stages, the total wall-time was 5325s. The above indicates that if we accept a solution that is at most 20% greater than optimal, then we could have saved 20 minutes.
This suggests a general strategy to approach computationally intensive problems like this: determine what level of optimality guarantee you need and request your solver stop when it finds that.
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
def find_minimal_area_ellipse(pts: np.ndarray, num_to_include: int) -> None:
"""Find the minimal area ellipse containing at least `num_to_include` points
Args:
pts (np.ndarray): Points array, must be of size (N,2)
num_to_include (int): Minimal number of points to include in the ellipse
"""
assert pts.shape[1] == 2
# Big number used to deactivate constraints
off_offset = 1000
# Number of points
N = pts.shape[0]
# off[i]=1 indicates a point doesn't need to be in the ellipse
# off[i]=0 indicates a point must be in the ellipse
off = cp.Variable(N, boolean=True)
# Surrogate for the area of the ellipse
u = cp.Variable()
# Parameters describing the ellipse. We don't use a
# matrix to prevent cvxpy from indavertently
# generating an SDP problem.
x = cp.Variable()
z = cp.Variable()
y = cp.Variable()
# Offset of the ellipse
b = cp.Variable(2)
# Objective value
obj = cp.Minimize(-u)
# Constrains u to be proportional to the volume of the ellipse
constraints = [x >= cp.quad_over_lin(cp.hstack((z,u)), y)]
# Matrix describing the ellipsoid
emat = cp.bmat([[x,z],[z,y]])
# Transpose matrix for ease of use
pts = pts.T
# Ensure each point is bounded by the ellipse (subject to `off`)
for i in range(N):
constraints.append(cp.norm(emat@pts[:,i] + b) <= 1 + off_offset*off[i])
# Ensure that enough points are included in the ellipse
constraints += [ cp.sum(off) <= N - num_to_include ]
# Set up the problem
prob = cp.Problem(obj, constraints)
# Solve the problem
optval = prob.solve(solver='XPRESS', verbose=True)
print(f"Optimum value: {optval}")
# plot the ellipse and data
angles = np.linspace(0, 2*np.pi, 200)
rhs = np.row_stack((np.cos(angles) - b.value[0], np.sin(angles) - b.value[1]))
ellipse = np.linalg.solve([[x.value, z.value], [z.value, y.value]], rhs)
off_values = off.value.astype(bool)
pts_excluded = pts[:, off_values]
pts_included = pts[:, np.logical_not(off_values)]
plt.scatter(pts_excluded[0,:], pts_excluded[1,:])
plt.scatter(pts_included[0,:], pts_included[1,:])
plt.plot(ellipse[0,:].T, ellipse[1,:].T, c="red")
plt.xlabel('Dimension 1'); plt.ylabel('Dimension 2')
plt.title('Minimum Volume Ellipsoid')
plt.show()
def sanity_check() -> None:
"""
Generates a dataset in which the correct ellipse
should contain extra points
"""
N = 100
Npack = 90
# Generate a big cloud of data
pts = np.random.normal(loc=0, scale=1, size=(N,2))
# Make most of the data cluster in a small area
pts[0:Npack,:] = np.random.normal(loc=0, scale=0.01, size=(Npack,2))
find_minimal_area_ellipse(pts, 95)
def display_capability(N, num_include: int) -> None:
"""
Solves a problem defined by a number of points `N`
and the number of points `num_include` that should
be inside the ellipse
"""
pts = np.random.random((N,2))
find_minimal_area_ellipse(pts, num_include)
# Run some problems
np.random.seed(27182)
sanity_check()
display_capability(100, 80)
np.random.seed(27182)
display_capability(100, 20)
np.random.seed(27182)
display_capability(50, 20)
xpress==9.4.0 (works fine for ==9.3.1). The reason seems to be a bug in the deprecation of the old addVariable method.
$\endgroup$
Commented
Apr 25 at 16:37