{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Noncommutative variables\n",
    "\n",
    "## Example 2.11\n",
    "\n",
    "We consider the Example 2.11 of [BKP16] in which the polynomial with noncommutative variables\n",
    "$(x * y + x^2)^2 = x^4 + x^3y + xyx^2 + xyxy$ is tested to be sum-of-squares.\n",
    "\n",
    "[BKP16] Sabine Burgdorf, Igor Klep, and Janez Povh.\n",
    "*Optimization of polynomials in non-commuting variables*.\n",
    "Berlin: Springer, 2016."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ x^{4} + x^{3}y + xyx^{2} + xyxy $$"
      ],
      "text/plain": [
       "x⁴ + x³y + xyx² + xyxy"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using DynamicPolynomials\n",
    "@ncpolyvar x y\n",
    "p = (x * y + x^2)^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "MathOptInterface.OptimizerWithAttributes(CSDP.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute,Any}[MathOptInterface.Silent() => true])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import CSDP\n",
    "using JuMP\n",
    "optimizer_constructor = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Newton polytope method has not been adapted to the noncommutative case yet,\n",
    "so we force the choice of certificate to `MaxDegree` instead of `Newton`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$ (1)x^{4} + (1)x^{3}y + (1)xyx^{2} + (1)xyxy \\text{ is SOS} $"
      ],
      "text/plain": [
       "(1)x⁴ + (1)x³y + (1)xyx² + (1)xyxy is SOS"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using SumOfSquares\n",
    "model = Model(optimizer_constructor)\n",
    "con_ref = @constraint(model, p in SOSCone())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "optimize!(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that both the monomials `xy` and `yx` are considered separately, this is a difference with the commutative version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "MonomialBasis{Monomial{false},MonomialVector{false}}(Monomial{false}[x², xy])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "certificate_basis(con_ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the solution correctly uses the monomial `xy` instead of `yx`. We also identify that only the monomials `x^2` and `xy` would be needed. This would be dectected by the Newton chip method of [Section 2.3, BKP16]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 SymMatrix{Float64}:\n",
       " 1.0  1.0\n",
       " 1.0  1.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gram_matrix(con_ref).Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When asking for the SOS decomposition, the numerically small entries makes the solution less readable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.0000000000000002*x^2 - x*y)^2 + (-6.265166515512128e-9*x^2 + 6.265166515512127e-9*x*y)^2"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sos_decomposition(con_ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They are however easily discarded by using a nonzero tolerance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.0000000000000002*x^2 - x*y)^2"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sos_decomposition(con_ref, 1e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2.2\n",
    "\n",
    "We consider now the Example 2.2 of [BKP16] in which the polynomial with noncommutative variables\n",
    "$(x + x^{10}y^{20}x^{10})^2$ is tested to be sum-of-squares."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ x^{10}y^{20}x^{20}y^{20}x^{10} + x^{11}y^{20}x^{10} + x^{10}y^{20}x^{11} + x^{2} $$"
      ],
      "text/plain": [
       "x¹⁰y²⁰x²⁰y²⁰x¹⁰ + x¹¹y²⁰x¹⁰ + x¹⁰y²⁰x¹¹ + x²"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using DynamicPolynomials\n",
    "@ncpolyvar x y\n",
    "n = 10\n",
    "p = (x + x^n * y^(2n) * x^n)^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$ (1)x^{10}y^{20}x^{20}y^{20}x^{10} + (1)x^{11}y^{20}x^{10} + (1)x^{10}y^{20}x^{11} + (1)x^{2} \\text{ is SOS} $"
      ],
      "text/plain": [
       "(1)x¹⁰y²⁰x²⁰y²⁰x¹⁰ + (1)x¹¹y²⁰x¹⁰ + (1)x¹⁰y²⁰x¹¹ + (1)x² is SOS"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using SumOfSquares\n",
    "model = Model(optimizer_constructor)\n",
    "con_ref = @constraint(model, p in SOSCone())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "optimize!(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only two monomials were considered for the basis of the gram matrix thanks to the Augmented Newton chip method detailed in [Section 2.4, BKP16]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "MonomialBasis{Monomial{false},MonomialVector{false}}(Monomial{false}[x¹⁰y²⁰x¹⁰, x])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "certificate_basis(con_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 SymMatrix{Float64}:\n",
       " 1.0  1.0\n",
       " 1.0  1.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gram_matrix(con_ref).Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.0000000000000002*x^10*y^20*x^10 - x)^2"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sos_decomposition(con_ref, 1e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.3.1",
   "language": "julia",
   "name": "julia-1.3"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.3.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}