{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the symmetric polynomial matrix (taken from [Example 3.77, BPT13])\n",
    "$$P(x) =\n",
    "\\begin{bmatrix}\n",
    "    x^2 - 2x + 2 & x\\\\\n",
    "    x            & x^2\n",
    "\\end{bmatrix}.$$\n",
    "We could like to know whether $P(x)$ is positive semidefinite for all $x \\in \\mathbb{R}$.\n",
    "\n",
    "[BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.\n",
    "*Semidefinite Optimization and Convex Algebraic Geometry*.\n",
    "Society for Industrial and Applied Mathematics, **2012**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Polynomial{true,Int64},2}:\n",
       " x² - 2x + 2  x \n",
       " x            x²"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using DynamicPolynomials\n",
    "@polyvar x\n",
    "P = [x^2 - 2x + 2 x\n",
    "            x     x^2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A sufficient condition for a symmetric polynomial matrix $P(x)$ to be positive semidefinite for all $x$ is to the existence of a matrix $M(x)$ such that $P(x) = M^\\top(x) M(x)$. If such matrix $M$ exists, we say that the matrix is an \\emph{sos matrix} (see [Definition 3.76, BPT13]).\n",
    "While determining whether $P(x)$ is positive semidefinite for all $x$, is NP-hard (checking nonnegativity of a polynomial is reduced to this problem for $1 \\times 1$ matrices), checking whether $P(x)$ is an sos matrix is an sos program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using SumOfSquares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first need to pick an SDP solver, see [here](http://www.juliaopt.org/JuMP.jl/dev/installation/#Getting-Solvers-1) for a list of the available choices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using CSDP\n",
    "factory = with_optimizer(CSDP.Optimizer, printlevel=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using MosekTools\n",
    "factory = with_optimizer(Mosek.Optimizer, QUIET=true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OPTIMAL::TerminationStatusCode = 1"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = SOSModel(factory)\n",
    "mat_cref = @constraint(model, P in PSDCone())\n",
    "optimize!(model)\n",
    "termination_status(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While the reformulation of sos matrix to sos polynomial is rather simple, as explained in the \"Sum-of-Squares reformulation\" section below, there is a technical subtelty about the Newton polytope that if not handled correctly may result in an SDP of large size with bad numerical behavior. For this reason, it is recommended to model sos *matrix* constraints as such as will be shown in this notebook and not do the formulation manually unless there is a specific reason to do so.\n",
    "\n",
    "As we can verify as follows, only 4 monomials are used using the sos *matrix* constraint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element MonomialVector{true}:\n",
       " x##367\n",
       " x##369\n",
       " ##367 \n",
       " ##369 "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "certificate_monomials(mat_cref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sum-of-Squares reformulation\n",
    "\n",
    "One way to obtain the reduction to an sos program is to create an intermediate vector of variable $y$ and check whether the polynomial $p(x, y) = y^\\top P(x) y$ is sos.\n",
    "However, special care is required when approximating the Newton polytope of $p(x, y)$. Indeed, for instance if the entries of $P(x)$ are quadratic forms then the Newton polytope of $p(x, y)$ is the cartesian product between the Newton polytope of $y^\\top y$ and the Newton polytope of $x^\\top x$. In other words, $p(x, y)$ belongs to a family of quartic forms called biquadratic forms. This fact is important when generating the semidefinite program so that only bilinear monomials are used. So if the cheap outer approximation is used (instead of the exact polyhedral computation) for the newton polytope then it is important to use a multipartite computation approximation of the newton polytope. The multipartie exact approach may perform worse compared to the unipartite exact in certain cases though. Consider for instance the polynomial matrix $\\mathrm{Diag}(x_1^1, x_2^2)$ for which $p(x, y) = x_1^2y_1^2 + x_2^2y_2^2$. For this polynomial, only the monomials $x_1y_1$ and $x_2y_2$ are needed in the SDP reformulation while the multipartite approach, as it will compute the Newton polytope as a cartesian product, will not see the dependence between $x$ and $y$ in the presence of monomials and will also select the monomials $x_1y_2$ and $x_2y_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ x^{2}y_{1}^{2} + x^{2}y_{2}^{2} - 2xy_{1}^{2} + 2xy_{1}y_{2} + 2y_{1}^{2} $$"
      ],
      "text/plain": [
       "x²y₁² + x²y₂² - 2xy₁² + 2xy₁y₂ + 2y₁²"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@polyvar y[1:2]\n",
    "p = vec(y)' * P * vec(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$ (1)x^{2}y_{1}^{2} + (1)x^{2}y_{2}^{2} + (-2)xy_{1}^{2} + (2)xy_{1}y_{2} + (2)y_{1}^{2} \\text{ is SOS} $"
      ],
      "text/plain": [
       "(1)x²y₁² + (1)x²y₂² + (-2)xy₁² + (2)xy₁y₂ + (2)y₁² is SOS"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = SOSModel(factory)\n",
    "cref = @constraint(model, p in SOSCone())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OPTIMAL::TerminationStatusCode = 1"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optimize!(model)\n",
    "termination_status(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see below, 6 monomials (instead of the 4 monomials with the sos *matrix* constraint) have been used as the unipartite cheap outer approximation was used for the Newton Polytope."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element MonomialVector{true}:\n",
       " xy₁ \n",
       " xy₂ \n",
       " y₁y₂\n",
       " x   \n",
       " y₁  \n",
       " y₂  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "certificate_monomials(cref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.1.0",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}