notebooks/duality.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Duality properties and sums-of-squares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A simple problem with equality and inequality constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider a problem with both equality and inequality constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{array}{ll}\n",
    "\\text{minimize}   & x + z\\\\\n",
    "\\text{subject to} & x - z \\geq 0\\\\\n",
    "                  & x^2 + z^2 = 1.\n",
    "\\end{array}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Polyopt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x, z = variables([\"x\", \"z\"]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = x + z;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "g = x - z;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "h = x^2 + z^2 - 1;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "prob = momentprob(2, f, [g], [h]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Open file 'polyopt.task'\n",
      "Problem\n",
      "  Name                   :                 \n",
      "  Objective sense        : max             \n",
      "  Type                   : CONIC (conic optimization problem)\n",
      "  Constraints            : 15              \n",
      "  Cones                  : 0               \n",
      "  Scalar variables       : 7               \n",
      "  Matrix variables       : 2               \n",
      "  Integer variables      : 0               \n",
      "\n",
      "Optimizer started.\n",
      "Conic interior-point optimizer started.\n",
      "Presolve started.\n",
      "Linear dependency checker started.\n",
      "Linear dependency checker terminated.\n",
      "Eliminator started.\n",
      "Freed constraints in eliminator : 0\n",
      "Eliminator terminated.\n",
      "Eliminator - tries                  : 1                 time                   : 0.00            \n",
      "Lin. dep.  - tries                  : 1                 time                   : 0.00            \n",
      "Lin. dep.  - number                 : 0               \n",
      "Presolve terminated. Time: 0.00    \n",
      "Optimizer  - threads                : 4               \n",
      "Optimizer  - solved problem         : the primal      \n",
      "Optimizer  - Constraints            : 15\n",
      "Optimizer  - Cones                  : 1\n",
      "Optimizer  - Scalar variables       : 8                 conic                  : 8               \n",
      "Optimizer  - Semi-definite variables: 2                 scalarized             : 27              \n",
      "Factor     - setup time             : 0.00              dense det. time        : 0.00            \n",
      "Factor     - ML order time          : 0.00              GP order time          : 0.00            \n",
      "Factor     - nonzeros before factor : 120               after factor           : 120             \n",
      "Factor     - dense dim.             : 0                 flops                  : 3.36e+03        \n",
      "ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  \n",
      "0   2.0e+00  1.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  0.01  \n",
      "1   6.4e-01  3.2e-01  4.1e-01  2.86e-01   -4.653334109e-01  -6.237304963e-01  3.2e-01  0.01  \n",
      "2   1.4e-01  6.8e-02  3.2e-01  1.01e+00   -1.250925203e+00  -1.213892948e+00  6.8e-02  0.01  \n",
      "3   1.7e-02  8.3e-03  1.3e-01  1.13e+00   -1.389157294e+00  -1.384192283e+00  8.3e-03  0.01  \n",
      "4   1.6e-03  8.2e-04  3.8e-02  1.11e+00   -1.412333100e+00  -1.411934407e+00  8.2e-04  0.01  \n",
      "5   3.9e-04  2.0e-04  1.8e-02  1.05e+00   -1.413859444e+00  -1.413774675e+00  2.0e-04  0.01  \n",
      "6   1.1e-04  5.5e-05  8.8e-03  1.01e+00   -1.414055001e+00  -1.414038061e+00  5.5e-05  0.01  \n",
      "7   3.1e-05  1.5e-05  4.3e-03  1.01e+00   -1.414185489e+00  -1.414182344e+00  1.5e-05  0.01  \n",
      "8   7.8e-06  3.9e-06  2.1e-03  1.00e+00   -1.414201695e+00  -1.414201132e+00  3.9e-06  0.01  \n",
      "9   1.8e-06  9.2e-07  1.0e-03  1.00e+00   -1.414211646e+00  -1.414211554e+00  9.2e-07  0.01  \n",
      "10  4.4e-07  2.2e-07  4.8e-04  1.00e+00   -1.414212912e+00  -1.414212896e+00  2.2e-07  0.02  \n",
      "11  1.1e-07  5.5e-08  2.4e-04  1.00e+00   -1.414213443e+00  -1.414213440e+00  5.5e-08  0.02  \n",
      "12  2.9e-08  1.4e-08  1.2e-04  1.00e+00   -1.414213518e+00  -1.414213517e+00  1.4e-08  0.02  \n",
      "13  7.3e-09  3.9e-09  6.0e-05  9.99e-01   -1.414213555e+00  -1.414213555e+00  3.6e-09  0.02  \n",
      "Interior-point optimizer terminated. Time: 0.02. \n",
      "\n",
      "Optimizer terminated. Time: 0.03    \n",
      "\n",
      "Interior-point solution summary\n",
      "  Problem status  : PRIMAL_AND_DUAL_FEASIBLE\n",
      "  Solution status : OPTIMAL\n",
      "  Primal.  obj: -1.4142135548e+00   nrm: 1e+00    Viol.  con: 7e-09    var: 0e+00    barvar: 0e+00  \n",
      "  Dual.    obj: -1.4142135547e+00   nrm: 1e+00    Viol.  con: 0e+00    var: 2e-09    barvar: 4e-09  \n"
     ]
    }
   ],
   "source": [
    "X, Z, t, y, solsta = solve_mosek(prob);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We inspect the solution,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15×2 Array{Polyopt.Poly{Float64},2}:\n",
       " 1.0      1.0      \n",
       " z        -0.707154\n",
       " z^2      0.500067 \n",
       " z^3      -0.353624\n",
       " z^4      0.250067 \n",
       " x        -0.707059\n",
       " x*z      0.5      \n",
       " x*z^2    -0.353577\n",
       " x*z^3    0.250033 \n",
       " x^2      0.499933 \n",
       " x^2*z    -0.35353 \n",
       " x^2*z^2  0.25     \n",
       " x^3      -0.353482\n",
       " x^3*z    0.249967 \n",
       " x^4      0.249933 "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[prob.basis y]"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "and verify optimality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-7.409961533255682e-11"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t - Polyopt.evalpoly(f, y[[6,2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Polyopt.evalpoly(g, y[[6,2]]) >= 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "abs(Polyopt.evalpoly(h, y[[6,2]])) < 1e-5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Duality interpretation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The returned solution consists objective value of the relaxation,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.4142135547775996"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "of two semidefinite matrices $X_1$ and $X_2$, corresponding to $f(x,z)$ and $g(x,z)$,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Array{Float64,2},1}:\n",
       " [0.762453 0.549465 … -0.139263 -0.281168; 0.549465 0.949308 … -0.00651993 -0.080656; … ; -0.139263 -0.00651993 … 0.839117 -0.0583538; -0.281168 -0.080656 … -0.0583538 0.519432]\n",
       " [0.307644 0.312502 0.122507; 0.312502 0.474474 -0.0325033; 0.122507 -0.0325033 0.205813]                                                                                        "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a symmetric indefinite matrix $Z_1$ corresponding to $h(x,z$),"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{Array{Float64,2},1}:\n",
       " [-0.65176 -0.104357 -0.205958; -0.104357 -0.585633 0.0583538; -0.205958 0.0583538 -0.519432]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and $y$ from which we extracted the coefficients for the optimal monomial solution,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15-element Array{Float64,1}:\n",
       "  1.0     \n",
       " -0.707154\n",
       "  0.500067\n",
       " -0.353624\n",
       "  0.250067\n",
       " -0.707059\n",
       "  0.5     \n",
       " -0.353577\n",
       "  0.250033\n",
       "  0.499933\n",
       " -0.35353 \n",
       "  0.25    \n",
       " -0.353482\n",
       "  0.249967\n",
       "  0.249933"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dual interpretation expresses $f(x,z) - t = v_1^T X_1 v_1 + g(x,z) v_2^T X_2 v_2 + h(x,z) v_2^T Z_1 v_2$, in other words we express $f(x,z)-t$ as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$f(x,z) - t = s_0(x,z) + g(x,z)s_1(x,z) - h(x,z) w(x, z)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $s_0(x,z)$ and $s_1(x,z)$ are sums-of-squares, but $w(x,z)$ is not. Let us verify the expression,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Array{Polyopt.Poly{Int64},1}:\n",
       " 1  \n",
       " z  \n",
       " z^2\n",
       " x  \n",
       " x*z\n",
       " x^2"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v1 = monomials(2, [x,z])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Polyopt.Poly{Int64},1}:\n",
       " 1\n",
       " z\n",
       " x"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v2 = monomials(1, [x,z])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4142135547775996+z+x"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f-t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4142135583531192+0.9999999925133197*z+0.9999999999634588*x"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "truncate( dot(v1, X[1]*v1) + g*dot(v2, X[2]*v2) + h*dot(v2, Z[1]*v2), 1e-8 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}