New Portfolio Optimization Cookbook chapters and notebooks

We are happy to share that our Portfolio Optimization Cookbook has been updated and new chapters has been added. We have also added new notebooks with code examples.
The additions cover among other things:

• CVaR
• EVaR
• Gaussian mixture return model
• Risk budgeting
• Robust optimization
• Multiperiod optimization

As usual the new concepts are implemented in Mosek Fusion for Python

Enjoy! 

SwissQuant Math Challenge on LASSO Regression: the MOSEK way

The Swiss company SwissQuant publishes monthly a math challenge in which they call for people to elaborate and discuss problems of different kind, somehow related to the company business. This month they focus on (quadratic) convex optimization, asking people to proposed specialized algorithm for the OLS-LASSO problem that can beat the established commercial code for quadratic programming as MOSEK. You can find details here Math Challenge March 2014.

Clearly we are glad they appreciate somehow our product, but also we feel challenged in tackling the problem itself. MOSEK is not specialized for OLS-LASSO, but still it can be very efficient we thought, especially in dual form. Moreover, we are always eager to find good test problems, for the core solver but also for our APIs.

The OLS-LASSO problem

So let's start from the formulation given in the challenge, where we are given $T$ samples in dimension $m$:

\begin{equation}\label{main}\tag{OLS-LASSO}
  \begin{aligned}
    \min_{u,b} & \frac{1}{2} || u ||^2 + \sum_{j=1}^{m} \lambda_j |b_j|\\
    &\\
    s.t.&\\
    & u = y - X b
  \end{aligned}
\end{equation}

where $b,\lambda\in \mathbb{R}^m, u,y\in \mathbb{R}^T$. We reformulate problem \ref{main} so to have a linear objective function:

\begin{equation*}
  \begin{aligned}
    \min_{u,b,z,w} &\quad z + \lambda^T w\\
    &\\
    s.t.&\\
    & u + X b = y\\
    & s=1\\
    & (s,z,u) \in Q_R^{T+2}\\
    & w_i \geq |b_i|& \forall i=1,\ldots,m.
  \end{aligned}
\end{equation*}

where $Q_R$ is the rotated cone. There are two way to to get rid of the absolute value:

(1) Linear constraints 

\[ w_i\geq |b_i| \Rightarrow \left\{ \begin{array}{l}  w \geq b \\ w \geq -b \\\end{array}\right.\]

(2) Second-order conic constraint

\[ w_i \geq |b_i| \Leftrightarrow w_i\geq \sqrt{b_i^2} \equiv (w_i,b_i)\in Q^2\]

We will see which primal form works the best. But it's the dual that we want to really see at work. Using the standard procedure outlined in our post about duality, we derive the dual of problem (PQ). It is easier than PL, believe us.

First let's fix a partition of the variable set as $[s,z,u,w,b]$, and clearly states the basic vector and
matrices:


(1) cost vector: $c^T= [0 \quad 1 \quad 0_T \quad \lambda^T \quad 0_m]$;
(2) right-hand side: $b=^T[1 \quad y^T]^T$;
(2) matrix coefficients: $A= \left[ \begin{array}{ccccc} 1 & 0&0_T&0_m&0_m\\ 0&0&I_T& 0_m&X\end{array}\right]$;


We then introduce as many dual variables as primal linear constraints:
$\mu_0$ for the first constraint, and $\mu\in \mathbb{R}^T$ for the
second one. We also include additional slack variables $\gamma$
partitioned in the same way as the primal variables, i.e. $\gamma^T=[\gamma_s,\gamma_z,\gamma_u,\gamma_w,\gamma_b]$.

The dual problem is

\begin{equation*}
  \begin{aligned}
    \max_{\mu_0, \mu, \gamma} & \quad \mu_0 + y^T \mu \\
    &\\
    s.t.&\\
    & c - A^T \left[ \begin{array}{c}\mu_0\\\mu\end{array} \right]= \gamma\\
    & (\gamma_s,\gamma_z,\gamma_u) \in Q_R^{T+2}\\
    & (\gamma_w,\gamma_b) \in \Pi_{i=1}^m Q^2.
  \end{aligned}
\end{equation*}

After some calculations, we obtain:

\begin{equation}\tag{D}
  \begin{aligned}
    \max_{\mu_0,\mu,\gamma_z} &  \quad \mu_0 + y^T \mu \\
    &\\
    s.t.&\\
    &  X^T \mu \leq  \lambda\\
    & \gamma_z=  1\\
    & \left(-\mu_0, \gamma_z, -\mu \right) \in Q_R^{2+T}.
  \end{aligned}
\end{equation}

Try to do the same with PL, you'll get the same model.

The results

We implement these three models using the MOSEK Fusion Python API and run a batch of tests for $m=[10,25,75,100]$ and $T$ from $10000$ to $100000$. To make it easy, $X,\lambda$ and $y$ are all uniformly generated in $[0,1]$. The code that interface with the model looks like


Here you can see the running time for problem PL in green, PQ in red and D in blue. In the first  plot we fix $m$ and let $T$ varies:

and then if we fix $T$ either to $10000$ or $100000$ and we get

Again, the dual is again a clear winner, as in a previous post! No way the primal can even get close to the dual.

But it is also interesting to see how  PQ is as good as PL. Red and green lines are almost indistinguishable. Only for larger $m$, i.e. when the number of cones increases, we can see a significant difference. Otherwise the difference is very very small.

How to write the dual problem in a nutshell

Duality in linear programming  is often mentioned and sometimes even used (see [2]). But in our experience, many practitioners unfortunately have a quite blurred idea  about the dual problem and the meaning of it. Things are even worse when considering conic programming, despite the fact that most of the theory remains unchanged. There are at least three good reasons to have a good understanding conic duality:

• The dual problem could be faster to solve (see our post), even if the complexity is the same.
• Sensitivity analysis uses dual information (see our post), just think about shadow prices.
• Infeasibility and unboundedness certificates are stated in terms of duality, using the Farkas' Lemma.

Duality theory for conic optimization is rich and well-established and deserves much more space than a blog post (see for instance [1]). As first step, we will show how to easily derive the dual formulation of a given conic optimization problem. For sake of clarity, we will not consider semidefinite optimization problems, even if the conic framework readily extends to them.

Let assume we are given an optimization problem, which we refer as the primal. If we are lucky, the primal can be reformulated as a conic optimization problem by means of some transformations (see [1],[3]). An hopefully small additional effort is usually needed to cast our conic problem in (primal) standard form. Let the primal variables be reordered and partitioned such that  $x=(x_1,x_2,\ldots,x_k)^T\in \mathbb{R}^n$ with $x_i\in \mathbb{R}^{n_i}$. Then the standard primal form is:

\begin{equation}
\begin{aligned}
  &\min \quad c^Tx &&\\
  & s.t. &&\\
  &&Ax - b = 0 &\\
  &&x_i \in \mathit{K_i^{n_i}} &\quad i=1,\ldots\,k
\end{aligned}
\end{equation}

Note how each block variable $x_i$ belongs to a cone $\mathit{K_i^{n_i}}$; moreover, we can  recover the standard LP formulation just setting $k=1$ and $n_1=n$. In short we may write $x\in \Pi_{i=1}^k K_i^{n_i}$. For our purpose, $\mathit{K_i^{n_i}}$ can be

- the positive orthant $\mathbb{R}^{n_i}$
- a Lorentz cone $L_i^{n_i}$ of dimension $n_i$
- a rotated Lorentz cone $R_i^{n_i}$ of dimension $n_i$

We introduce as many dual variables $y$ as cones in Problem (1) and obtain the dual problem as follows:

\begin{equation}
\begin{aligned}
  &\max \quad b^Ty \\
  & s.t. \\
  &&  A^T y -c \in\Pi_{i=1}^k \mathit{K}_i^*\\
\end{aligned}
\end{equation}

where $\mathit{K}_i^*$ is the dual cone of $\mathit{K}_i$. In practice, the cones we work with are all self-dual, i.e. they are all dual to themselves, and this means roughly speaking $\mathit{K}_i^*\equiv \mathit{K}_i$. Introducing slack variables Problem (2) reads:

\begin{equation}
\begin{aligned}
  &\max \quad b^Ty \\
  & s.t. \\
  && c- A^T y + s= 0\\
   &&s_i \in \mathit{K}^*_i  &\quad i=1,\ldots\,k
\end{aligned}
\end{equation}

which can be converted again in standard form by noticing that $\min x = - \max -x$ and transforming free variables in non negative ones (see for instance [2]).

Summarizing, we recommend to follow these steps:

1. if possible recast the given problem in conic form
2. transform the conic form in standard form
3. derive the dual formulation
4. simplify the dual if possible

In practice, step 1 requires the biggest effort, while the others follow quite directly. Indeed, the crucial point is to certified the given problem can be expressed in conic form, and we do this actually constructing such a representation. Once we get the conic formulation, most of the job is done. So, working in conic form has also the interesting side-effect of a much simpler access to dual formulation and information.

Check this post to see a practical example! Others will follow in the next days.

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

[2] Papadimitriou, C. H., & Steiglitz, K. (1998). Combinatorial optimization: algorithms and complexity. Courier Dover Publications.

[3] MOSEK modeling manual , 2013, Mosek ApS

Conic duality in practice: when dual beats the primal...

Despite its simplicity, the problem we have been discussing in the last two posts (see the first and second post) is a nice source of examples and insights about conic programming. The formulation we obtained 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^{(n+1)}.
\end{aligned}
\end{equation}

Recall that $p_0,p_i$ are row vectors and $e$ is a vector with all ones. We might wonder whether the dual problem of (1) is harder or easier to solve. Using Conic Programming duality this can be done pretty easily.

So let's use the duality theory of Conic Programming to derive the dual problem of (1).  First, to derive the dual we restate (1) in a slightly different way:

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

Following [1], this is our primal formulation and the construction of the dual problem should be quite straightforward! Consider the standard primal form (see pp. 69 in [1]):

\begin{equation}
\begin{aligned}
\min c^T x&\\
 s.t.&\\
 &\quad A_i x_i -b_i \in K_i & i=1,\ldots,k,
\end{aligned}
\end{equation}

where $K_i$'s are cones and $x_i=(r_o , p_0)^T$: problem (2) follows by setting $ c=( 1 , 0_{n})^T$ and

\[
A_i= I_{n+1},\quad b_i= \begin{pmatrix} 0 \\ p_i^T \end{pmatrix},\quad K_i=Q^{(n+1)} , \quad i=1,\ldots,k..
\]

To formulate the dual,  for each cone $i$ we introduce a vector $y_i \in R^{n+1}$ of dual variables such that $y_i$ must belong to the self-dual cone of $Q^{(n+1)}$, i.e. $Q^{(n+1)}$ itself. The dual problem takes the form

\begin{equation}
\begin{aligned}
\max & \sum_{i=1}^k b_i^T y_i\\
 s.t.&\\
 &\sum_{i=1}^k y_i = c\\
 &y_i\in Q^{(n+1)} & i=1,\ldots,k.
\end{aligned}
\end{equation}

which is equivalent to

\begin{equation}
\begin{aligned}
\max & \mathrm{Tr} \sum_{i=1}^k B Y^T \\
 s.t.&\\
 & e_k^T Y = c\\
 &Y \in \Pi_{i=1}^k Q^{(n+1)}
\end{aligned}
\end{equation}

where \[ Y=  \begin{pmatrix} y_1^T \\ \vdots \\ y_k^T \end{pmatrix}, B=  \begin{pmatrix} p_1^T \\ \vdots \\ p_k^T \end{pmatrix}, e_k= \begin{pmatrix} 1\\ \vdots \\ 1\end{pmatrix}\in R^k.\]

Let's implement the dual formulation (5) using the MOSEK Fusion API. The code looks like

The code is very simple, compact and self-commenting. Note how we also exploit the fact that the first column of $B$ is composed by zeros. What about performance? The dual beats the primal in every sense! Look at the plot below (tests have been done on my desktop machine Pentium(R) Dual-Core  CPU E5700  @ 3.00GHz under Linux) where we report the running time needed to build the optimization model and to solve it.

The dual is much faster both in terms of time spent in building the model, which is indeed much simpler, and also in terms of solver execution time! Both model building scale almost linearly, but also the solution time does. If the former is expected, the latter might come as a surprise! Note however, that model building time matters when tackling polynomially solvable problems....

Final take-away message:

- Solving the dual instead of the primal can be much more efficient!

- The conic formulation helps to derive the dual problem!

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

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

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.

A nice example of conic optimization modeling

Few days ago, my boss Erling Andersen suggested me to look at this problem to learn the MOSEK Fusion API:

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

The intuitive formulation is to minimize the radius r of the sphere centered in $p_0 \in R^n$, i.e.

\begin{align}
& \min r \\
& s.t. \\
& & || p_i - p_0||_2 \leq r, & i=1,\ldots,k
\end{align}

The inspiration came from a recent book¹ where the problem, probably taken from the GAMS repository², is presented as a difficult test for nonlinear local optimizers.

At the same time (what a coincidence!), O'Neil and Puget published interesting post about the Chebyshev centers of polygons...which turns out to be equivalent to what I was looking at! As Puget pointed out, the problem is known to be simple: Megiddo³ showed back in 1982 that it can be solved in linear time by a specialized algorithm.

Anyway, we will use it as an example of conic programming modeling example. First we introduce auxiliary variables such that  $w_i = p_i - p_0$, yielding

\begin{align}
& \min r \\
& s.t.\\
& & || w_i ||_2 \leq r, & i=1,\ldots,k\\
& & w_i = p_i - p_0, & i=1,\ldots,k
\end{align}

Constraints (6) define quadratic (Lorentz) cones $Q^{(n+1)}$, i.e.

\begin{equation} || w_i ||_2   \leq r    \Rightarrow (r, w_i) \in  Q^{(k+1)}.\end{equation}

For the second order cones to be disjointed, we introduce auxiliary variables $r_i, i=1,\ldots,n$ that must be equal to $r$. The final formulation of our problem is

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

which is a fairly simple conic quadratic optimization problem!

In the next post we will show how to easily implement this problem using the Mosek Fusion API.

¹ Neculai, Andrei, "Nonlinear Optimization Application Using the GAMS Technology", Springer 2013
² http://www.gams.com/testlib/libhtml/circlen.htm
³ Megiddo, Nimrod. "Linear-time algorithms for linear programming in R3 and related problems." Foundations of Computer Science, 1982. SFCS'08. 23rd Annual Symposium on. IEEE, 1982.

Regression techniques for Portfolio Optimisation using MOSEK

A new paper with the title: "Regression techniques for Portfolio Optimisation using MOSEK" has been published. The paper illustrates how to use MOSEK for real-world portfolio optimization problems.