{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "We want to compute $A_4^*$. In other words, we are looking for the coefficients of a convex polynomial\n", "$$p(x, y) = \\sum_{i=0}^{2d} c_i x^i y^{2d-i}$$\n", "that maximizes\n", "$$\\frac{c_4}{70}.$$\n", "and satisfies\n", "$p(1, 0) + p(0, 1) = 2$, or equivalently $c_0 + c_8 = 1$.\n", "\n", "In the following we perform a series of symmetry reduction techniques that will allow us to solve this problem analytically.\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "code_folding": [ 3, 17 ] }, "outputs": [], "source": [ "# Some helpful functions\n", "\n", "# The following function creates an N x N symmetric symbolic matrix\n", "def symb_matrix(N, base_name='Q_%d%d', sym=False):\n", " var_swapper = lambda i,j: (i, j)\n", " if sym:\n", " var_swapper = lambda i,j: (min(i,j), max(i,j))\n", " \n", " Q_varnames = ','.join([base_name % var_swapper( i,j) for (i,j) in cartesian_product([range(N),range(N)])])\n", " Q_vars = var(Q_varnames) \n", " Q = matrix(SR,N,N,Q_vars)\n", " \n", " return Q\n", "\n", "\n", "# the following function gives a full rank description of\n", "# matrices `Q` that satisfy the equations in `eqn`\n", "def eliminate_var_matrix(Q, eqn, name, additional_vars=[]):\n", " all_vars = list(Q.variables())\n", " all_vars += additional_vars\n", " all_vars = map(SR, all_vars)\n", "\n", " sols = solve(map(lambda eqn_i: SR(eqn_i) == 0, eqn), all_vars)\n", " sols = sols[0]\n", " \n", " # write Q as sum ri * Fi, where Fi are constant matrices\n", " r_F = [(term.rhs(), jacobian(Q, term.lhs())[0,0]) for term in sols]\n", " Fr = sum([ri * Fi for ri, Fi in r_F])\n", " \n", " # substitue the variables ri with their full rank descirption\n", " var_to_sub = sum([list(term.rhs().variables()) for term in sols], [])\n", " var_sub = {v: var(name+str(i)) for i,v in enumerate(var_to_sub)}\n", " return Fr.subs(var_sub)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reducing the number of coefficients in $p$\n", "\n", "Note that if $p$ is such a polynomial, then so is $p(-x, -y)$, so we\n", "can assume that $p$ is even (i.e., $c_i = 0$ for $i=1,3,5,7$). Furthermore, if $p$ is a solution, so is\n", "$p(y, x)$, so we can assume that $p$ is symmetric (i.e., $c_i = c_{8-i}$). All in all, without\n", "loss of generality we assume\n", "\n", "$$p(x, y) = x^8 + y^8 + \\alpha x^4y^4 + \\beta (x^2y^6 + x^6y^2) \\text{ for some $\\alpha, \\beta \\in \\mathbb R$}.$$\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}b x^{6} y^{2} + a x^{4} y^{4} + b x^{2} y^{6} + x^{8} + y^{8}$ " ], "text/plain": [ "b*x^6*y^2 + a*x^4*y^4 + b*x^2*y^6 + x^8 + y^8" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%display latex\n", "d = 4\n", "R. = QQ['a,b,x,y,u,v']\n", "vars = [x,y,u,v]\n", "q = x^8 + y^8 + a*x^4*y^4 + b*(x^2*y^6 + x^6*y^2)\n", "q" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Formulating the problem as an SDP\n", "Because of Theorem 3.6, imposing the condition on the convexity of $p$ we can formulated as a search problem for an $8 \\times 8$ positive semidefinite matrix $Q$ that satisfies\n", " $${\\mathbf u}^T\\nabla^2q({\\mathbf x}){\\mathbf u} = z^TQz \\quad \\forall {\\mathbf x},{\\mathbf u} \\in \\mathbb R^2$$\n", " where\n", "$$z^T({\\mathbf x}, {\\mathbf u}) = (u x^{2d}, ux^{2d-1}y, \\ldots, v y^{2d}, v x^{2d}, v x^{2d-1}y, \\ldots, v y^{2d}).$$\n", "\n", "Once we have $Q$, we can recover the polynomial $q$ via Euler's identity as follows:\n", "\n", "$$q(\\mathbf x) = \\frac1{8 \\cdot 7} z(\\mathbf x, \\mathbf x)^T Q z(\\mathbf x, \\mathbf x)$$\n", "\n", "\n", "Therefore, the objective function ca be written exclusively in terms of the matrix $Q$, and we can rewrite our optimization problem (*) completely in terms of the matrix $Q$ as follows:\n", "\n", "$$\\max_{Q \\succeq 0} Q \\text{ s.t. } Q \\in \\mathcal L$$\n", "where $\\mathcal L$ is the following linear space\n", "$$\\{Q \\; | \\; \\exists q = x^8 + y^8 + \\alpha x^4y^4 + \\beta (x^2y^6 + x^6y^2) \\text{ s.t. } {\\mathbf u}^T\\nabla^2q({\\mathbf x}){\\mathbf u} = z^TQz.\\}$$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "code_folding": [] }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrrrrrr}\n", "Q_{00} & Q_{01} & Q_{02} & Q_{03} & Q_{04} & Q_{05} & Q_{06} & Q_{07} \\\\\n", "Q_{01} & Q_{11} & Q_{12} & Q_{13} & Q_{14} & Q_{15} & Q_{16} & Q_{17} \\\\\n", "Q_{02} & Q_{12} & Q_{22} & Q_{23} & Q_{24} & Q_{25} & Q_{26} & Q_{27} \\\\\n", "Q_{03} & Q_{13} & Q_{23} & Q_{33} & Q_{34} & Q_{35} & Q_{36} & Q_{37} \\\\\n", "Q_{04} & Q_{14} & Q_{24} & Q_{34} & Q_{44} & Q_{45} & Q_{46} & Q_{47} \\\\\n", "Q_{05} & Q_{15} & Q_{25} & Q_{35} & Q_{45} & Q_{55} & Q_{56} & Q_{57} \\\\\n", "Q_{06} & Q_{16} & Q_{26} & Q_{36} & Q_{46} & Q_{56} & Q_{66} & Q_{67} \\\\\n", "Q_{07} & Q_{17} & Q_{27} & Q_{37} & Q_{47} & Q_{57} & Q_{67} & Q_{77}\n", "\\end{array}\\right)$ " ], "text/plain": [ "[Q_00 Q_01 Q_02 Q_03 Q_04 Q_05 Q_06 Q_07]\n", "[Q_01 Q_11 Q_12 Q_13 Q_14 Q_15 Q_16 Q_17]\n", "[Q_02 Q_12 Q_22 Q_23 Q_24 Q_25 Q_26 Q_27]\n", "[Q_03 Q_13 Q_23 Q_33 Q_34 Q_35 Q_36 Q_37]\n", "[Q_04 Q_14 Q_24 Q_34 Q_44 Q_45 Q_46 Q_47]\n", "[Q_05 Q_15 Q_25 Q_35 Q_45 Q_55 Q_56 Q_57]\n", "[Q_06 Q_16 Q_26 Q_36 Q_46 Q_56 Q_66 Q_67]\n", "[Q_07 Q_17 Q_27 Q_37 Q_47 Q_57 Q_67 Q_77]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "Q = symb_matrix(8, 'Q_%d%d', sym=True)\n", "Q" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "code_folding": [] }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrrrrrr}\n", "-w_{53} + 2 \\, w_{62} - 2 \\, w_{63} - 28 \\, w_{68} + w_{71} + w_{73} - 2 \\, w_{74} + 28 & w_{62} - w_{63} - 12 \\, w_{68} + w_{72} & 0 & w_{65} - w_{67} + w_{70} & -\\frac{1}{2} \\, w_{48} + \\frac{1}{2} \\, w_{56} + w_{57} & -\\frac{1}{2} \\, w_{48} - \\frac{1}{2} \\, w_{56} & \\frac{1}{2} \\, w_{44} - \\frac{1}{2} \\, w_{54} + w_{55} & -\\frac{1}{2} \\, w_{44} - \\frac{1}{2} \\, w_{54} \\\\\n", "w_{62} - w_{63} - 12 \\, w_{68} + w_{72} & w_{53} & w_{65} & -w_{67} - w_{70} & w_{48} & -w_{57} & w_{44} & -w_{55} \\\\\n", "0 & w_{65} & w_{53} - 2 \\, w_{62} + 2 \\, w_{63} + 28 \\, w_{68} - w_{71} - w_{73} + 2 \\, w_{74} + 28 & w_{62} & -\\frac{1}{2} \\, w_{44} + \\frac{1}{2} \\, w_{54} - w_{55} & -\\frac{1}{2} \\, w_{44} - \\frac{1}{2} \\, w_{54} & \\frac{1}{2} \\, w_{48} - \\frac{1}{2} \\, w_{56} - w_{57} & -\\frac{1}{2} \\, w_{48} - \\frac{1}{2} \\, w_{56} \\\\\n", "w_{65} - w_{67} + w_{70} & -w_{67} - w_{70} & w_{62} & -w_{53} - 2 \\, w_{63} - 2 \\, w_{72} + 12 \\, w_{74} & w_{54} & w_{55} & w_{56} & w_{57} \\\\\n", "-\\frac{1}{2} \\, w_{48} + \\frac{1}{2} \\, w_{56} + w_{57} & w_{48} & -\\frac{1}{2} \\, w_{44} + \\frac{1}{2} \\, w_{54} - w_{55} & w_{54} & -4 \\, w_{62} + 2 \\, w_{63} + 54 \\, w_{68} - w_{71} - 2 \\, w_{72} & w_{63} & -2 \\, w_{65} + w_{67} - w_{70} & w_{67} \\\\\n", "-\\frac{1}{2} \\, w_{48} - \\frac{1}{2} \\, w_{56} & -w_{57} & -\\frac{1}{2} \\, w_{44} - \\frac{1}{2} \\, w_{54} & w_{55} & w_{63} & 2 \\, w_{68} - w_{73} & w_{70} & 0 \\\\\n", "\\frac{1}{2} \\, w_{44} - \\frac{1}{2} \\, w_{54} + w_{55} & w_{44} & \\frac{1}{2} \\, w_{48} - \\frac{1}{2} \\, w_{56} - w_{57} & w_{56} & -2 \\, w_{65} + w_{67} - w_{70} & w_{70} & w_{71} & w_{72} \\\\\n", "-\\frac{1}{2} \\, w_{44} - \\frac{1}{2} \\, w_{54} & -w_{55} & -\\frac{1}{2} \\, w_{48} - \\frac{1}{2} \\, w_{56} & w_{57} & w_{67} & 0 & w_{72} & w_{73}\n", "\\end{array}\\right)$ " ], "text/plain": [ "[-w53 + 2*w62 - 2*w63 - 28*w68 + w71 + w73 - 2*w74 + 28 w62 - w63 - 12*w68 + w72 0 w65 - w67 + w70 -1/2*w48 + 1/2*w56 + w57 -1/2*w48 - 1/2*w56 1/2*w44 - 1/2*w54 + w55 -1/2*w44 - 1/2*w54]\n", "[ w62 - w63 - 12*w68 + w72 w53 w65 -w67 - w70 w48 -w57 w44 -w55]\n", "[ 0 w65 w53 - 2*w62 + 2*w63 + 28*w68 - w71 - w73 + 2*w74 + 28 w62 -1/2*w44 + 1/2*w54 - w55 -1/2*w44 - 1/2*w54 1/2*w48 - 1/2*w56 - w57 -1/2*w48 - 1/2*w56]\n", "[ w65 - w67 + w70 -w67 - w70 w62 -w53 - 2*w63 - 2*w72 + 12*w74 w54 w55 w56 w57]\n", "[ -1/2*w48 + 1/2*w56 + w57 w48 -1/2*w44 + 1/2*w54 - w55 w54 -4*w62 + 2*w63 + 54*w68 - w71 - 2*w72 w63 -2*w65 + w67 - w70 w67]\n", "[ -1/2*w48 - 1/2*w56 -w57 -1/2*w44 - 1/2*w54 w55 w63 2*w68 - w73 w70 0]\n", "[ 1/2*w44 - 1/2*w54 + w55 w44 1/2*w48 - 1/2*w56 - w57 w56 -2*w65 + w67 - w70 w70 w71 w72]\n", "[ -1/2*w44 - 1/2*w54 -w55 -1/2*w48 - 1/2*w56 w57 w67 0 w72 w73]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# vector of monomial\n", "z = vector((x^3*u - y^3*v, x*y^2*u - x^2*y*v,\n", " x^3*u + y^3*v, x*y^2*u + x^2*y*v, \n", " x^2*y*u - x*y^2*v, y^3*u - x^3*v,\n", " x^2*y*u + x*y^2*v, y^3*u + x^3*v))\n", "\n", "# hessian of q\n", "Hq = jacobian(jacobian(q, (x, y)), (x, y))\n", "\n", "# Q is in L if the following polynomial has all its coefficients equal to 0\n", "diff = z * Q * z - vector([u,v]) * Hq * vector([u, v])\n", "eqn_L = QQ[Q.variables() + (a,b)][x,y,u,v](diff).coefficients()\n", "\n", "\n", "Q_L = eliminate_var_matrix(Q, eqn_L, 'w', [a])\n", "Q_L" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{1}{70} w_{74}$ " ], "text/plain": [ "1/70*w74" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# write the objective function in terms of the new paramterization\n", "q_from_Q = (z*Q_L*z).subs({u:x, v:y}) / (2*d*(2*d-1))\n", "q_from_Q = QQ[Q_L.variables()][x,y](q_from_Q)\n", "objective = q_from_Q.coefficient({xi: 4 for xi in q_from_Q.variables()}) / binomial(2*d, d)\n", "objective = objective.base_ring()(objective)\n", "objective" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Symmetry reduction techniques applied to $Q$\n", "\n", "Let $\\rho: \\mathbb R^4 \\mapsto \\mathbb R^4$ be any linear transformation of the variables $(x,y,u,v)$ that leaves the coefficient of $x^4y^4$ in $p$ invariant and leaves polynomial $q(x,y,u,v) := {\\mathbf u}^T\\nabla^2q({\\mathbf x}){\\mathbf u}$ invariant, i.e.\n", "$$q(x,y,u,v) = q \\circ \\rho(x,y,u,v).$$\n", "The function $\\rho$ naturally acts linearly on the vector $z(\\mathbf x, \\mathbf y)$. If we call the $8 \\times 8$ matrix of this linear operator $M_\\rho$, then \n", "\n", "$$Q \\in \\mathcal L \\iff M_\\rho^T Q M_\\rho \\in \\mathcal L.$$\n", " \n", " \n", "Let $G$ be a group of such transformations $\\rho$, then $Q$ is an optimal solution for (*) if and only if\n", "\n", "$$Q^* = \\frac1{|G|} \\sum_{\\rho \\in G} M_\\rho^T Q M_\\rho$$\n", "\n", "is also a solution. \n", "\n", "We will take $G$ to be the group generated by the following transformations\n", "\n", "\\begin{align*}\n", " (x,y,u,v) & \\mapsto (y,x,v,u)\\\\\n", " \\cdot \\quad & \\mapsto (-x,y,-u,v)\\\\\n", " \\cdot \\quad & \\mapsto (-x,-y,u,v)\\\\\n", " \\cdot \\quad & \\mapsto (-x,y,u,-v).\n", "\\end{align*}\n", "\n", "The matrix $Q^*$ will have a very structured sparsity pattern.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To make our computations exact, we work over the ring $R := \\mathbb Q[x,y,u,v]$ of polynomials in the variables $x,y,u,v$ with rational coefficients." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(\\begin{array}{rrrr}\n", "0 & 1 & 0 & 0 \\\\\n", "1 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 1 \\\\\n", "0 & 0 & 1 & 0\n", "\\end{array}\\right), \\left(\\begin{array}{rrrr}\n", "-1 & 0 & 0 & 0 \\\\\n", "0 & 1 & 0 & 0 \\\\\n", "0 & 0 & -1 & 0 \\\\\n", "0 & 0 & 0 & 1\n", "\\end{array}\\right), \\left(\\begin{array}{rrrr}\n", "-1 & 0 & 0 & 0 \\\\\n", "0 & -1 & 0 & 0 \\\\\n", "0 & 0 & 1 & 0 \\\\\n", "0 & 0 & 0 & 1\n", "\\end{array}\\right), \\left(\\begin{array}{rrrr}\n", "-1 & 0 & 0 & 0 \\\\\n", "0 & 1 & 0 & 0 \\\\\n", "0 & 0 & 1 & 0 \\\\\n", "0 & 0 & 0 & -1\n", "\\end{array}\\right), \\left(\\begin{array}{rrrr}\n", "1 & 0 & 0 & 0 \\\\\n", "0 & -1 & 0 & 0 \\\\\n", "0 & 0 & 1 & 0 \\\\\n", "0 & 0 & 0 & -1\n", "\\end{array}\\right), \\left(\\begin{array}{rrrr}\n", "1 & 0 & 0 & 0 \\\\\n", "0 & -1 & 0 & 0 \\\\\n", "0 & 0 & -1 & 0 \\\\\n", "0 & 0 & 0 & 1\n", "\\end{array}\\right)\\right]$ " ], "text/plain": [ "[\n", "[0 1 0 0] [-1 0 0 0] [-1 0 0 0] [-1 0 0 0] [ 1 0 0 0]\n", "[1 0 0 0] [ 0 1 0 0] [ 0 -1 0 0] [ 0 1 0 0] [ 0 -1 0 0]\n", "[0 0 0 1] [ 0 0 -1 0] [ 0 0 1 0] [ 0 0 1 0] [ 0 0 1 0]\n", "[0 0 1 0], [ 0 0 0 1], [ 0 0 0 1], [ 0 0 0 -1], [ 0 0 0 -1],\n", "\n", "[ 1 0 0 0]\n", "[ 0 -1 0 0]\n", "[ 0 0 -1 0]\n", "[ 0 0 0 1]\n", "]" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "symmetries =[[y,x,v,u],\n", " [-x,y,-u,v],\n", " [-x,-y,u,v],\n", " [-x,y,u,-v],\n", " [x, -y, u, -v],\n", " [x, -y, -u, v]]\n", "symmetries_rho = map(lambda sym_i: matrix(QQ, jacobian(sym_i, vars)), symmetries)\n", "symmetries_rho" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The group G has 16 elements\n" ] } ], "source": [ "G = gap.Group(symmetries_rho)\n", "G = gap.Elements(G).sage()\n", "G = [matrix(QQ, g_elem_i) for g_elem_i in G]\n", "print(\"The group G has %d elements\"% len(G))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "First two elements of M_rho\n" ] }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(\\begin{array}{rrrrrrrr}\n", "1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\n", "\\end{array}\\right), \\left(\\begin{array}{rrrrrrrr}\n", "-1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & -1\n", "\\end{array}\\right)\\right]$ " ], "text/plain": [ "[\n", "[1 0 0 0 0 0 0 0] [-1 0 0 0 0 0 0 0]\n", "[0 1 0 0 0 0 0 0] [ 0 -1 0 0 0 0 0 0]\n", "[0 0 1 0 0 0 0 0] [ 0 0 -1 0 0 0 0 0]\n", "[0 0 0 1 0 0 0 0] [ 0 0 0 -1 0 0 0 0]\n", "[0 0 0 0 1 0 0 0] [ 0 0 0 0 -1 0 0 0]\n", "[0 0 0 0 0 1 0 0] [ 0 0 0 0 0 -1 0 0]\n", "[0 0 0 0 0 0 1 0] [ 0 0 0 0 0 0 -1 0]\n", "[0 0 0 0 0 0 0 1], [ 0 0 0 0 0 0 0 -1]\n", "]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# action of G on the vector of monomials z\n", "rho_acting_on_z = [z.subs({vi: subi for vi,subi in zip(vars, g*vector(vars))}) for g in G]\n", "\n", "# matrix representation of that action\n", "def compute_M_rho(rho_z):\n", " return matrix(QQ,\n", " [ [1 if zj == zi else -1 if zj == -zi else 0 for zi in z]\n", " for zj in rho_z])\n", "\n", "M_rho = map(compute_M_rho, rho_acting_on_z)\n", "print(\"First two elements of M_rho\")\n", "M_rho[:2]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Blocks on the diagonal of Q:\n" ] }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(\\begin{array}{rr}\n", "-w_{53} + 2 \\, w_{62} - 2 \\, w_{63} - 28 \\, w_{68} + w_{71} + w_{73} - 2 \\, w_{74} + 28 & w_{62} - w_{63} - 12 \\, w_{68} + w_{72} \\\\\n", "w_{62} - w_{63} - 12 \\, w_{68} + w_{72} & w_{53}\n", "\\end{array}\\right), \\left(\\begin{array}{rr}\n", "w_{53} - 2 \\, w_{62} + 2 \\, w_{63} + 28 \\, w_{68} - w_{71} - w_{73} + 2 \\, w_{74} + 28 & w_{62} \\\\\n", "w_{62} & -w_{53} - 2 \\, w_{63} - 2 \\, w_{72} + 12 \\, w_{74}\n", "\\end{array}\\right), \\left(\\begin{array}{rr}\n", "-4 \\, w_{62} + 2 \\, w_{63} + 54 \\, w_{68} - w_{71} - 2 \\, w_{72} & w_{63} \\\\\n", "w_{63} & 2 \\, w_{68} - w_{73}\n", "\\end{array}\\right), \\left(\\begin{array}{rr}\n", "w_{71} & w_{72} \\\\\n", "w_{72} & w_{73}\n", "\\end{array}\\right)\\right]$ " ], "text/plain": [ "[\n", "[-w53 + 2*w62 - 2*w63 - 28*w68 + w71 + w73 - 2*w74 + 28 w62 - w63 - 12*w68 + w72] [w53 - 2*w62 + 2*w63 + 28*w68 - w71 - w73 + 2*w74 + 28 w62] [-4*w62 + 2*w63 + 54*w68 - w71 - 2*w72 w63]\n", "[ w62 - w63 - 12*w68 + w72 w53], [ w62 -w53 - 2*w63 - 2*w72 + 12*w74], [ w63 2*w68 - w73],\n", "\n", "[w71 w72]\n", "[w72 w73]\n", "]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Thanks to the symmetry, Q^* is block diagonal\n", "Q_star = sum(M.T *Q_L*M for M in M_rho) / len(M_rho)\n", "print(\"Blocks on the diagonal of Q:\")\n", "[Q_star[i:i+2, i:i+2] for i in range(0, 8, 2)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The matrix $Q^*$ has therefore the form\n", "$Q^* = diag(Q_1, Q_2, Q_3, Q_4)$,\n", "where the $Q_i$ are $2 \\times 2$ symmetric matrices. $Q^* \\succeq 0$ if and only if $Q_1,\\ldots,Q_4 \\succeq 0$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Dual \n", "\n", "Let us now derive the dual of the previous problem" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rr|rr|rr|rr}\n", "\\mathit{D0}_{00} & \\mathit{D0}_{01} & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "\\mathit{D0}_{01} & \\mathit{D0}_{11} & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "\\hline\n", " 0 & 0 & \\mathit{D1}_{00} & \\mathit{D1}_{01} & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & \\mathit{D1}_{01} & \\mathit{D1}_{11} & 0 & 0 & 0 & 0 \\\\\n", "\\hline\n", " 0 & 0 & 0 & 0 & \\mathit{D2}_{00} & \\mathit{D2}_{01} & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & \\mathit{D2}_{01} & \\mathit{D2}_{11} & 0 & 0 \\\\\n", "\\hline\n", " 0 & 0 & 0 & 0 & 0 & 0 & \\mathit{D3}_{00} & \\mathit{D3}_{01} \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & \\mathit{D3}_{01} & \\mathit{D3}_{11}\n", "\\end{array}\\right)$ " ], "text/plain": [ "[D0_00 D0_01| 0 0| 0 0| 0 0]\n", "[D0_01 D0_11| 0 0| 0 0| 0 0]\n", "[-----------+-----------+-----------+-----------]\n", "[ 0 0|D1_00 D1_01| 0 0| 0 0]\n", "[ 0 0|D1_01 D1_11| 0 0| 0 0]\n", "[-----------+-----------+-----------+-----------]\n", "[ 0 0| 0 0|D2_00 D2_01| 0 0]\n", "[ 0 0| 0 0|D2_01 D2_11| 0 0]\n", "[-----------+-----------+-----------+-----------]\n", "[ 0 0| 0 0| 0 0|D3_00 D3_01]\n", "[ 0 0| 0 0| 0 0|D3_01 D3_11]" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#DQ is the dual variable\n", "DQi = [symb_matrix(2, 'D'+str(i)+'_%d%d', sym=True) for i in range(4)]\n", "DQ = block_diagonal_matrix(DQi)\n", "DQ" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us compute the lagrangian\n", "\n", "$$L(Q, D) = \\langle Q, D \\rangle + \\text{objective function}$$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(- \\mathit{D0}_{00} + \\mathit{D0}_{11} + \\mathit{D1}_{00} - \\mathit{D1}_{11}\\right) w_{53} + \\left(2 \\mathit{D0}_{00} + 2 \\mathit{D0}_{01} - 2 \\mathit{D1}_{00} + 2 \\mathit{D1}_{01} - 4 \\mathit{D2}_{00}\\right) w_{62} + \\left(-2 \\mathit{D0}_{00} - 2 \\mathit{D0}_{01} + 2 \\mathit{D1}_{00} - 2 \\mathit{D1}_{11} + 2 \\mathit{D2}_{00} + 2 \\mathit{D2}_{01}\\right) w_{63} + \\left(-28 \\mathit{D0}_{00} - 24 \\mathit{D0}_{01} + 28 \\mathit{D1}_{00} + 54 \\mathit{D2}_{00} + 2 \\mathit{D2}_{11}\\right) w_{68} + \\left(\\mathit{D0}_{00} - \\mathit{D1}_{00} - \\mathit{D2}_{00} + \\mathit{D3}_{00}\\right) w_{71} + \\left(2 \\mathit{D0}_{01} - 2 \\mathit{D1}_{11} - 2 \\mathit{D2}_{00} + 2 \\mathit{D3}_{01}\\right) w_{72} + \\left(\\mathit{D0}_{00} - \\mathit{D1}_{00} - \\mathit{D2}_{11} + \\mathit{D3}_{11}\\right) w_{73} + \\left(-2 \\mathit{D0}_{00} + 2 \\mathit{D1}_{00} + 12 \\mathit{D1}_{11} + \\frac{1}{70}\\right) w_{74} + 28 \\mathit{D0}_{00} + 28 \\mathit{D1}_{00}$ " ], "text/plain": [ "(-D0_00 + D0_11 + D1_00 - D1_11)*w53 + (2*D0_00 + 2*D0_01 - 2*D1_00 + 2*D1_01 - 4*D2_00)*w62 + (-2*D0_00 - 2*D0_01 + 2*D1_00 - 2*D1_11 + 2*D2_00 + 2*D2_01)*w63 + (-28*D0_00 - 24*D0_01 + 28*D1_00 + 54*D2_00 + 2*D2_11)*w68 + (D0_00 - D1_00 - D2_00 + D3_00)*w71 + (2*D0_01 - 2*D1_11 - 2*D2_00 + 2*D3_01)*w72 + (D0_00 - D1_00 - D2_11 + D3_11)*w73 + (-2*D0_00 + 2*D1_00 + 12*D1_11 + 1/70)*w74 + 28*D0_00 + 28*D1_00" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hessian = lambda pp: jacobian(jacobian(pp, (x,y)), (x,y))\n", "dot = lambda u,v: vector(u) * vector(v)\n", "mdot = lambda A, B: dot(A.list(), B.list())\n", "\n", "\n", "lagrangian = mdot(DQ, Q_star) + objective\n", "var_primal = Q_star.variables()\n", "var_dual = DQ.variables()\n", "lagrangian = QQ[var_dual][var_primal] (lagrangian)\n", "lagrangian" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrrrrrr}\n", "z_{28} & z_{29} & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "z_{29} & -\\frac{7}{6} \\, z_{27} + \\frac{7}{6} \\, z_{28} - \\frac{1}{840} & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & z_{27} & z_{27} - z_{28} - z_{29} + 2 \\, z_{30} & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & z_{27} - z_{28} - z_{29} + 2 \\, z_{30} & -\\frac{1}{6} \\, z_{27} + \\frac{1}{6} \\, z_{28} - \\frac{1}{840} & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & z_{30} & -\\frac{7}{6} \\, z_{27} + \\frac{7}{6} \\, z_{28} + z_{29} - z_{30} - \\frac{1}{840} & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & -\\frac{7}{6} \\, z_{27} + \\frac{7}{6} \\, z_{28} + z_{29} - z_{30} - \\frac{1}{840} & -14 \\, z_{27} + 14 \\, z_{28} + 12 \\, z_{29} - 27 \\, z_{30} & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & z_{27} - z_{28} + z_{30} & -\\frac{1}{6} \\, z_{27} + \\frac{1}{6} \\, z_{28} - z_{29} + z_{30} - \\frac{1}{840} \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & -\\frac{1}{6} \\, z_{27} + \\frac{1}{6} \\, z_{28} - z_{29} + z_{30} - \\frac{1}{840} & -13 \\, z_{27} + 13 \\, z_{28} + 12 \\, z_{29} - 27 \\, z_{30}\n", "\\end{array}\\right)$ " ], "text/plain": [ "[ z28 z29 0 0 0 0 0 0]\n", "[ z29 -7/6*z27 + 7/6*z28 - 1/840 0 0 0 0 0 0]\n", "[ 0 0 z27 z27 - z28 - z29 + 2*z30 0 0 0 0]\n", "[ 0 0 z27 - z28 - z29 + 2*z30 -1/6*z27 + 1/6*z28 - 1/840 0 0 0 0]\n", "[ 0 0 0 0 z30 -7/6*z27 + 7/6*z28 + z29 - z30 - 1/840 0 0]\n", "[ 0 0 0 0 -7/6*z27 + 7/6*z28 + z29 - z30 - 1/840 -14*z27 + 14*z28 + 12*z29 - 27*z30 0 0]\n", "[ 0 0 0 0 0 0 z27 - z28 + z30 -1/6*z27 + 1/6*z28 - z29 + z30 - 1/840]\n", "[ 0 0 0 0 0 0 -1/6*z27 + 1/6*z28 - z29 + z30 - 1/840 -13*z27 + 13*z28 + 12*z29 - 27*z30]" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Find a full rank description of the dual varialbe DQ\n", "lagrang_coeffs = jacobian(lagrangian, var_primal)[0]\n", "lagrang_coeffs = list(lagrang_coeffs)\n", "DQ_star = eliminate_var_matrix(DQ, lagrang_coeffs, 'z')\n", "DQ_star" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# use KKT conditions to transform the SDP to a system of polynomial equations" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "R = QQ[DQ_star.variables() + Q_star.variables()]\n", "\n", "# KKT equations\n", "KKT_eqn = (DQ_star * Q_star).list()\n", "KKT_eqn = list(set(KKT_eqn))\n", "\n", "# ideal generated by the KKT equations\n", "I = map(lambda p: R(p), KKT_eqn)\n", "I = R*I" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ideal I has 17 equations in 12 variables\n" ] } ], "source": [ "print(\"Ideal I is generated by %d equations in %d variables\" % (len(I.gens()), len(R.gens())))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(9 w_{74}^{14} - 2070 w_{74}^{13} + 170888 w_{74}^{12} - 5575088 w_{74}^{11} + 24286192 w_{74}^{10} + 1799298144 w_{74}^{9} - 13606662400 w_{74}^{8} - 140207761920 w_{74}^{7} + 434444465920 w_{74}^{6} + 3454598698496 w_{74}^{5} - 69628131328 w_{74}^{4} - 12362877284352 w_{74}^{3} - 5808164827136 w_{74}^{2} + 2696323768320 w_{74}\\right)\\Bold{Q}[z_{27}, z_{28}, z_{29}, z_{30}, w_{53}, w_{62}, w_{63}, w_{68}, w_{71}, w_{72}, w_{73}, w_{74}]$ " ], "text/plain": [ "Ideal (9*w74^14 - 2070*w74^13 + 170888*w74^12 - 5575088*w74^11 + 24286192*w74^10 + 1799298144*w74^9 - 13606662400*w74^8 - 140207761920*w74^7 + 434444465920*w74^6 + 3454598698496*w74^5 - 69628131328*w74^4 - 12362877284352*w74^3 - 5808164827136*w74^2 + 2696323768320*w74) of Multivariate Polynomial Ring in z27, z28, z29, z30, w53, w62, w63, w68, w71, w72, w73, w74 over Rational Field" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# eliminate all variables except those present in the objective function\n", "I_a = I.elimination_ideal([v for v in R.gens() if SR(v) not in objective.variables()])\n", "I_a" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The ideal I_a has 1 generator\n" ] } ], "source": [ "print(\"The ideal I_a has %d generator\" % len(I_a.gens()))" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}6884147200 \\, {\\left(42875 \\, t^{3} - 40425 \\, t^{2} - 2975 \\, t + 13\\right)} {\\left(3675 \\, t^{2} - 2870 \\, t - 37\\right)} {\\left(35 \\, t + 1\\right)} {\\left(35 \\, t - 1\\right)} {\\left(35 \\, t - 3\\right)} {\\left(35 \\, t - 19\\right)} {\\left(15 \\, t + 1\\right)} {\\left(5 \\, t + 1\\right)} {\\left(5 \\, t - 1\\right)} {\\left(t - 1\\right)} t$ " ], "text/plain": [ "6884147200*(42875*t^3 - 40425*t^2 - 2975*t + 13)*(3675*t^2 - 2870*t - 37)*(35*t + 1)*(35*t - 1)*(35*t - 3)*(35*t - 19)*(15*t + 1)*(5*t + 1)*(5*t - 1)*(t - 1)*t" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "var('t')\n", "p = SR(I_a.gens()[0]).subs({SR(objective.variables()[0]): 70*t})\n", "factor(p)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}t^{3} - \\frac{33}{35} \\, t^{2} - \\frac{17}{245} \\, t + \\frac{13}{42875}$ " ], "text/plain": [ "t^3 - 33/35*t^2 - 17/245*t + 13/42875" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "minimal_poly = p.factor_list()[0][0]\n", "minimal_poly /= minimal_poly.coefficient(t^3)\n", "minimal_poly" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(-\\frac{1}{2} \\, {\\left(i \\, \\sqrt{3} + 1\\right)} {\\left(\\frac{256}{55125} i \\, \\sqrt{3} + \\frac{256}{6125}\\right)}^{\\frac{1}{3}} - \\frac{32 \\, {\\left(-i \\, \\sqrt{3} + 1\\right)}}{525 \\, {\\left(\\frac{256}{55125} i \\, \\sqrt{3} + \\frac{256}{6125}\\right)}^{\\frac{1}{3}}} + \\frac{11}{35}, 1\\right), \\left(-\\frac{1}{2} \\, {\\left(\\frac{256}{55125} i \\, \\sqrt{3} + \\frac{256}{6125}\\right)}^{\\frac{1}{3}} {\\left(-i \\, \\sqrt{3} + 1\\right)} - \\frac{32 \\, {\\left(i \\, \\sqrt{3} + 1\\right)}}{525 \\, {\\left(\\frac{256}{55125} i \\, \\sqrt{3} + \\frac{256}{6125}\\right)}^{\\frac{1}{3}}} + \\frac{11}{35}, 1\\right), \\left({\\left(\\frac{256}{55125} i \\, \\sqrt{3} + \\frac{256}{6125}\\right)}^{\\frac{1}{3}} + \\frac{64}{525 \\, {\\left(\\frac{256}{55125} i \\, \\sqrt{3} + \\frac{256}{6125}\\right)}^{\\frac{1}{3}}} + \\frac{11}{35}, 1\\right)\\right]$ " ], "text/plain": [ "[(-1/2*(I*sqrt(3) + 1)*(256/55125*I*sqrt(3) + 256/6125)^(1/3) - 32/525*(-I*sqrt(3) + 1)/(256/55125*I*sqrt(3) + 256/6125)^(1/3) + 11/35,\n", " 1),\n", " (-1/2*(256/55125*I*sqrt(3) + 256/6125)^(1/3)*(-I*sqrt(3) + 1) - 32/525*(I*sqrt(3) + 1)/(256/55125*I*sqrt(3) + 256/6125)^(1/3) + 11/35,\n", " 1),\n", " ((256/55125*I*sqrt(3) + 256/6125)^(1/3) + 64/525/(256/55125*I*sqrt(3) + 256/6125)^(1/3) + 11/35,\n", " 1)]" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "minimal_poly.roots()" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(-0.07246205065830763, 1\\right), \\left(0.0041380873389235545, 1\\right), \\left(1.0111811061765266, 1\\right)\\right]$ " ], "text/plain": [ "[(-0.07246205065830763, 1),\n", " (0.0041380873389235545, 1),\n", " (1.0111811061765266, 1)]" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "minimal_poly.roots(ring=RDF)" ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "SageMath 8.6", "language": "", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.16" } }, "nbformat": 4, "nbformat_minor": 2 }