polynomial_optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{\\Re}{\\mathbb{R}} \\newcommand{\\EE}{\\mathbb{E}}$\n",
    "# Polynomial Optimization and the Generalized Moment Problem\n",
    "\n",
    "## Polynomial Optimization\n",
    "\n",
    "The polynomial optimization problem is to\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\text{minimize} \\quad &f(x)\\\\\n",
    "    \\text{subject to} \\quad &g_i(x) \\geq 0,\\quad i=1,\\ldots,N\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "where $x \\in \\Re^d$ and $f(x), g_i(x) \\in \\Re[x]$ are polynomials. We implement the techniques developed by Parrilo et al. and Lasserre et al. where this problem is relaxed to a semidefinite program. \n",
    "\n",
    "In the first part of this worksheet, we show how to solve polynomial optimization problems using mompy where several example problems in the [GloptiPoly paper](http://homepages.laas.fr/henrion/papers/gloptipoly3.pdf) and the [solution extraction paper](http://homepages.laas.fr/henrion/papers/extract.pdf) are reproduced.\n",
    "\n",
    "In the second part, we give the simplest example of estimating a mixture model by solving the Generalized Moment Problem (GMP). Some extensions can be found in [the extra examples worksheet](extra_examples.ipynb). This software was used for our paper [Estimating mixture models via mixtures of polynomials](http://papers.nips.cc/paper/5702-estimating-mixture-models-via-mixtures-of-polynomials.pdf) where we look at some more elaborate settings for the mixture model problem.\n",
    "\n",
    "Other implementation for solving the polynomial optimization problem / and Generalized Moment Problem are [GloptiPoly](http://homepages.laas.fr/henrion/software/gloptipoly/) and [SOSTOOLS](http://www.cds.caltech.edu/sostools/) in MATLAB, and [ncpol2sdpa](https://pypi.python.org/pypi/ncpol2sdpa/) in python described [here](http://arxiv.org/abs/1308.6029).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# we are dependent on numpy, sympy and cvxopt.\n",
    "import numpy as np\n",
    "import cvxopt\n",
    "import mompy as mp\n",
    "\n",
    "# just some basic settings and setup\n",
    "mp.cvxsolvers.options['show_progress'] = False\n",
    "from IPython.display import display, Markdown, Math, display_markdown\n",
    "sp.init_printing()\n",
    "\n",
    "def print_problem(obj, constraints = None, moment_constraints = None):\n",
    "    display_markdown(mp.problem_to_str(obj,constraints,moment_constraints, False), raw=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1) First example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "Minimizing $\\mathcal{L}(x^{2} + y^{2})$\n",
       "\n",
       "subject to \t $x + y \\geq 4$, \t $x + y \\leq 4$, \t\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the maximum degree appearing in the problem is 2\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=1,\t rank=1,\t size=3,\t obj=8.000\n",
      "lost 0.0000001\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{x: array([ 1.99999971]), y: array([ 1.99999971])}"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x,y = sp.symbols('x,y')\n",
    "f =  x**2 + y**2\n",
    "gs = [x+y>=4, x+y<=4]\n",
    "print_problem(f, gs)\n",
    "sol = mp.solvers.solve_GMP(f, gs)\n",
    "mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2) Unconstrained optimization: the six hump camel back function\n",
    "\n",
    "A plot of this function can be found at [library of simutations](http://www.sfu.ca/~ssurjano/camel6.html). MATLAB got the solution $x^* = [0.0898   -0.7127]$ and corresponding optimal values of $f(x^*) = -1.0316$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "Minimizing $\\mathcal{L}(\\frac{x_{1}^{6}}{3} - 2.1 x_{1}^{4} + 4 x_{1}^{2} + x_{1} x_{2} + 4 x_{2}^{4} - 4 x_{2}^{2})$\t subject to no constraints"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the maximum degree appearing in the problem is 6\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=1,\t rank=3,\t size=10,\t obj=-1.032\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{x1: array([ 0.08984255, -0.08984255]), x2: array([-0.71265247,  0.71265247])}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1,x2 = sp.symbols('x1:3')\n",
    "f = 4*x1**2+x1*x2-4*x2**2-2.1*x1**4+4*x2**4+x1**6/3\n",
    "print_problem(f)\n",
    "sol = mp.solvers.solve_GMP(f, rounds=1)\n",
    "mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3) Multiple rounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Generally more rounds are needed to get the correct solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "Minimizing $\\mathcal{L}(- \\left(x_{1} - 1\\right)^{2} - \\left(x_{1} - x_{2}\\right)^{2} - \\left(x_{2} - 3\\right)^{2})$\n",
       "\n",
       "subject to \t $- \\left(x_{1} - 1\\right)^{2} + 1 \\geq 0$, \t $- \\left(x_{1} - x_{2}\\right)^{2} + 1 \\geq 0$, \t $- \\left(x_{2} - 3\\right)^{2} + 1 \\geq 0$, \t\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the maximum degree appearing in the problem is 2\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=1,\t rank=3,\t size=3,\t obj=-3.000\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=2,\t rank=6,\t size=6,\t obj=-3.000\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=3,\t rank=7,\t size=10,\t obj=-2.000\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=4,\t rank=8,\t size=15,\t obj=-2.000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{x1: array([ 2.00000811,  2.00000032,  0.99999978]),\n",
       " x2: array([ 3.00000234,  2.00000006,  2.00000007])}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1,x2,x3 = sp.symbols('x1:4')\n",
    "f = -(x1 - 1)**2 - (x1 - x2)**2 - (x2 - 3)**2\n",
    "gs = [1 - (x1 - 1)**2 >= 0, 1 - (x1 - x2)**2 >= 0, 1 - (x2 - 3)**2 >= 0]\n",
    "print_problem(f, gs)\n",
    "sol = mp.solvers.solve_GMP(f, gs, rounds=4)\n",
    "mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yet another example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "Minimizing $\\mathcal{L}(- 2 x_{1} + x_{2} - x_{3})$\n",
       "\n",
       "subject to \t $x_{1} \\geq 0$, \t $x_{1} \\leq 2$, \t $x_{2} \\geq 0$, \t $x_{3} \\geq 0$, \t $x_{3} \\leq 3$, \t $x_{1} + x_{2} + x_{3} \\leq 4$, \t $3 x_{2} + x_{3} \\leq 6$, \t $4 x_{1}^{2} - 4 x_{1} x_{2} + 4 x_{1} x_{3} - 20 x_{1} + 2 x_{2}^{2} - 2 x_{2} x_{3} + 9 x_{2} + 2 x_{3}^{2} - 13 x_{3} + 24 \\geq 0$, \t\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the maximum degree appearing in the problem is 2\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=1,\t rank=4,\t size=4,\t obj=-6.000\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=2,\t rank=8,\t size=10,\t obj=-5.692\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=3,\t rank=15,\t size=20,\t obj=-4.068\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=4,\t rank=17,\t size=35,\t obj=-4.000\n",
      "{x3: array([  3.48007643e-07,   2.99999775e+00]), x1: array([ 1.99999986,  0.50000288]), x2: array([ -1.10400714e-08,   1.46357718e-05])}\n"
     ]
    }
   ],
   "source": [
    "x1,x2,x3 = sp.symbols('x1:4')\n",
    "f = -2*x1 + x2 - x3\n",
    "gs = [0<=x1, x1<=2, x2>=0, x3>=0, x3<=3, x1+x2+x3<=4, 3*x2+x3<=6,\\\n",
    "      24-20*x1+9*x2-13*x3+4*x1**2 - 4*x1*x2+4*x1*x3+2*x2**2-2*x2*x3+2*x3**2>=0]; hs = [];\n",
    "print_problem(f, gs)\n",
    "sol = mp.solvers.solve_GMP(f, gs, hs, rounds=4)\n",
    "print mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4) Motzkin polynomial\n",
    "The Motzkin polynomial is non-negative, but cannot be expressed in sum of squares. It attains global minimum of 0 at $|x_1| = |x_2| = \\sqrt{3}/3$ (4 points). The first few relaxations are unbounded and might take a while for cvxopt to realize this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "Minimizing $\\mathcal{L}(x_{1}^{2} x_{2}^{2} \\left(x_{1}^{2} + x_{2}^{2} - 1\\right) + 0.037037037037037)$\t subject to no constraints"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the maximum degree appearing in the problem is 6\n",
      "slacks\n",
      "1e-06\n",
      "status: unknown\n",
      "round=1,\t rank=10,\t size=10,\t obj=-1119264077.823\n",
      "slacks\n",
      "1e-06\n",
      "status: unknown\n",
      "round=2,\t rank=15,\t size=15,\t obj=-408995.200\n",
      "slacks\n",
      "1e-06\n",
      "status: unknown\n",
      "round=3,\t rank=21,\t size=21,\t obj=-5690.864\n",
      "slacks\n",
      "1e-06\n",
      "status: unknown\n",
      "round=4,\t rank=28,\t size=28,\t obj=-918.009\n",
      "slacks\n",
      "1e-06\n",
      "status: unknown\n",
      "round=5,\t rank=36,\t size=36,\t obj=-67.576\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=6,\t rank=13,\t size=45,\t obj=0.000\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=7,\t rank=16,\t size=55,\t obj=0.000\n"
     ]
    }
   ],
   "source": [
    "x1,x2 = sp.symbols('x1:3')\n",
    "f = x1**2 * x2**2 * (x1**2 + x2**2 - 1) + 1./27\n",
    "print_problem(f)\n",
    "sol = mp.solvers.solve_GMP(f, rounds=7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{x1: array([-0.57735018, -0.5773504 ,  0.57735014,  0.57735044]), x2: array([ 0.57735087, -0.57735087, -0.57735112,  0.57735112])}\n"
     ]
    }
   ],
   "source": [
    "print mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 4, maxdeg=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generalized Moment Problem (GMP)\n",
    "The GMP is to\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\text{minimize} \\quad &f(x)\\\\\n",
    "    \\text{subject to} \\quad &g_i(x) \\geq 0,\\quad i=1,\\ldots,N\\\\\n",
    "    \\quad &\\mathcal{L}(h_j(x)) \\geq 0,\\quad j=1,\\ldots,M\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "where $x \\in \\Re^d$ and $f(x), g_i(x), h_j(x) \\in \\Re[x]$ are polynomials, and $\\mathcal{L}(\\cdot): \\Re[x] \\mapsto \\Re$ is a linear functional. Loosely speaking, we would like a measure $\\mu$ to exist so that $\\mathcal{L}(h) = \\int h\\ d\\mu$.\n",
    "\n",
    "\n",
    "For details see:\n",
    "\n",
    "\n",
    "Jean B. Lasserre, 2011, Moments, Positive Polynomials and Their Applications\n",
    "\n",
    "Jean B. Lasserre, 2008, A semidefinite programming approach to the generalized problem of moments\n",
    " \n",
    "\n",
    "### Estimating mixtures of exponential distributions\n",
    "We are interested in estimating the parameters of a mixture model from some observed moments, possibly with constraints on the parameters. Suppose that we have a mixture of 2 exponential distributions with density function $p(x; \\beta) \\sim \\exp(-x/\\beta)/\\beta$ in for $x \\in [0,\\infty)$. Suppose that $\\pi_1 = \\pi_2 = 0.5$ and $\\beta_1 = 1, \\beta_2 = 2$. We are allowed to observe moments of the data $\\EE[X^n] = {n!}(\\pi_1 \\beta_1^n + \\pi_2 \\beta_2^n)$ for $n=1,\\ldots,6$. For more details on estimating mixture model using this method, see our paper [Estimating mixture models via mixtures of polynomials](http://papers.nips.cc/paper/5702-estimating-mixture-models-via-mixtures-of-polynomials.pdf) and see [the extra examples worksheet](extra_examples.ipynb) for another easy example on mixture of Gaussians."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "Minimizing $\\mathcal{L}(\\beta^{4} + \\beta^{2} + 1)$\n",
       "\n",
       "subject to \t$\\mathcal{L}(\\beta - 1.5) = 0$, \t$\\mathcal{L}(\\beta^{2} - 2.5) = 0$, \t$\\mathcal{L}(\\beta^{3} - 4.5) = 0$, \t$\\mathcal{L}(\\beta^{4} - 8.5) = 0$, \t"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the maximum degree appearing in the problem is 4\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=1,\t rank=2,\t size=3,\t obj=12.000\n",
      "lost 0.0000000\n",
      "{beta: array([ 1.,  2.])}\n"
     ]
    }
   ],
   "source": [
    "beta = sp.symbols('beta')\n",
    "beta0 = [1,2]; pi0 = [0.5,0.5]\n",
    "hs = [beta**m - (pi0[0]*beta0[0]**m+pi0[1]*beta0[1]**m) for m in range(1,5)]\n",
    "f = sum([beta**(2*i) for i in range(3)])\n",
    "# note that hs are the LHS of h==0, whereas gs are sympy inequalities\n",
    "print_problem(f, None, hs)\n",
    "\n",
    "sol = mp.solvers.solve_GMP(f, None, hs)\n",
    "print mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 2, tol = 1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we try to solve the problem with insufficient moment conditions, but extra constraints on the parameters themselves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "Minimizing $\\mathcal{L}(\\beta^{4} + \\beta^{2} + 1)$\n",
       "\n",
       "subject to \t $\\beta \\geq 1$, \t\n",
       "$\\mathcal{L}(\\beta - 1.5) = 0$, \t$\\mathcal{L}(\\beta^{2} - 2.5) = 0$, \t"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the maximum degree appearing in the problem is 4\n",
      "slacks\n",
      "1e-06\n",
      "status: optimal\n",
      "round=1,\t rank=2,\t size=3,\t obj=12.000\n",
      "lost 0.0000001\n",
      "{beta: array([ 0.99999876,  1.99999896])}\n"
     ]
    }
   ],
   "source": [
    "gs=[beta>=1]\n",
    "f = sum([beta**(2*i) for i in range(3)])\n",
    "hs_sub = hs[0:2]\n",
    "print_problem(f, gs, hs_sub)\n",
    "\n",
    "sol = mp.solvers.solve_GMP(f, gs, hs_sub)\n",
    "print mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 2, tol = 1e-3)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}