minimize_circle_ellipse_size_points.md

minimize_circle_ellipse_size_points.py

import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt

np.random.seed(0)
points = np.random.rand(8, 2)
r = cp.Variable()
constraints = [cp.norm(d) <= r for d in points]
problem = cp.Problem(cp.Minimize(r),constraints)
problem.solve()

circle1 = plt.Circle((0, 0), r.value, color='C1', fill=False, label="circle centered at zero")
plt.gca().add_patch(circle1)

r = cp.Variable()
x = cp.Variable(2)
constraints = [cp.norm(d-x) <= r for d in points]
problem = cp.Problem(cp.Minimize(r),constraints)
problem.solve()

circle2 = plt.Circle(x.value, r.value, color='C2', fill=False, label="circle")
plt.gca().add_patch(circle2)
plt.scatter(*points.T)

# ellipses inspired by:
# http://web.cvxr.com/cvx/examples/cvxbook/Ch08_geometric_probs/html/min_vol_elp_finite_set.html
# https://github.com/rmsandu/Ellipsoid-Fit/blob/main/max_inner_ellipsoid_v2.py
# https://notebook.community/stephenbeckr/convex-optimization-class/CVX_demo/cvxpy_intro

n, m = points.shape
A = cp.Variable((m,m), PSD=True)
cost = cp.Maximize(cp.log_det(A))
constraints = [ cp.norm(A@points[i,:]) <= 1 for i in range(n) ]
prob = cp.Problem(cost,constraints)
prob.solve()

theta = np.linspace(0, 2 * np.pi, 200)
sphere_pts = np.c_[np.cos(theta), np.sin(theta)]
ellipse = np.linalg.solve(A.value, sphere_pts.T)

plt.plot(ellipse[0, :], ellipse[1, :], c='C3', label="ellipse centered at zero")


n, m = points.shape
A = cp.Variable((m,m), PSD=True)
b = cp.Variable((m,1))
cost = cp.Maximize(cp.log_det(A))
# constraints = [ cp.norm(A@points[i,:] + b) <= 1 for i in range(n) ]
constraints = [ cp.norm([email protected] + b, 2, 0) <= 1 ]
prob = cp.Problem(cost,constraints)
prob.solve()

theta = np.linspace(0, 2 * np.pi, 200)
sphere_pts = np.c_[np.cos(theta)-b.value[0], np.sin(theta)-b.value[1]]
ellipse = np.linalg.solve(A.value, sphere_pts.T)

plt.plot(ellipse[0, :], ellipse[1, :], c='C4', label="ellipse")

plt.title("Minimal circles and ellipses")
plt.axis("equal")
plt.legend()
plt.show()