$\newcommand{\Re}{\mathbb{R}} \newcommand{\EE}{\mathbb{E}}$
The polynomial optimization problem is to
$ \begin{align} \text{minimize} \quad &f(x)\\ \text{subject to} \quad &g_i(x) \geq 0,\quad i=1,\ldots,N \end{align} $
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.
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 and the solution extraction paper are reproduced.
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. This software was used for our paper Estimating mixture models via mixtures of polynomials where we look at some more elaborate settings for the mixture model problem.
Other implementation for solving the polynomial optimization problem / and Generalized Moment Problem are GloptiPoly and SOSTOOLS in MATLAB, and ncpol2sdpa in python described here.
In [1]:
# we are dependent on numpy, sympy and cvxopt.
import numpy as np
import cvxopt
import mompy as mp
# just some basic settings and setup
mp.cvxsolvers.options['show_progress'] = False
from IPython.display import display, Markdown, Math, display_markdown
sp.init_printing()
def print_problem(obj, constraints = None, moment_constraints = None):
display_markdown(mp.problem_to_str(obj,constraints,moment_constraints, False), raw=True)
In [2]:
x,y = sp.symbols('x,y')
f = x**2 + y**2
gs = [x+y>=4, x+y<=4]
print_problem(f, gs)
sol = mp.solvers.solve_GMP(f, gs)
mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 1)
Out[2]:
A plot of this function can be found at library of simutations. MATLAB got the solution $x^* = [0.0898 -0.7127]$ and corresponding optimal values of $f(x^*) = -1.0316$
In [3]:
x1,x2 = sp.symbols('x1:3')
f = 4*x1**2+x1*x2-4*x2**2-2.1*x1**4+4*x2**4+x1**6/3
print_problem(f)
sol = mp.solvers.solve_GMP(f, rounds=1)
mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 2)
Out[3]:
Generally more rounds are needed to get the correct solutions.
In [4]:
x1,x2,x3 = sp.symbols('x1:4')
f = -(x1 - 1)**2 - (x1 - x2)**2 - (x2 - 3)**2
gs = [1 - (x1 - 1)**2 >= 0, 1 - (x1 - x2)**2 >= 0, 1 - (x2 - 3)**2 >= 0]
print_problem(f, gs)
sol = mp.solvers.solve_GMP(f, gs, rounds=4)
mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 3)
Out[4]:
Yet another example
In [5]:
x1,x2,x3 = sp.symbols('x1:4')
f = -2*x1 + x2 - x3
gs = [0<=x1, x1<=2, x2>=0, x3>=0, x3<=3, x1+x2+x3<=4, 3*x2+x3<=6,\
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 = [];
print_problem(f, gs)
sol = mp.solvers.solve_GMP(f, gs, hs, rounds=4)
print mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 2)
In [6]:
x1,x2 = sp.symbols('x1:3')
f = x1**2 * x2**2 * (x1**2 + x2**2 - 1) + 1./27
print_problem(f)
sol = mp.solvers.solve_GMP(f, rounds=7)
In [7]:
print mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 4, maxdeg=3)
The GMP is to
$ \begin{align} \text{minimize} \quad &f(x)\\ \text{subject to} \quad &g_i(x) \geq 0,\quad i=1,\ldots,N\\ \quad &\mathcal{L}(h_j(x)) \geq 0,\quad j=1,\ldots,M \end{align} $
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$.
For details see:
Jean B. Lasserre, 2011, Moments, Positive Polynomials and Their Applications
Jean B. Lasserre, 2008, A semidefinite programming approach to the generalized problem of moments
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 and see the extra examples worksheet for another easy example on mixture of Gaussians.
In [8]:
beta = sp.symbols('beta')
beta0 = [1,2]; pi0 = [0.5,0.5]
hs = [beta**m - (pi0[0]*beta0[0]**m+pi0[1]*beta0[1]**m) for m in range(1,5)]
f = sum([beta**(2*i) for i in range(3)])
# note that hs are the LHS of h==0, whereas gs are sympy inequalities
print_problem(f, None, hs)
sol = mp.solvers.solve_GMP(f, None, hs)
print mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 2, tol = 1e-3)
Now we try to solve the problem with insufficient moment conditions, but extra constraints on the parameters themselves.
In [9]:
gs=[beta>=1]
f = sum([beta**(2*i) for i in range(3)])
hs_sub = hs[0:2]
print_problem(f, gs, hs_sub)
sol = mp.solvers.solve_GMP(f, gs, hs_sub)
print mp.extractors.extract_solutions_lasserre(sol['MM'], sol['x'], 2, tol = 1e-3)