{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Term sparsity"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Adapted from**: Example 3.5 of [WML20b]\n",
    "\n",
    "[WML20a] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre.\n",
    "*TSSOS: A Moment-SOS hierarchy that exploits term sparsity*.\n",
    "arXiv preprint arXiv:1912.08899 (2020).\n",
    "\n",
    "[WML20b] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre.\n",
    "*Chordal-TSSOS: a moment-SOS hierarchy that exploits term sparsity with chordal extension*.\n",
    "arXiv preprint arXiv:2003.03210 (2020)."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(DynamicPolynomials.PolyVar{true}[x₁, x₂, x₃],)"
     },
     "metadata": {},
     "execution_count": 1
    }
   ],
   "cell_type": "code",
   "source": [
    "using DynamicPolynomials\n",
    "@polyvar x[1:3]"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "We would like to find the minimum value of the polynomial"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "2x₁²x₂² + 142x₂²x₃² - 2x₁²x₂ + 18x₂²x₃ - 54x₂x₃² + x₁² - 2x₁x₂ + 3x₂² - 2x₂x₃ + 6x₃²",
      "text/latex": "$$ 2x_{1}^{2}x_{2}^{2} + 142x_{2}^{2}x_{3}^{2} - 2x_{1}^{2}x_{2} + 18x_{2}^{2}x_{3} - 54x_{2}x_{3}^{2} + x_{1}^{2} - 2x_{1}x_{2} + 3x_{2}^{2} - 2x_{2}x_{3} + 6x_{3}^{2} $$"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "poly = x[1]^2 - 2x[1]*x[2] + 3x[2]^2 - 2x[1]^2*x[2] + 2x[1]^2*x[2]^2 - 2x[2]*x[3] + 6x[3]^2 + 18x[2]^2*x[3] - 54x[2]*x[3]^2 + 142x[2]^2*x[3]^2"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "The minimum value of the polynomial can be found to be zero."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter:  7 Ap: 9.00e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  3.1796787e-13 \n",
      "Success: SDP solved\n",
      "Primal objective value: 0.0000000e+00 \n",
      "Dual objective value: 0.0000000e+00 \n",
      "Relative primal infeasibility: 6.80e-09 \n",
      "Relative dual infeasibility: 2.83e-50 \n",
      "Real Relative Gap: 0.00e+00 \n",
      "XZ Relative Gap: 8.00e-50 \n",
      "DIMACS error measures: 1.02e-08 0.00e+00 2.83e-50 0.00e+00 0.00e+00 8.00e-50\n",
      "CSDP 6.2.0\n",
      "Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00 \n",
      "Iter:  1 Ap: 9.71e-01 Pobj: -2.4909930e+02 Ad: 8.28e-01 Dobj:  7.2430886e+02 \n",
      "Iter:  2 Ap: 9.99e-01 Pobj: -1.6816440e+02 Ad: 9.32e-01 Dobj:  1.6018610e+02 \n",
      "Iter:  3 Ap: 9.96e-01 Pobj: -6.7498290e+00 Ad: 9.24e-01 Dobj:  6.1399335e+01 \n",
      "Iter:  4 Ap: 1.00e+00 Pobj: -1.4722921e+00 Ad: 9.01e-01 Dobj:  8.9329470e+00 \n",
      "Iter:  5 Ap: 1.00e+00 Pobj: -5.0478524e-01 Ad: 7.90e-01 Dobj:  2.7911688e+00 \n",
      "Iter:  6 Ap: 1.00e+00 Pobj: -1.5409216e-01 Ad: 7.76e-01 Dobj:  8.7431136e-01 \n",
      "Iter:  7 Ap: 1.00e+00 Pobj: -4.8559418e-02 Ad: 8.12e-01 Dobj:  2.4102522e-01 \n",
      "Iter:  8 Ap: 1.00e+00 Pobj: -1.2896403e-02 Ad: 8.39e-01 Dobj:  6.6196298e-02 \n",
      "Iter:  9 Ap: 1.00e+00 Pobj: -1.7033049e-03 Ad: 1.00e+00 Dobj:  1.0847092e-02 \n",
      "Iter: 10 Ap: 1.00e+00 Pobj: -2.4388656e-04 Ad: 1.00e+00 Dobj:  7.6100068e-04 \n",
      "Iter: 11 Ap: 9.96e-01 Pobj: -7.3346824e-06 Ad: 1.00e+00 Dobj: -4.5795945e-05 \n",
      "Iter: 12 Ap: 9.97e-01 Pobj: -2.2313479e-07 Ad: 1.00e+00 Dobj: -7.5985091e-05 \n",
      "Iter: 13 Ap: 1.00e+00 Pobj: -9.6177081e-09 Ad: 1.00e+00 Dobj: -4.9851670e-06 \n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "-3.045195207107554e-10"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "import CSDP\n",
    "solver = CSDP.Optimizer\n",
    "using SumOfSquares\n",
    "function sos_min(sparsity)\n",
    "    model = Model(solver)\n",
    "    @variable(model, t)\n",
    "    @objective(model, Max, t)\n",
    "    con_ref = @constraint(model, poly - t in SOSCone(), sparsity = sparsity)\n",
    "    optimize!(model)\n",
    "    return value(t), moment_matrix(con_ref)\n",
    "end\n",
    "\n",
    "bound, ν = sos_min(Sparsity.NoPattern())\n",
    "bound"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "We find the corresponding minimizer `(0, 0, 0)` by matching the moments\n",
    "of the moment matrix with a dirac measure centered at this minimizer."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Atomic measure on the variables x[1], x[2], x[3] with 1 atoms:\n at [2.0703595321702734e-6, 1.2367907094787687e-6, 2.2362397702348983e-8] with weight 1.0000000014903554"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "extractatoms(ν, 1e-6)"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can see below that the basis contained 6 monomials hence we needed to use 6x6 PSD matrix variables."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₁x₂, x₂x₃, x₁, x₂, x₃, 1])"
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "cell_type": "code",
   "source": [
    "ν.basis"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "Using the monomial/term sparsity method of [WML20a] based on cluster completion, we find the same bound."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter: 14 Ap: 9.68e-01 Pobj: -3.0451952e-10 Ad: 9.70e-01 Dobj: -2.0275242e-07 \n",
      "Success: SDP solved\n",
      "Primal objective value: -3.0451952e-10 \n",
      "Dual objective value: -2.0275242e-07 \n",
      "Relative primal infeasibility: 1.62e-13 \n",
      "Relative dual infeasibility: 1.55e-09 \n",
      "Real Relative Gap: -2.02e-07 \n",
      "XZ Relative Gap: 2.86e-09 \n",
      "DIMACS error measures: 1.75e-13 0.00e+00 2.56e-09 0.00e+00 -2.02e-07 2.86e-09\n",
      "CSDP 6.2.0\n",
      "Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00 \n",
      "Iter:  1 Ap: 9.71e-01 Pobj: -2.4909930e+02 Ad: 8.28e-01 Dobj:  7.2430886e+02 \n",
      "Iter:  2 Ap: 9.99e-01 Pobj: -1.6816440e+02 Ad: 9.32e-01 Dobj:  1.6018610e+02 \n",
      "Iter:  3 Ap: 9.96e-01 Pobj: -6.7498290e+00 Ad: 9.24e-01 Dobj:  6.1399335e+01 \n",
      "Iter:  4 Ap: 1.00e+00 Pobj: -1.4722921e+00 Ad: 9.01e-01 Dobj:  8.9329470e+00 \n",
      "Iter:  5 Ap: 1.00e+00 Pobj: -5.0478524e-01 Ad: 7.90e-01 Dobj:  2.7911688e+00 \n",
      "Iter:  6 Ap: 1.00e+00 Pobj: -1.5409216e-01 Ad: 7.76e-01 Dobj:  8.7431136e-01 \n",
      "Iter:  7 Ap: 1.00e+00 Pobj: -4.8559418e-02 Ad: 8.12e-01 Dobj:  2.4102522e-01 \n",
      "Iter:  8 Ap: 1.00e+00 Pobj: -1.2896403e-02 Ad: 8.39e-01 Dobj:  6.6196298e-02 \n",
      "Iter:  9 Ap: 1.00e+00 Pobj: -1.7033049e-03 Ad: 1.00e+00 Dobj:  1.0847092e-02 \n",
      "Iter: 10 Ap: 1.00e+00 Pobj: -2.4388656e-04 Ad: 1.00e+00 Dobj:  7.6100068e-04 \n",
      "Iter: 11 Ap: 9.96e-01 Pobj: -7.3346824e-06 Ad: 1.00e+00 Dobj: -4.5795945e-05 \n",
      "Iter: 12 Ap: 9.97e-01 Pobj: -2.2313479e-07 Ad: 1.00e+00 Dobj: -7.5985091e-05 \n",
      "Iter: 13 Ap: 1.00e+00 Pobj: -9.6177081e-09 Ad: 1.00e+00 Dobj: -4.9851670e-06 \n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "-3.045195207107554e-10"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "bound, ν = sos_min(Sparsity.Monomial())\n",
    "bound"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "Which is not suprising as no sparsity reduction could be performed."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "1-element Vector{MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}}:\n MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₁x₂, x₂x₃, x₁, x₂, x₃, 1])"
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "[sub.basis for sub in ν.sub_moment_matrices]"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "cell_type": "markdown",
   "source": [
    "Using the monomial/term sparsity method of [WML20b] based on chordal completion, the lower bound is smaller than 0."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter: 14 Ap: 9.68e-01 Pobj: -3.0451952e-10 Ad: 9.70e-01 Dobj: -2.0275242e-07 \n",
      "Success: SDP solved\n",
      "Primal objective value: -3.0451952e-10 \n",
      "Dual objective value: -2.0275242e-07 \n",
      "Relative primal infeasibility: 1.62e-13 \n",
      "Relative dual infeasibility: 1.55e-09 \n",
      "Real Relative Gap: -2.02e-07 \n",
      "XZ Relative Gap: 2.86e-09 \n",
      "DIMACS error measures: 1.75e-13 0.00e+00 2.56e-09 0.00e+00 -2.02e-07 2.86e-09\n",
      "CSDP 6.2.0\n",
      "Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00 \n",
      "Iter:  1 Ap: 9.17e-01 Pobj: -3.5364812e+02 Ad: 6.71e-01 Dobj:  5.8532180e+02 \n",
      "Iter:  2 Ap: 8.26e-01 Pobj: -3.6133495e+02 Ad: 9.36e-01 Dobj:  1.9736716e+02 \n",
      "Iter:  3 Ap: 9.70e-01 Pobj: -1.1800069e+02 Ad: 8.23e-01 Dobj:  1.4747690e+02 \n",
      "Iter:  4 Ap: 1.00e+00 Pobj: -1.2898462e+01 Ad: 8.32e-01 Dobj:  5.5286769e+01 \n",
      "Iter:  5 Ap: 1.00e+00 Pobj: -2.4900158e+00 Ad: 9.07e-01 Dobj:  8.7767866e+00 \n",
      "Iter:  6 Ap: 1.00e+00 Pobj: -1.2196157e+00 Ad: 8.79e-01 Dobj:  2.3717906e+00 \n",
      "Iter:  7 Ap: 1.00e+00 Pobj: -4.3509371e-01 Ad: 6.82e-01 Dobj:  1.1617211e+00 \n",
      "Iter:  8 Ap: 1.00e+00 Pobj: -1.7554721e-01 Ad: 7.67e-01 Dobj:  3.9978327e-01 \n",
      "Iter:  9 Ap: 1.00e+00 Pobj: -6.1256598e-02 Ad: 8.61e-01 Dobj:  9.4968927e-02 \n",
      "Iter: 10 Ap: 1.00e+00 Pobj: -2.0602943e-02 Ad: 1.00e+00 Dobj:  8.3003823e-03 \n",
      "Iter: 11 Ap: 1.00e+00 Pobj: -8.1086444e-03 Ad: 1.00e+00 Dobj:  2.1149126e-04 \n",
      "Iter: 12 Ap: 1.00e+00 Pobj: -5.0556265e-03 Ad: 1.00e+00 Dobj: -2.8106237e-03 \n",
      "Iter: 13 Ap: 1.00e+00 Pobj: -3.7777679e-03 Ad: 1.00e+00 Dobj: -3.4791027e-03 \n",
      "Iter: 14 Ap: 1.00e+00 Pobj: -3.5703395e-03 Ad: 1.00e+00 Dobj: -3.5995200e-03 \n",
      "Iter: 15 Ap: 1.00e+00 Pobj: -3.5520509e-03 Ad: 1.00e+00 Dobj: -3.6089795e-03 \n",
      "Iter: 16 Ap: 1.00e+00 Pobj: -3.5512394e-03 Ad: 1.00e+00 Dobj: -3.5545003e-03 \n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "-0.0035512009904666852"
     },
     "metadata": {},
     "execution_count": 8
    }
   ],
   "cell_type": "code",
   "source": [
    "bound, ν = sos_min(Sparsity.Monomial(ChordalCompletion()))\n",
    "bound"
   ],
   "metadata": {},
   "execution_count": 8
  },
  {
   "cell_type": "markdown",
   "source": [
    "However, this bound was obtained with an SDP with 4 matrices of size 3x3."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "4-element Vector{MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}}:\n MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₂x₃, x₂, x₃])\n MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₂x₃, x₂, 1])\n MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₁, x₂, 1])\n MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x₁x₂, x₁, 1])"
     },
     "metadata": {},
     "execution_count": 9
    }
   ],
   "cell_type": "code",
   "source": [
    "[sub.basis for sub in ν.sub_moment_matrices]"
   ],
   "metadata": {},
   "execution_count": 9
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.2"
  },
  "kernelspec": {
   "name": "julia-1.7",
   "display_name": "Julia 1.7.2",
   "language": "julia"
  }
 },
 "nbformat": 4
}