{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Minimization of a univariate polynomial"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Contributed by**: Benoît Legat"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using DynamicPolynomials\n",
    "using SumOfSquares\n",
    "import CSDP"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "Consider the problem of finding both the minimum value of `p = x^4 - 4x^3 - 2x^2 + 12x + 3` as well as its minimizers."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can use SumOfSquares.jl to find such these values as follows.\n",
    "We first define the polynomial using DynamicPolynomials."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3 + 12x - 2x² - 4x³ + x⁴",
      "text/latex": "$$ 3 + 12x - 2x^{2} - 4x^{3} + x^{4} $$"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "@polyvar x\n",
    "p = x^4 - 4x^3 - 2x^2 + 12x + 3"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "Secondly, we create a Sum-of-Squares program searching for the maximal lower bound `σ` of the polynomial."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "σ",
      "text/latex": "$ σ $"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "model = SOSModel(CSDP.Optimizer)\n",
    "@variable(model, σ)\n",
    "@constraint(model, cref, p >= σ)\n",
    "@objective(model, Max, σ)"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "Thirdly, solve the program and find `σ = -6` as lower bound:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success: SDP solved\n",
      "Primal objective value: 0.0000000e+00 \n",
      "Dual objective value: 0.0000000e+00 \n",
      "Relative primal infeasibility: 3.29e-17 \n",
      "Relative dual infeasibility: 5.00e-11 \n",
      "Real Relative Gap: 0.00e+00 \n",
      "XZ Relative Gap: 1.06e-10 \n",
      "DIMACS error measures: 4.93e-17 0.00e+00 5.00e-11 0.00e+00 0.00e+00 1.06e-10\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: 8.10e-01 Pobj: -1.3830330e+01 Ad: 7.29e-01 Dobj: -8.2765454e+00 \n",
      "Iter:  2 Ap: 9.00e-01 Pobj: -1.5490537e+01 Ad: 9.00e-01 Dobj: -1.3831579e-01 \n",
      "Iter:  3 Ap: 1.00e+00 Pobj: -7.5281681e+00 Ad: 9.00e-01 Dobj: -3.8241528e+00 \n",
      "Iter:  4 Ap: 1.00e+00 Pobj: -6.0076491e+00 Ad: 9.00e-01 Dobj: -5.7732642e+00 \n",
      "Iter:  5 Ap: 1.00e+00 Pobj: -6.0000210e+00 Ad: 1.00e+00 Dobj: -5.9999957e+00 \n",
      "Iter:  6 Ap: 1.00e+00 Pobj: -6.0000045e+00 Ad: 1.00e+00 Dobj: -6.0000015e+00 \n",
      "Iter:  7 Ap: 1.00e+00 Pobj: -6.0000003e+00 Ad: 1.00e+00 Dobj: -6.0000001e+00 \n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "* Solver : CSDP\n\n* Status\n  Result count       : 1\n  Termination status : OPTIMAL\n  Message from the solver:\n  \"Problem solved to optimality.\"\n\n* Candidate solution (result #1)\n  Primal status      : FEASIBLE_POINT\n  Dual status        : FEASIBLE_POINT\n  Objective value    : -6.00000e+00\n  Dual objective value : -6.00000e+00\n\n* Work counters\n  Solve time (sec)   : 9.92060e-04\n"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "optimize!(model)\n",
    "solution_summary(model)"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can look at the certificate that `σ = -6` is a lower bound:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(-3.000000004964362 - 1.999999998797403*x + 0.9999999995176428*x^2)^2"
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "cell_type": "code",
   "source": [
    "sos_dec = sos_decomposition(cref, 1e-4)"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "Indeed, `p + 6 = (x^2 - 2x - 3)^2` so `p ≥ -6`.\n",
    "\n",
    "## Extraction of minimizers\n",
    "\n",
    "We can now find the minimizers from the moment matrix:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3×3 SymMatrix{Float64}:\n 1.0        0.0669168   3.13383\n 0.0669168  3.13383     6.46842\n 3.13383    6.46842    22.3383"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "ν = moment_matrix(cref)\n",
    "ν.Q"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "This matrix is the convex combination of the moment matrices corresponding to two atomic measures at `-1` and `3`\n",
    "which allows us to conclude that `-1` and `3` are global minimizers."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "2-element Vector{Float64}:\n -1.0000001975826238\n  2.9999999542477047"
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "η = atomic_measure(ν, 1e-4)\n",
    "minimizers = [η.atoms[1].center; η.atoms[2].center]"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "cell_type": "markdown",
   "source": [
    "Below are more details on what we mean by convex combination.\n",
    "The moment matrix of the atomic measure at the first minimizer is:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3×3 SymMatrix{Float64}:\n  1.0  -1.0   1.0\n -1.0   1.0  -1.0\n  1.0  -1.0   1.0"
     },
     "metadata": {},
     "execution_count": 8
    }
   ],
   "cell_type": "code",
   "source": [
    "η1 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[1])), ν.basis.monomials)\n",
    "η1.Q"
   ],
   "metadata": {},
   "execution_count": 8
  },
  {
   "cell_type": "markdown",
   "source": [
    "The moment matrix of the atomic measure at the second minimizer is:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3×3 SymMatrix{Float64}:\n 1.0   3.0   9.0\n 3.0   9.0  27.0\n 9.0  27.0  81.0"
     },
     "metadata": {},
     "execution_count": 9
    }
   ],
   "cell_type": "code",
   "source": [
    "η2 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[2])), ν.basis.monomials)\n",
    "η2.Q"
   ],
   "metadata": {},
   "execution_count": 9
  },
  {
   "cell_type": "markdown",
   "source": [
    "And the moment matrix is the convex combination of both:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3×3 Matrix{Float64}:\n 1.0        0.0669169   3.13383\n 0.0669169  3.13383     6.46842\n 3.13383    6.46842    22.3383"
     },
     "metadata": {},
     "execution_count": 10
    }
   ],
   "cell_type": "code",
   "source": [
    "Q12 = η1.Q * η.atoms[1].weight + η2.Q * η.atoms[2].weight"
   ],
   "metadata": {},
   "execution_count": 10
  },
  {
   "cell_type": "markdown",
   "source": [
    "Another way to see this (by linearity of the expectation) is that `ν` is the moment matrix\n",
    "of the convex combination of the two atomic measures."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Changing the polynomial basis\n",
    "\n",
    "The monomial basis used by default can leave a problem quite ill-conditioned for the solver.\n",
    "Let's try to use another basis instead:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter:  8 Ap: 9.00e-01 Pobj: -6.0000000e+00 Ad: 1.00e+00 Dobj: -6.0000000e+00 \n",
      "Success: SDP solved\n",
      "Primal objective value: -6.0000000e+00 \n",
      "Dual objective value: -6.0000000e+00 \n",
      "Relative primal infeasibility: 6.64e-12 \n",
      "Relative dual infeasibility: 3.14e-10 \n",
      "Real Relative Gap: 2.23e-09 \n",
      "XZ Relative Gap: 3.25e-09 \n",
      "DIMACS error measures: 7.25e-12 0.00e+00 5.34e-10 0.00e+00 2.23e-09 3.25e-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.00e-01 Pobj: -1.8383751e+01 Ad: 7.29e-01 Dobj: -1.0343874e+01 \n",
      "Iter:  2 Ap: 9.00e-01 Pobj: -1.5927039e+01 Ad: 9.00e-01 Dobj: -2.3550363e+00 \n",
      "Iter:  3 Ap: 9.00e-01 Pobj: -7.7938220e+00 Ad: 9.00e-01 Dobj: -3.9182690e+00 \n",
      "Iter:  4 Ap: 9.00e-01 Pobj: -6.1761570e+00 Ad: 9.00e-01 Dobj: -5.7791841e+00 \n",
      "Iter:  5 Ap: 9.00e-01 Pobj: -6.0169870e+00 Ad: 9.00e-01 Dobj: -5.9771917e+00 \n",
      "Iter:  6 Ap: 9.00e-01 Pobj: -6.0016653e+00 Ad: 1.00e+00 Dobj: -5.9999542e+00 \n",
      "Iter:  7 Ap: 1.00e+00 Pobj: -6.0002842e+00 Ad: 1.00e+00 Dobj: -5.9997933e+00 \n",
      "Iter:  8 Ap: 9.00e-01 Pobj: -6.0000531e+00 Ad: 1.00e+00 Dobj: -5.9999938e+00 \n",
      "Iter:  9 Ap: 1.00e+00 Pobj: -6.0000089e+00 Ad: 1.00e+00 Dobj: -5.9999987e+00 \n",
      "Iter: 10 Ap: 1.00e+00 Pobj: -6.0000012e+00 Ad: 1.00e+00 Dobj: -6.0000000e+00 \n",
      "Iter: 11 Ap: 9.00e-01 Pobj: -6.0000002e+00 Ad: 1.00e+00 Dobj: -6.0000000e+00 \n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "* Solver : CSDP\n\n* Status\n  Result count       : 1\n  Termination status : OPTIMAL\n  Message from the solver:\n  \"Problem solved to optimality.\"\n\n* Candidate solution (result #1)\n  Primal status      : FEASIBLE_POINT\n  Dual status        : FEASIBLE_POINT\n  Objective value    : -6.00000e+00\n  Dual objective value : -6.00000e+00\n\n* Work counters\n  Solve time (sec)   : 1.50084e-03\n"
     },
     "metadata": {},
     "execution_count": 11
    }
   ],
   "cell_type": "code",
   "source": [
    "model = SOSModel(CSDP.Optimizer)\n",
    "@variable(model, σ)\n",
    "@constraint(model, cheby_cref, p >= σ, basis = ChebyshevBasisFirstKind)\n",
    "@objective(model, Max, σ)\n",
    "optimize!(model)\n",
    "solution_summary(model)"
   ],
   "metadata": {},
   "execution_count": 11
  },
  {
   "cell_type": "markdown",
   "source": [
    "Although the gram matrix in the monomial basis:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g.basis = MonomialBasis([1, x, x²])\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3×3 SymMatrix{Float64}:\n  9.0   6.0  -3.0\n  6.0   4.0  -2.0\n -3.0  -2.0   1.0"
     },
     "metadata": {},
     "execution_count": 12
    }
   ],
   "cell_type": "code",
   "source": [
    "g = gram_matrix(cref)\n",
    "@show g.basis\n",
    "g.Q"
   ],
   "metadata": {},
   "execution_count": 12
  },
  {
   "cell_type": "markdown",
   "source": [
    "looks different from the gram matrix in the Chebyshev basis:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cheby_g.basis = ChebyshevBasisFirstKind([1.0, x, -1.0 + 2.0x²])\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3×3 SymMatrix{Float64}:\n  6.25   5.0  -1.25\n  5.0    4.0  -1.0\n -1.25  -1.0   0.25"
     },
     "metadata": {},
     "execution_count": 13
    }
   ],
   "cell_type": "code",
   "source": [
    "cheby_g = gram_matrix(cheby_cref)\n",
    "@show cheby_g.basis\n",
    "cheby_g.Q"
   ],
   "metadata": {},
   "execution_count": 13
  },
  {
   "cell_type": "markdown",
   "source": [
    "they both yields the same Sum-of-Squares decomposition:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(-3.000000002974674 - 1.9999999995036628*x + 0.9999999997868867*x^2)^2"
     },
     "metadata": {},
     "execution_count": 14
    }
   ],
   "cell_type": "code",
   "source": [
    "cheby_sos_dec = sos_decomposition(cheby_cref, 1e-4)"
   ],
   "metadata": {},
   "execution_count": 14
  },
  {
   "cell_type": "markdown",
   "source": [
    "The gram matrix in the Chebyshev basis can be understood as follows.\n",
    "To express the polynomial $-x^2 + 2x + 3$ in the Chebyshev basis, we start by\n",
    " substituting $x$ into $\\cos(\\theta)$ to obtain\n",
    "$-\\cos(\\theta)^2 + 2\\cos(\\theta) + 3$.\n",
    "We now express it as a combination of $\\cos(n\\theta)$ for $n = 0, 1, 2$:\n",
    "$-(2\\cos(\\theta) - 1) /2 + 2 \\cos(\\theta) + 5/2.$\n",
    "Therefore, the coefficients in the Chebyshev basis is:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3-element Vector{Float64}:\n  2.5\n  2.0\n -0.5"
     },
     "metadata": {},
     "execution_count": 15
    }
   ],
   "cell_type": "code",
   "source": [
    "cheby_coefs = [5/2, 2, -1/2]"
   ],
   "metadata": {},
   "execution_count": 15
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can indeed observe that we obtain the same matrix as `cheby_g.Q`"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3×3 Matrix{Float64}:\n  6.25   5.0  -1.25\n  5.0    4.0  -1.0\n -1.25  -1.0   0.25"
     },
     "metadata": {},
     "execution_count": 16
    }
   ],
   "cell_type": "code",
   "source": [
    "cheby_coefs * cheby_coefs'"
   ],
   "metadata": {},
   "execution_count": 16
  },
  {
   "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.10.3"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}