{ "cells": [ { "cell_type": "markdown", "source": [ "# On the importance of Dualization" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "using DynamicPolynomials\n", "using SumOfSquares" ], "metadata": {}, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "Sum-of-Squares programs are usually solved by SemiDefinite Programming solvers (SDPs).\n", "These programs can be represented into two different formats:\n", "Either the *standard conic form*, also known as *kernel form*:\n", "$$\n", "\\begin{aligned}\n", " \\min\\limits_{Q \\in \\mathbb{S}_n} & \\langle C, Q \\rangle\\\\\n", " \\text{subject to:} & \\langle A_i, Q \\rangle = b_i, \\quad i=1,2,\\ldots,m\\\\\n", " & Q \\succeq 0,\n", "\\end{aligned}\n", "$$\n", "or the *geometric conic form*, also known as *image form*:\n", "$$\n", "\\begin{aligned}\n", " \\max\\limits_{y \\in \\mathbb{R}^m} & \\langle b, y \\rangle\\\\\n", " \\text{subject to:} & C \\succeq \\sum_{i=1}^m A_i y_i\\\\\n", " & y\\ \\mathsf{free},\n", "\\end{aligned}\n", "$$" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "In this tutorial, we investigate in which of these two forms a Sum-of-Squares\n", "constraint should be written into.\n", "Consider the simple example of trying to determine whether the following univariate\n", "polynomial is a Sum-of-Squares:" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "------------------------------------------------------------------\n", "\t SCS v3.2.4 - Splitting Conic Solver\n", "\t(c) Brendan O'Donoghue, Stanford University, 2012\n", "------------------------------------------------------------------\n", "problem: variables n: 6, constraints m: 11\n", "cones: \t z: primal zero / dual free vars: 5\n", "\t s: psd vars: 6, ssize: 1\n", "settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07\n", "\t alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1\n", "\t max_iters: 100000, normalize: 1, rho_x: 1.00e-06\n", "\t acceleration_lookback: 10, acceleration_interval: 10\n", "\t compiled with openmp parallelization enabled\n", "lin-sys: sparse-direct-amd-qdldl\n", "\t nnz(A): 12, nnz(P): 0\n", "------------------------------------------------------------------\n", " iter | pri res | dua res | gap | obj | scale | time (s)\n", "------------------------------------------------------------------\n", " 0| 1.30e+01 1.87e+00 2.95e+01 1.47e+01 1.00e-01 2.57e-04 \n", " 50| 6.01e-07 2.84e-07 1.95e-06 -9.76e-07 1.00e-01 4.16e-04 \n", "------------------------------------------------------------------\n", "status: solved\n", "timings: total: 4.18e-04s = setup: 1.06e-04s + solve: 3.11e-04s\n", "\t lin-sys: 1.96e-05s, cones: 1.81e-04s, accel: 3.81e-06s\n", "------------------------------------------------------------------\n", "objective = -0.000001\n", "------------------------------------------------------------------\n" ] } ], "cell_type": "code", "source": [ "import SCS\n", "@polyvar x\n", "p = (x + 1)^2 * (x + 2)^2\n", "model_scs = Model(SCS.Optimizer)\n", "con_ref = @constraint(model_scs, p in SOSCone())\n", "optimize!(model_scs)" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "As we can see in the log, SCS reports `6` variables and `11` constraints.\n", "We can also choose to dualize the problem before it is\n", "passed to SCS as follows:" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "------------------------------------------------------------------\n", "\t SCS v3.2.4 - Splitting Conic Solver\n", "\t(c) Brendan O'Donoghue, Stanford University, 2012\n", "------------------------------------------------------------------\n", "problem: variables n: 5, constraints m: 6\n", "cones: \t s: psd vars: 6, ssize: 1\n", "settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07\n", "\t alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1\n", "\t max_iters: 100000, normalize: 1, rho_x: 1.00e-06\n", "\t acceleration_lookback: 10, acceleration_interval: 10\n", "\t compiled with openmp parallelization enabled\n", "lin-sys: sparse-direct-amd-qdldl\n", "\t nnz(A): 6, nnz(P): 0\n", "------------------------------------------------------------------\n", " iter | pri res | dua res | gap | obj | scale | time (s)\n", "------------------------------------------------------------------\n", " 0| 1.97e+02 1.05e+01 4.39e+03 -2.20e+03 1.00e-01 1.23e-04 \n", " 50| 6.29e-06 1.81e-06 2.22e-05 -1.11e-05 1.00e-01 2.71e-04 \n", "------------------------------------------------------------------\n", "status: solved\n", "timings: total: 2.72e-04s = setup: 4.09e-05s + solve: 2.31e-04s\n", "\t lin-sys: 1.37e-05s, cones: 1.29e-04s, accel: 2.91e-06s\n", "------------------------------------------------------------------\n", "objective = -0.000011\n", "------------------------------------------------------------------\n" ] } ], "cell_type": "code", "source": [ "using Dualization\n", "model_dual_scs = Model(dual_optimizer(SCS.Optimizer))\n", "@objective(model_dual_scs, Max, 0.0)\n", "con_ref = @constraint(model_dual_scs, p in SOSCone())\n", "optimize!(model_dual_scs)" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "This time, SCS reports `5` variables and `6` constraints." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## Bridges operating behind the scenes\n", "\n", "The difference comes from the fact that, when designing the JuMP interface of\n", "SCS, it was decided that the model would be read in the image form.\n", "SCS therefore declares that it only supports free variables, represented in\n", "JuMP as variables in `MOI.Reals` and affine semidefinite constraints,\n", "represented in JuMP as\n", "`MOI.VectorAffineFunction`-in-`MOI.PositiveSemidefiniteConeTriangle`\n", "constraints.\n", "On the other hand, SumOfSquares gave the model in kernel form so the\n", "positive semidefinite (PSD) variables were reformulated as free variables\n", "constrained to be PSD using an affine PSD constraints.\n", "\n", "This transformation is done transparently without warning but it can be\n", "inspected using `print_active_bridges`.\n", "As shown below, we can see\n", "`Unsupported variable: MOI.PositiveSemidefiniteConeTriangle` and\n", "`adding as constraint`\n", "indicating that PSD variables are not supported and they are added as free\n", "variables.\n", "Then we have `Unsupported constraint: MOI.VectorOfVariables-in-MOI.PositiveSemidefiniteConeTriangle`\n", "indicating that SCS does not support constraining variables in the PSD cone\n", "so it will just convert it into affine expressions in the PSD cone.\n", "Of course, this is equivalent but it means that SCS will not exploit this\n", "particular structure of the problem hence solving might be less efficient." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " * Supported objective: MOI.ScalarAffineFunction{Float64}\n", " * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-SumOfSquares.SOSPolynomialSet{FullSpace, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}\n", " | bridged by:\n", " | SumOfSquares.Bridges.Constraint.SOSPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, FullSpace, Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}, Vector{Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}}}}, Union{SumOfSquares.EmptyCone, SumOfSquares.PositiveSemidefinite2x2ConeTriangle, MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle}, MOI.PositiveSemidefiniteConeTriangle, MonomialBasis, MonomialBasis, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}\n", " | may introduce:\n", " | * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-PolyJuMP.ZeroPolynomialSet{FullSpace, MonomialBasis, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}\n", " | | bridged by:\n", " | | PolyJuMP.Bridges.Constraint.ZeroPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}\n", " | | may introduce:\n", " | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Zeros\n", " | * Unsupported variable: SumOfSquares.EmptyCone\n", " | | adding as constraint:\n", " | | * Supported variable: MOI.Reals\n", " | | * Unsupported constraint: MOI.VectorOfVariables-in-SumOfSquares.EmptyCone\n", " | | | bridged by:\n", " | | | SumOfSquares.Bridges.Constraint.EmptyBridge{Float64, MOI.VectorOfVariables}\n", " | | | may introduce:\n", " | * Unsupported variable: MOI.Nonnegatives\n", " | | adding as constraint:\n", " | | * Supported variable: MOI.Reals\n", " | | * Unsupported constraint: MOI.VectorOfVariables-in-MOI.Nonnegatives\n", " | | | bridged by:\n", " | | | MOIB.Constraint.FunctionConversionBridge{Float64, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables, MOI.Nonnegatives}\n", " | | | may introduce:\n", " | | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Nonnegatives\n", " | * Unsupported variable: SumOfSquares.PositiveSemidefinite2x2ConeTriangle\n", " | | bridged by:\n", " | | SumOfSquares.Bridges.Variable.PositiveSemidefinite2x2Bridge{Float64}\n", " | | may introduce:\n", " | | * Unsupported variable: MOI.RotatedSecondOrderCone\n", " | | | adding as constraint:\n", " | | | * Supported variable: MOI.Reals\n", " | | | * Unsupported constraint: MOI.VectorOfVariables-in-MOI.RotatedSecondOrderCone\n", " | | | | bridged by:\n", " | | | | MOIB.Constraint.RSOCtoSOCBridge{Float64, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables}\n", " | | | | may introduce:\n", " | | | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.SecondOrderCone\n", " | * Unsupported variable: MOI.PositiveSemidefiniteConeTriangle\n", " | | adding as constraint:\n", " | | * Supported variable: MOI.Reals\n", " | | * Unsupported constraint: MOI.VectorOfVariables-in-MOI.PositiveSemidefiniteConeTriangle\n", " | | | bridged by:\n", " | | | MOIB.Constraint.SetDotScalingBridge{Float64, MOI.PositiveSemidefiniteConeTriangle, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables}\n", " | | | may introduce:\n", " | | | * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}\n", " | | | | bridged by:\n", " | | | | SCS.ScaledPSDConeBridge{Float64, MOI.VectorAffineFunction{Float64}}\n", " | | | | may introduce:\n", " | | | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-SCS.ScaledPSDCone\n" ] } ], "cell_type": "code", "source": [ "print_active_bridges(model_scs)" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "With the dual version, we can see that variables in the PSD cone are supported\n", "directly hence we don't need that extra conversion." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " * Supported objective: MOI.ScalarAffineFunction{Float64}\n", " * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-SumOfSquares.SOSPolynomialSet{FullSpace, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}\n", " | bridged by:\n", " | SumOfSquares.Bridges.Constraint.SOSPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, FullSpace, Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}, Vector{Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}}}}, Union{SumOfSquares.EmptyCone, SumOfSquares.PositiveSemidefinite2x2ConeTriangle, MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle}, MOI.PositiveSemidefiniteConeTriangle, MonomialBasis, MonomialBasis, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}\n", " | may introduce:\n", " | * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-PolyJuMP.ZeroPolynomialSet{FullSpace, MonomialBasis, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}\n", " | | bridged by:\n", " | | PolyJuMP.Bridges.Constraint.ZeroPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}\n", " | | may introduce:\n", " | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Zeros\n", " | * Unsupported variable: SumOfSquares.EmptyCone\n", " | | adding as constraint:\n", " | | * Supported variable: MOI.Reals\n", " | | * Unsupported constraint: MOI.VectorOfVariables-in-SumOfSquares.EmptyCone\n", " | | | bridged by:\n", " | | | SumOfSquares.Bridges.Constraint.EmptyBridge{Float64, MOI.VectorOfVariables}\n", " | | | may introduce:\n", " | * Supported variable: MOI.Nonnegatives\n", " | * Unsupported variable: SumOfSquares.PositiveSemidefinite2x2ConeTriangle\n", " | | bridged by:\n", " | | SumOfSquares.Bridges.Variable.PositiveSemidefinite2x2Bridge{Float64}\n", " | | may introduce:\n", " | | * Supported variable: MOI.RotatedSecondOrderCone\n", " | * Supported variable: MOI.PositiveSemidefiniteConeTriangle\n" ] } ], "cell_type": "code", "source": [ "print_active_bridges(model_dual_scs)" ], "metadata": {}, "execution_count": 5 }, { "cell_type": "markdown", "source": [ "## In more details\n", "\n", "Consider a polynomial\n", "$$\n", "p(x) = \\sum_{\\alpha} p_\\alpha x^\\alpha,\n", "$$\n", "a vector of monomials `b(x)` and the set\n", "$$\n", "\\mathcal{A}_\\alpha = \\{\\,(\\beta, \\gamma) \\in b(x)^2 \\mid x^\\beta x^\\gamma = x^\\alpha\\,\\}\n", "$$\n", "The constraint encoding the existence of a PSD matrix `Q` such that `p(x) = b(x)' * Q * b(x)`\n", "can be written in standard conic form as follows:\n", "$$\n", "\\begin{aligned}\n", " \\langle \\sum_{(\\beta, \\gamma) \\in \\mathcal{A}_\\alpha} e_\\beta e_\\gamma^\\top, Q \\rangle & = p_\\alpha, \\quad\\forall \\alpha\\\\\n", " Q & \\succeq 0\n", "\\end{aligned}\n", "$$\n", "Given an arbitrary choice of elements in each set $\\mathcal{A}_\\alpha$:\n", "$(\\beta_\\alpha, \\gamma_\\alpha) \\in \\mathcal{A}_\\alpha$.\n", "It can also equivalently be written in the geometric conic form as follows:\n", "$$\n", "\\begin{aligned}\n", " p_\\alpha e_{\\beta_\\alpha} e_{\\gamma_\\alpha}^\\top +\n", " \\sum_{(\\beta, \\gamma) \\in \\mathcal{A}_\\alpha \\setminus (\\beta_\\alpha, \\gamma_\\alpha)}\n", " y_{\\beta,\\gamma} (e_\\beta e_\\gamma - e_{\\beta_\\alpha} e_{\\gamma_\\alpha}^\\top)^\\top\n", " & \\succeq 0\\\\\n", " y_{\\beta,\\gamma} & \\text{ free}\n", "\\end{aligned}\n", "$$\n", "\n", "## Should I dualize or not ?\n", "\n", "Let's study the evolution of the dimensions `m` and `n` of the semidefinite\n", "program in two extreme examples and then try to extrapolate from these.\n", "\n", "### Univariate case\n", "\n", "Suppose `p` is a univariate polynomial of degree $2d$.\n", "Then `n` will be equal to `d(d + 1)/2` for both the standard and geometric conic forms.\n", "On the other hand, `m` will be equal to `2d + 1` for the standard conic form and\n", "`d(d + 1) / 2 - (2d + 1)` for the geometric form case.\n", "So `m` grows **linearly** for the kernel form but **quadratically** for the image form!\n", "\n", "### Quadratic case\n", "\n", "Suppose `p` is a quadratic form of `d` variables.\n", "Then `n` will be equal to `d` for both the standard and geometric conic forms.\n", "On the other hand, `m` will be equal to `d(d + 1)/2` for the standard conic form and\n", "`0` for the geometric form case.\n", "So `m` grows **quadratically** for the kernel form but is zero for the image form!\n", "\n", "### In general\n", "\n", "In general, if $s_d$ is the dimension of the space of polynomials of degree `d` then\n", "$m = s_{2d}$ for the kernel form and $m = s_{d}(s_{d} + 1)/2 - s_{2d}$ for the image form.\n", "As a rule of thumb, the kernel form will have a smaller `m` if `p` has a low number of variables\n", "and low degree and vice versa.\n", "Of course, you can always try with and without Dualization and see which one works best." ], "metadata": {} }, { "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 }