Smallest sphere containing points using Fusion

In the previous post we formulated the fairly intuitive problem

How to find the smallest sphere that encloses a set of $k$ points $p_i \in R^n$.

as a second order conic problem:

$$\begin{equation}   \begin{aligned} \min r_0&&&\\ s.t.&&&\\ & r_i = r_0,&\quad& i=1,\ldots,k\\ & w_i = p_i - p_0,&\quad& i=1,\ldots,k\\ & (r_i,w_i) \in Q^{(n+1)}, &\quad& i=1,\ldots,k. \end{aligned} \end{equation}$$


We will show in this post how the MOSEK Fusion API for conic optimization allows for a very compact prototyping. Let's start!

First, we assume that the points are represented as row vectors, i.e. $p_0=(p_0^1,\ldots,p_0^n)$. The given points are organized in a matrix $P$ by row, i.e. $P\in R^{k\times n}$

We reformulate the problem using a different, yet equivalent, formulation using the fact, see [1], that an affine transformation preserve the conic representability of a set, and thus we can write:

$$\begin{equation} r= r_0\mathbf{e}, \quad w = P - \mathbf{e}p_0. \end{equation}$$

being $\mathbf{e}\in R^n$ a vector with all components equal to one. The formulation we are looking for is

$$\begin{equation} \begin{aligned} \min r_0&\\  s.t.&\\  &\quad ( r_0\mathbf{e}, P-\mathbf{e}p_0^T) \in \Pi_{i=1}^k Q_k^{(n+1)}. \end{aligned} \end{equation}$$

Now let's input the model using the Python Fusion API, but for the other supported languages it would look like much the same. Assuming the given points are stored in a matrix p, the following few steps suffices (click here for the whole code):

1- create an optimization model

M = Model("minimal sphere enclosing a set of points")

view raw gistfile1.txt hosted with ❤ by GitHub

2- declare the variables (note $p_0$ is a row vector)

	r = M.variable("r", Domain.unbounded())
	p0 = M.variable("p_0", NDSet(1,n), Domain.unbounded())

view raw gistfile1.txt hosted with ❤ by GitHub

3- define the constraints in one shot (here np is shorthand for numpy)

	M.constraint('norm',\
	Expr.hstack([
	Expr.mul(np.ones(k),r),
	Expr.sub(Expr.mul(np.ones( (k,1) ), p0 ), DenseMatrix(p) )
	]),
	Domain.inQCone(k,n+1))

view raw gistfile1.txt hosted with ❤ by GitHub

 Note how we "stack" columns side by side obtaining a $k\times (n+1)$ matrix of expressions, which will be interpreted line by line to obtain the $k$ ice-cream cones.

4- define the objective function

M.objective('min_radius', ObjectiveSense.Minimize, r)

view raw gistfile1.txt hosted with ❤ by GitHub

5- solve !

M.solve()

view raw gistfile1.txt hosted with ❤ by GitHub

Does it work? Of course it does....here a solution for $n=1000000, k=2$.

Pretty neat right? Where is the magic and what are the pros and cons?

The magic is the ability of Fusion to convert the expressions in standard conic constraints, i.e. the first model we formulate: all the additional variables will be generated automatically by Fusion!

How does it look like the same code with the standard Pyhton API? Just take a look at the constraint declaration. First we need some offset in the variable set:

	numvar= 1 + n + k*n + k
	offset_r= 0
	offset_p= 1
	offset_w= 1+n
	offset_t= 1+n+n*k
	task.appendvars(numvar)
	task.putvarboundslice(0,numvar,[mosek.boundkey.fr for i in range(numvar)], numpy.zeros(numvar), numpy.zeros(numvar))

view raw gistfile1.txt hosted with ❤ by GitHub

Then a loop to create all the linear and conic constraints:

	numcon= 2*k + k
	task.appendcons(numcon)
	for i in range(k):#for each point
	task.putarow(i, [0,i+offset_t], [1., -1,]) #delta coord#radius aux
	task.putconbound(i, mosek.boundkey.fx, 0. ,0.)
	for j in range(n):
	task.putconbound(k+ i*n + j, mosek.boundkey.fx, p[i][j] , p[i][j])
	task.putarow(k+i*n + j,[offset_p + j, offset_w+i*n + j],[1.0,-1.0])
	task.appendcone(mosek.conetype.quad,0., np.hstack([ offset_t+i, [ offset_w+i*n + j for j in range(n)]] ))

view raw gistfile1.txt hosted with ❤ by GitHub

The code is still manageable one might say, but it is partially because Python allows for very compact notation. But the point is: could you easily deduce the model implemented by the code? And moreover, what if an offset goes wrong? To you the answers....

Fusion vs. Standard API

Pros:
- a cleaner and shorter code
- easier for the solver to analyze the structure of the model
- avoid the use of indexes and offset to input the problem data
- closer to a modeling language

Cons:
- require more modeling skills
- slower execution during the creation of the model, not the solver execution, at least for large models

We believe usually the pros exceed the cons, especially in the long run when users get used to the syntax and modeling philosophy.

Give Fusion a try!

[1] Ben-Tal, A., & Nemirovski, A. (2001). Lectures on modern convex optimization: analysis, algorithms, and engineering applications (Vol. 2). Siam.