9
$\begingroup$

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.

$\endgroup$

1 Answer 1

13
$\begingroup$

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.

Mathematical Programming as a Technique

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.

Hierarchy of convex programming

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.

Ellipsoids are Positive Semi-definite Matrices

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

Ellipsoids have volumes

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})$.

Convert SDP to SOCP

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.

Point Constraints

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ᵢ

Including only some of the points

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

Mixed Integer Programming

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.

  • CPLEX and MOSEK are expensive and don't seem to have accessible academic/trial versions.
  • Gurobi is expensive. It does have academic and trial licenses, but getting these set up is annoying.
  • SCIP is free. However, it can be somewhat challenging to install and, in my tests, was much slower than XPRESS.
  • XPRESS is expensive, but it has a community license that allows you to experiment with it on smaller problem sizes. Most importantly, we can install XPRESS with that community license simply by typing pip3 install xpress. In my testing, XPRESS also solved the problem much faster than our alternative, SCIP.

Implementation

A working implementation is below using XPRESS as a solver and cvxpy as a domain-specific language.

Theory

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

n choose k

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.

Results

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.

Sanity check

We can then run for various problem sizes:

N=100, Include=80

N=100, Include=80

N=100, Include=20

N=100, Include=20

N=100, Include=10

N=100, Include=10

N=50, Include=20

N=50, Include=20

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.

Implementation (for real this time)

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)
$\endgroup$
2
  • 5
    $\begingroup$ Excellent write-up! Thanks for sharing this knowledge. $\endgroup$ Dec 4, 2021 at 21:12
  • 3
    $\begingroup$ There are not that many answers on this forum that go into such extreme detail. Fantastic job! $\endgroup$ Dec 6, 2021 at 4:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.