Skip to content

Instantly share code, notes, and snippets.

@silgon

silgon/minimize_circle_ellipse_size_points.md

Last active June 8, 2022 16:34
Show Gist options
  • Save silgon/3800dcf531d4558c56fef80088b33526 to your computer and use it in GitHub Desktop.
Save silgon/3800dcf531d4558c56fef80088b33526 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
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(A@points.T + 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment