Minimum encompassing circle and ellipse
1. Minimum encompassing circle
This looks like a very easy problem to formulate and solve, not worth spending much time on. I think there are some interesting angles, you may not have thought about before.
Given \(n\) points \({\color{darkblue}p}_i\) (with their \(x\) and \(y\) coordinates), find the smallest circle identified by it center coordinates \({\color{darkred}c}\) and its radius \({\color{darkred}r}\), by minimizing \({\color{darkred}r}\), under the constraint that all points \({\color{darkblue}p}_i\) are inside the circle.
| Smallest circle model |
|---|
| \[\begin{align}\min_{{\color{darkred}r},{\color{darkred}c}}\> & {\color{darkred}r} \\ & {\bf{dist}}({\color{darkblue}p}_i,{\color{darkred}c}) \le {\color{darkred}r} && \forall i\\ & {\color{darkred}r}\ge 0 \\ & i \in \{1,\dots,n\} \end{align}\] |
| Practical NLP Model |
|---|
| \[\begin{align}\min\> & {\color{darkred}R} \\ & \sum_{k}({\color{darkblue}p}_{i,k}-{\color{darkred}c}_{k})^2 \le {\color{darkred}R} && \forall i \\ & i \in \{1,\dots,n\} \\ & k \in \{x,y\} \end{align}\] |
Here \({\color{darkred}R} = {\color{darkred}r}^2\). Using \({\color{darkred}R}\) saves us a square root operation (which is a bit dangerous close to zero: its derivative is not defined there). A good starting point can be found as follows: \[\begin{align}&{\color{darkred}c}_k := \frac{1}{n}\sum_i {\color{darkblue}p}_{i,k} \\ &{\color{darkred}R} := \max_i \sum_k ({\color{darkblue}p}_{i,k}-{\color{darkred}c}_k)^2 \end{align}\]
---- 87 PARAMETER results initial nlp c .x 41.710 50.425 c .y 43.748 43.448 r . 69.349 62.342 time . 0.015 iterations. 2.000
 |
| Initial and final, optimal circle with \(n=100\) data points |
Data points
---- 20 PARAMETER p data points x y point1 17.175 84.327 point2 55.038 30.114 point3 29.221 22.405 point4 34.983 85.627 point5 6.711 50.021 point6 99.812 57.873 point7 99.113 76.225 point8 13.069 63.972 point9 15.952 25.008 point10 66.893 43.536 point11 35.970 35.144 point12 13.149 15.010 point13 58.911 83.089 point14 23.082 66.573 point15 77.586 30.366 point16 11.049 50.238 point17 16.017 87.246 point18 26.511 28.581 point19 59.396 72.272 point20 62.825 46.380 point21 41.331 11.770 point22 31.421 4.655 point23 33.855 18.210 point24 64.573 56.075 point25 76.996 29.781 point26 66.111 75.582 point27 62.745 28.386 point28 8.642 10.251 point29 64.125 54.531 point30 3.152 79.236 point31 7.277 17.566 point32 52.563 75.021 point33 17.812 3.414 point34 58.513 62.123 point35 38.936 35.871 point36 24.303 24.642 point37 13.050 93.345 point38 37.994 78.340 point39 30.003 12.548 point40 74.887 6.923 point41 20.202 0.507 point42 26.961 49.985 point43 15.129 17.417 point44 33.064 31.691 point45 32.209 96.398 point46 99.360 36.990 point47 37.289 77.198 point48 39.668 91.310 point49 11.958 73.548 point50 5.542 57.630 point51 5.141 0.601 point52 40.123 51.988 point53 62.888 22.575 point54 39.612 27.601 point55 15.237 93.632 point56 42.266 13.466 point57 38.606 37.463 point58 26.848 94.837 point59 18.894 29.751 point60 7.455 40.135 point61 10.169 38.389 point62 32.409 19.213 point63 11.237 59.656 point64 51.145 4.507 point65 78.310 94.575 point66 59.646 60.734 point67 36.251 59.407 point68 67.985 50.659 point69 15.925 65.689 point70 52.388 12.440 point71 98.672 22.812 point72 67.565 77.678 point73 93.245 20.124 point74 29.714 19.723 point75 24.635 64.648 point76 73.497 8.544 point77 15.035 43.419 point78 18.694 69.269 point79 76.297 15.481 point80 38.938 69.543 point81 84.581 61.272 point82 97.597 2.689 point83 18.745 8.712 point84 54.040 12.686 point85 73.400 11.323 point86 48.835 79.560 point87 49.205 53.356 point88 1.062 54.387 point89 45.113 97.533 point90 18.385 16.353 point91 2.463 17.782 point92 6.132 1.664 point93 83.565 60.166 point94 2.702 19.609 point95 95.071 33.554 point96 59.426 25.919 point97 64.063 15.525 point98 46.002 39.334 point99 80.546 54.099 point100 39.072 55.782
Conopt Log
CONOPT 4 53.1.0 6a03d3b9 Feb 16, 2026 WEI x86 64bit/MS Window C O N O P T version 4.39.0 Copyright (C) GAMS Software GmbH GAMS Development Corporation Will use up to 4 threads. The user model has 100 constraints and 3 variables with 300 Jacobian elements, 200 of which are nonlinear. The Hessian of the Lagrangian has 2 elements on the diagonal, 0 elements below the diagonal, and 2 nonlinear variables. Iter Phase Ninf Infeasibility RGmax NSB Step InItr MX OK 0 0 0.0000000000E+00 (Input point) There are 100 penalty and minimax constraints with 1 variables. Preprocessed model has 100 constraints and 3 variables with 300 Jacobian elements, 200 of which are nonlinear. Iter Phase Ninf Infeasibility RGmax NSB Step InItr MX OK 0.0000000000E+00 (Full preprocessed model) 0.0000000000E+00 (After scaling) 0.0000000000E+00 (After adjusting penalty variables) ** Feasible solution. Value of objective = 4809.27124686 Iter Phase Ninf Objective RGmax NSB Step InItr MX OK 1 3 3.8865468687E+03 1.9E+03 3 1.0E+00 2 T T 2 4 3.8865468687E+03 0.0E+00 0 ** Optimal solution. There are no superbasic variables.
As this is a convex quadratic problem, we can formulate a conic model. This would allow us to solve the model with a different set of solvers. This exercise is always a bit of a chore, as the format of a SOCP-constraint (SOCP=Second-Order Cone Program) is very strict. Solvers expect an SOCP constraint of the form: \[\sum_i {\color{darkred}x}_i^2 \le {\color{darkred}y}^2\] with \({\color{darkred}x}_i\) free variables and \({\color{darkred}y}\ge 0\). So, following this to the letter, such a conic model can look like:
| Conic model |
|---|
| \[\begin{align}\min\> & {\color{darkred}r} \\ & {\color{darkred}\Delta}_{i,k} = {\color{darkblue}p}_{i,k}-{\color{darkred}c}_k && \forall i \\ & \sum_k{\color{darkred}\Delta}_{i,k}^2 \le {\color{darkred}r}^2 && \forall i\\ & {\color{darkred}r}\ge 0\\ &{\color{darkred}\Delta}_{i,k} \>{\bf{free}}\\& i \in \{1,\dots,n\} \\ & k\in \{x,y\} \end{align}\] |
| Conic model v2 |
|---|
| \[\begin{align}\min\> & {\color{darkred}r} \\ & {\color{darkred}\Delta}_{i,k} = {\color{darkblue}p}_{i,k}-{\color{darkred}c}_k && \forall i \\ & {\color{darkred}s}_i = {\color{darkred}r} &&\forall i \\ & \sum_k{\color{darkred}\Delta}_{i,k}^2 \le {\color{darkred}s}_i^2 && \forall i\\ & {\color{darkred}r}\ge 0 \\ & {\color{darkred}s}_i \ge 0\\ &{\color{darkred}\Delta}_{i,k} \>{\bf{free}}\\& i \in \{1,\dots,n\} \\ & k\in \{x,y\} \end{align}\] |
This model gives the same results as the NLP model. For conic models we don't really need an initial point. But, the modeling is a bit less natural, as we need to shoehorn the model into the exact format expected by the solver. This requires extra variables and constraints.
Obviously, we can reduce the number of points by only considering the convex hull. For the above data set with \(n=100\) points, the convex hull only has 14 points.
 |
| Initial and optimal circle using \(n=14\) convex hull points |
Complete GAMS Model
$ontext
Given n data points, find the smallest circle that contains all the points.
$offtext
option nlp=conopt, qcp=cplex;
* 0: use all data points * 1: use only convex hull $set convexhull 0
* 0: no plot * 1: produce html plot $set htmlplot 0
*----------------------------------------------------- * data *-----------------------------------------------------
set i 'data points' /point1*point100/ k 'coordinates' /x,y/ ;
parameter p(i,k) 'data points'; p(i,k) = uniform(0,100); display p;
*----------------------------------------------------- * compute the convex hull *-----------------------------------------------------
set hull(i) 'convex hull';
$if %convexhull%==0 $goto skiphull
embeddedCode Python: import scipy as sp import numpy as np import gams.transfer as gt
print("Computing convex hull...")
i = list(gams.get("i")) p = gt.Container(gams.db)["p"].toDense() hull = sp.spatial.ConvexHull(p) h = [i[pt] for pt in hull.vertices] gams.set("hull",h) endEmbeddedCode hull
parameter numpoints(*) 'number of points'; numpoints('data points') = card(i); numpoints('convex hull') = card(hull); display numpoints,hull;
$label skiphull
set j(i) 'subset j of points i actually used in model'; $if %convexhull%==0 j(i) = yes; $if %convexhull%==1 j(hull) = yes;
*----------------------------------------------------- * model 1: nlp model *-----------------------------------------------------
variables c(k) 'center' r2 'squared radius' ;
equation inside(i) 'point i is inside circle';
inside(j).. sum(k, sqr(p(j,k)-c(k))) =l= r2;
* (very good) initial point c.l(k) = sum(j,p(j,k))/card(i); r2.l = smax(j,sum(k, sqr(p(j,k)-c.l(k))));
model m1 /inside/;
*----------------------------------------------------- * solve and reporting *-----------------------------------------------------
parameter results(*,*,*); results('c',k,'initial') = c.l(k); results('r','','initial') = sqrt(r2.l);
solve m1 minimizing r2 using nlp;
results('c',k,'nlp') = c.l(k); results('r','','nlp') = sqrt(r2.l); results('time','','nlp') = m1.resusd; results('iterations','','nlp') = m1.iterusd; display results;
*----------------------------------------------------- * model 2: conic model *-----------------------------------------------------
variable r 'radius' diff(i,k) 'p(i,k)-c(k)' ; r.lo = 0;
positive variable s(i) 'needed for mosek' ;
equations ediff(i,k) 'auxiliary constraint, needed for socp' eaux(i) 'needed for mosek' inside2(i) 'point i is inside circle' ;
ediff(j,k).. diff(j,k) =e= p(j,k)-c(k); eaux(j).. r =e= s(j); inside2(j).. sum(k, sqr(diff(j,k))) =l= sqr(s(j));
model m2 /ediff,eaux,inside2/;
*----------------------------------------------------- * solve and reporting *-----------------------------------------------------
solve m2 minimizing r using qcp;
results('c',k,'socp') = c.l(k); results('r','','socp') = r.l; results('time','','socp') = m2.resusd; results('iterations','','socp') = m2.iterusd; display results;
*--------------------------------------------------------------------- * visualization *---------------------------------------------------------------------
$set html plot.html $set data data.js
$if %htmlplot%==0 $goto skipplot
parameter init(k), initr; init(k) = sum(j,p(j,k))/card(j); initr = sqrt(smax(j,sum(k, sqr(p(j,k)-init(k)))));
file fdata /%data%/; put fdata;
* points put "px=["; loop(j,put p(j,'x'):0:4,","); put "];"/; put "py=["; loop(j,put p(j,'y'):0:4,","); put "];"/; put "xa0=",(init('x')-initr):0:4,";"/; put "ya0=",(init('y')-initr):0:4,";"/; put "xa1=",(init('x')+initr):0:4,";"/; put "ya1=",(init('y')+initr):0:4,";"/; put "xb0=",(c.l('x')-r.l):0:4,";"/; put "yb0=",(c.l('y')-r.l):0:4,";"/; put "xb1=",(c.l('x')+r.l):0:4,";"/; put "yb1=",(c.l('y')+r.l):0:4,";"/;
$onecho > %html% <html> <script src="https://cdn.plot.ly/plotly-3.4.0.min.js" charset="utf-8"></script> <script src="%data%" charset="utf-8"></script> <h1>Smallest Encompassing Circle</h1> <div id="plotDiv"></div> <script> var data = { x: px, y: py, mode: 'markers', type: 'scatter', name: 'data points' };
var layout = { autosize: false, width: 650, height: 600, showlegend: true, shapes: [{type:'circle', xref:'x', yref:'y', x0:xa0, y0:ya0, x1:xa1, y1:ya1, line:{color:'orange',width:2}, fillcolor:'rgba(0,0,0,0)' }, {type:'circle', xref:'x', yref:'y', x0:xb0, y0:yb0, x1:xb1, y1:yb1, line:{color:'purple',width:2}, fillcolor:'rgba(0,0,0,0)' } ]
} var trc = [data]; var options = {staticPlot: true, displayModeBar: false}; Plotly.newPlot('plotDiv', trc, layout, options); </script> </html> $offecho
executetool 'win32.ShellExecute "%html%"';
$label skipplot |
2. Minimum enclosing ellipse
An ellipse represented by its center \(c\) and a positive definite matrix \(M\), can be defined by the set \[\{x|(x-c)^T M (x-c)\le 1\}\] Some special cases: when \(M=I\) we have a circle with radius \(r=1\). Using \(M=rI\), we described a circle with radius \(r\).
The area of this ellipse is \[\mathit{Area}=\frac{\pi}{\sqrt{{{\bf{det}}}M}}\] where \({\bf{det}}M\) is the determinant of \(M\). There is no easy, closed formula for the circumference (it is a complex elliptic integral). There are some good approximations however.
A more convenient representation for an ellipse is[1]: \[\{x|\> ||A x+b|| \le 1\}\] again with \(A\) being positive definite. The area is \[\mathit{Area}=\frac{\pi}{{\bf{det}}A}\] Using this, the minimum enclosing ellipse problem can be stated as the following convex semi-definite programming problem:
| SDP Model |
|---|
| \[\begin{align}\max\> & \log {\bf{det}} {\color{darkred}A} \\ & || {\color{darkred}A} {\color{darkblue}p}_i + {\color{darkred}b} || \le 1 && \forall i\\ & {\color{darkred}A}\succeq 0 \end{align}\] |
The wavy great-than-or-equal symbol \(\succeq\) indicates that \({\color{darkred}A}\) must be positive semidefinite. A slight relaxation from our earlier "positive definite." I am not worrying about that. The variable \({\color{darkred}b}\) is not the center. However we can recover the center by \[c:=-{\color{darkred}A}^{-1} {\color{darkred}b}\] This must be done in post-processing. It is not so easy to incorporate this into the model as a constraint, as this would introduce a non-convexity.
Derivation. We can go from the first to the second representation by taking \(M = A^TA\) and forming \((x-c)^T A^TA (x-c) = \left(A(x-c)\right)^T\left(A(x-c)\right) = ||Ax-Ac||^2\). We easily recognize \(-Ac=b\).
The \(\log {\bf{det}} {\color{darkred}A}\) objective is guaranteed to be concave. This way the model is proven to be convex and we can simply use CVXPY to solve this. CVXPY will only accept problems it can prove to be convex. Sometimes it will refuse models that are convex, but is unable to verify convexity.
Example
First we generate \(n=50\) data points, drawn from a bivariate normal distribution. That should give us a data set that make sense to describe by an ellipse. Our data looks like:
The second step is optional: identify the convex hull. This will reduce the number of constraints in the optimization model.
The number of points we need to consider is now reduced from \(n=50\) to \(m=12\). We feed this into our optimization model, expressed in CVXPY:
d = 2 # dimension of the space A = cp.Variable((d, d), PSD=True) b = cp.Variable(d) constraints = [cp.norm(A @ points2[i] + b) <= 1 for i in range(m)] objective = cp.Maximize(cp.log_det(A)) prob = cp.Problem(objective, constraints) prob.solve(verbose=True) print(f"Optimal objective value: {prob.value}") print(f"A:\n{A.value}") print(f"b:\n{b.value}") center = -np.linalg.inv(A.value) @ b.value print(f"center:\n{center}")
As \({\color{darkred}A}\) is symmetric, we only have \(3+2=5\) variables in this CVXPY model. However, we can see from the log that the actual size of the model that is fed to solver, is much larger:
This behavior is quite common in when dealing with convex problems. Often auxiliary variables and constraints are needed.
After solving, we have:
The solution is:
Optimal objective value: -2.8896363124545594 A: [[ 0.2584063 -0.08959489] [-0.08959489 0.24621565]] b: [ 0.02996678 -0.1443177 ]
CVXPY notebook
DIY SOCP Formulation
Before we used CVXPY to do the heavy lifting and transform this into a SDP problem. Here, I'll try to develop a model we can solve with a solver like cplex: no SDP solver, but we can do a rotated second order cone.
For our 2d case, \({\color{darkred}A}\) is \(2\times 2\), so we can write \[{\bf{det}}{\color{darkred}A} = {\color{darkred}a}_{1,1} {\color{darkred}a}_{2,2} - {\color{darkred}a}_{1,2}^2\] The trick is to form: \[\max {\color{darkred}t},\> {\color{darkred}t}^2 + {\color{darkred}a}_{1,2}^2 \le {\color{darkred}a}_{1,1} {\color{darkred}a}_{2,2}\] This is rotated cone. We need to add the bounds \({\color{darkred}a}_{1,1}, {\color{darkred}a}_{2,2} \ge 0\). This is related to requiring that \({\color{darkred}A}\) is positive definite.
The constraint \(||{\color{darkred}A} {\color{darkblue}p}_i +{\color{darkred}b}|| \le 1\) is ok. We just need to follow the accepted input format and write something like: \({\color{darkred}\alpha}_x^2 +{\color{darkred}\alpha}_y^2 \le 1\) using auxiliary variables \({\color{darkred}\alpha}\).
When we feed this into GAMS/Cplex, we see:
---- 118 VARIABLE a11.L = 0.258 VARIABLE a12.L = -0.090 VARIABLE a22.L = 0.246 ---- 118 VARIABLE b.L x 0.030, y -0.144 ---- 118 VARIABLE t.L = 0.236 objective PARAMETER detA = 0.056 determinant of A PARAMETER area = 56.507
GAMS Model
$onText
min area = π/det(A) ||Ap + b || <= 1 A pos def
1. A is symmetric, so we just use a11,a12,a22
2. Objective replace obj by max sqrt(det(A)) then we can form max t, t^2 + a12^2 <= a11*a22, a11,a22>=0 this is a rotated second order cone
3. ||Ap + b|| <= 1 is ok: it is convex. Use intermediate variables to convey this to the solver. 4. A is pos.def. Constraints: a11,a22 >= 0 a11*a22 - a12^2 >= 0 This is already taken care of by 2.
$offText
option qcp = cplex; *---------------------------------------------------------------------- * data *----------------------------------------------------------------------
set i 'points' /point1*point50/ c 'coordinates' /x,y/ h(i) 'subset: convex hull' ;
parameter p(i,c) 'points';
$onEmbeddedCode Python:
# returns # points p(i,c) drawn from bivariate normal distribution # subset h(i) points that are part of the convex hull
import numpy as np import scipy as sp
seti = list(gams.get("i")) n = len(seti) setc = list(gams.get("c"))
seed = 1234 print(f"{n=} {seed=}")
μ = np.array([0, 0]) Σ = np.array([[5, 4.5], [4.5, 6]])
rng = np.random.default_rng(seed) points = rng.multivariate_normal(μ, Σ, size=n)
gams.set("p",[(seti[i],setc[j],points[i,j]) for i in range(n) for j in range(2)])
ipoints2 = sp.spatial.ConvexHull(points).vertices gams.set("h",[(seti[i]) for i in ipoints2])
$offEmbeddedCode p h
display p,h;
set j(i) 'used points'; *j(i) = yes; j(h) = yes;
*---------------------------------------------------------------------- * socp model *----------------------------------------------------------------------
variable t 'objective' * matrix A is 2x2 and symmetric a11,a12,a22 ap1(i),ap2(i) b(c) ; * diagonal elements of A should be >= 0 positive variable a11,a22;
equations rcone 'rotated cone for objective' q(i) '||Ap+b||<=1' defap1(i) 'first element of Ap+b' defap2(i) 'second element of Ap+b' ;
rcone.. sqr(t) + sqr(a12) =l= a11*a22;
defap1(i).. ap1(i) =e= a11*p(i,'x')+a12*p(i,'y') + b('x'); defap2(i).. ap2(i) =e= a12*p(i,'x')+a22*p(i,'y') + b('y'); q(i).. sqr(ap1(i)) + sqr(ap2(i)) =l= 1;
model m /all/;
solve m maximizing t using qcp;
*---------------------------------------------------------------------- * reporting *----------------------------------------------------------------------
scalar detA 'determinant of A' area ;
detA = a11.l*a22.l-sqr(a12.l); area = pi / detA;
display a11.l,a12.l,a22.l,b.l,t.l,detA,area; |
References
- Stephen Boyd, Lieven Vandenberghe, Convex Optimization, section 8.4.
No comments:
Post a Comment