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
2- declare the variables (note $p_0$ is a row vector)
3- define the constraints in one shot (here np is shorthand for numpy)
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
5- solve !
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:
Then a loop to create all the linear and conic constraints:
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.
