In [1]:
using Pkg
Pkg.activate(".")
# Pkg.instantiate()

  Activating project at `~/Mathematics/Teaching/2022/Computational Group Theory/2022 WS Karlsruhe/Symbolic SOS`

Optimization¶

In [2]:
using LinearAlgebra
using DynamicPolynomials
using SumOfSquares
using SCS
In [3]:
@polyvar x y
motzkin = x^4*y^2 + y^4*x^2 - 3*x^2*y^2 + 1
Out[3]:
$$ x^{4}y^{2} + x^{2}y^{4} - 3x^{2}y^{2} + 1 $$
In [4]:
m = SOSModel()
sos_constraint = @constraint m motzkin in SOSCone()
print(m)
$$ \begin{aligned} \text{feasibility}\\ \text{Subject to} \quad & (1)x^{4}y^{2} + (1)x^{2}y^{4} + (-3)x^{2}y^{2} + (1) \text{ is SOS}\\ \end{aligned} $$
In [5]:
set_optimizer(m, SCS.Optimizer)
optimize!(m)
@info termination_status(m)

------------------------------------------------------------------
	       SCS v3.2.1 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 10, constraints m: 20
cones: 	  z: primal zero / dual free vars: 10
	  s: psd vars: 10, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 100000, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct-amd-qdldl
	  nnz(A): 20, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 5.99e+00  2.00e-01  2.40e+00  1.20e+00  1.00e-01  3.76e-04 
    25| 2.80e+13  5.69e+11  3.53e+19  1.76e+19  1.00e-01  5.95e-04 
------------------------------------------------------------------
status:  infeasible
timings: total: 5.98e-04s = setup: 1.09e-04s + solve: 4.89e-04s
	 lin-sys: 4.04e-05s, cones: 2.73e-04s, accel: 3.15e-05s
------------------------------------------------------------------
objective = inf
------------------------------------------------------------------

[ Info: INFEASIBLE
In [6]:
m = Model()
@constraint m sos_constraint (x^2 + y^2 + 1)*motzkin in SOSCone()
print(m)
$$ \begin{aligned} \text{feasibility}\\ \text{Subject to} \quad & (1)x^{6}y^{2} + (2)x^{4}y^{4} + (1)x^{2}y^{6} + (-2)x^{4}y^{2} + (-2)x^{2}y^{4} + (-3)x^{2}y^{2} + (1)x^{2} + (1)y^{2} + (1) \text{ is SOS}\\ \end{aligned} $$
In [7]:
set_optimizer(m, SCS.Optimizer)
optimize!(m)
@info termination_status(m)

------------------------------------------------------------------
	       SCS v3.2.1 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 45, constraints m: 72
cones: 	  z: primal zero / dual free vars: 27
	  s: psd vars: 45, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 100000, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct-amd-qdldl
	  nnz(A): 90, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 3.00e+00  2.61e-01  1.87e+00  9.35e-01  1.00e-01  3.17e-04 
    50| 8.10e-07  2.31e-07  1.58e-06  7.88e-07  1.00e-01  1.39e-03 
------------------------------------------------------------------
status:  solved
timings: total: 1.39e-03s = setup: 1.53e-04s + solve: 1.24e-03s
	 lin-sys: 1.12e-04s, cones: 8.32e-04s, accel: 5.73e-05s
------------------------------------------------------------------
objective = 0.000001
------------------------------------------------------------------

[ Info: OPTIMAL
In [8]:
sos_dec = sos_decomposition(m[:sos_constraint], 1e-6);
res = polynomial(sos_dec) - (x^2 + y^2 + 1)*motzkin
Out[8]:
$$ -8.682134544635289e-8x^{6}y^{2} + 3.7611829760295325e-16x^{5}y^{3} - 2.389406228431312e-7x^{4}y^{4} - 1.6640317781952892e-16x^{3}y^{5} - 8.682134611248671e-8x^{2}y^{6} - 1.8919396178657707e-16x^{5}y^{2} - 1.9499281370497737e-16x^{4}y^{3} + 6.052739527431082e-16x^{3}y^{4} - 4.2021135102739127e-16x^{2}y^{5} - 4.390374863660895e-8x^{4}y^{2} - 4.440892098500631e-16x^{3}y^{3} - 4.390374730434132e-8x^{2}y^{4} - 4.258723168891021e-16x^{4}y - 3.768555672296405e-16x^{3}y^{2} - 1.3854979075951675e-16x^{2}y^{3} - 1.7327769156731698e-16xy^{4} + 5.523413228733438e-16x^{3}y + 3.040020368771934e-7x^{2}y^{2} + 6.68904004214329e-16xy^{3} + 2.5688464354345762e-16x^{2}y - 5.637859011547065e-16xy^{2} + 1.311590251962258e-7x^{2} - 9.992007221626391e-16xy + 1.3115902319782435e-7y^{2} + 6.978400064696924e-16x + 9.227568172392428e-16y - 6.583188838504839e-8 $$
In [9]:
norm(res, 1)
Out[9]:
1.1325427921717647e-6
In [10]:
@info "" rounded_solution=round(polynomial(sos_dec), digits=6) original_motzkin=(x^2 + y^2 + 1.0)*motzkin

┌ Info: 
│   rounded_solution = x⁶y² + 2.0x⁴y⁴ + x²y⁶ - 2.0x⁴y² - 2.0x²y⁴ - 3.0x²y² + x² + y² + 1.0
└   original_motzkin = x⁶y² + 2.0x⁴y⁴ + x²y⁶ - 2.0x⁴y² - 2.0x²y⁴ - 3.0x²y² + x² + y² + 1.0
In [11]:
G = gram_matrix(m[:sos_constraint]);
# @show norm(G.Q - round.(G.Q, digits = 2), 1)
round.(G.Q, digits=10)
Out[11]:
9×9 Matrix{Float64}:
  1.0        0.0       0.499997  -0.0       …   0.0        0.0       -0.0
  0.0        1.00001  -0.0       -0.0          -0.0        0.0       -1.0
  0.499997  -0.0       1.0       -0.0           0.0       -0.0       -0.0
 -0.0       -0.0      -0.0        0.999988      0.0       -0.999994   0.0
  0.0        0.0       0.0        0.0          -0.999994   0.0       -0.0
 -1.49999   -0.0      -1.49999    0.0       …  -0.0        0.0        0.0
  0.0       -0.0       0.0        0.0           1.0       -0.0       -0.0
  0.0        0.0      -0.0       -0.999994     -0.0        1.0        0.0
 -0.0       -1.0      -0.0        0.0          -0.0        0.0        1.0
In [12]:
g = G.basis.monomials' * round.(Int, 2G.Q) * G.basis.monomials
g
Out[12]:
$$ 2x^{6}y^{2} + 4x^{4}y^{4} + 2x^{2}y^{6} - 4x^{4}y^{2} - 4x^{2}y^{4} - 6x^{2}y^{2} + 2x^{2} + 2y^{2} + 2 $$
In [13]:
g - 2(x^2 + y^2 + 1)*motzkin
Out[13]:
$$ 0 $$
In [14]:
round.(Int, 2G.Q)
Out[14]:
9×9 Matrix{Int64}:
  2   0   1   0   0  -3   0   0   0
  0   2   0   0   0   0   0   0  -2
  1   0   2   0   0  -3   0   0   0
  0   0   0   2   0   0   0  -2   0
  0   0   0   0   2   0  -2   0   0
 -3   0  -3   0   0   6   0   0   0
  0   0   0   0  -2   0   2   0   0
  0   0   0  -2   0   0   0   2   0
  0  -2   0   0   0   0   0   0   2

Finding the denominator¶

In [15]:
m = SOSModel()
# X = DynamicPolynomials.monomials([x,y], 2:2)
X = [x^2, y^2, 1]
Out[15]:
3-element Vector{Term{true, Int64}}:
 x²
 y²
 1
In [16]:
@variable m denom Poly(X)
Out[16]:
$$ (_[1])x^{2} + (_[2])y^{2} + (_[3]) $$
In [17]:
@constraint m denom*motzkin ∈ SOSCone()
Out[17]:
$$ (_[1])x^{6}y^{2} + (_[1] + _[2])x^{4}y^{4} + (_[2])x^{2}y^{6} + (-3 _[1] + _[3])x^{4}y^{2} + (-3 _[2] + _[3])x^{2}y^{4} + (-3 _[3])x^{2}y^{2} + (_[1])x^{2} + (_[2])y^{2} + (_[3]) \text{ is SOS} $$
In [18]:
@constraint m sum(coefficients(denom)) >= 1
Out[18]:
$$ {\_}_{1} + {\_}_{2} + {\_}_{3} \geq 1.0 $$
In [19]:
print(m)
$$ \begin{aligned} \text{feasibility}\\ \text{Subject to} \quad & {\_}_{1} + {\_}_{2} + {\_}_{3} \geq 1.0\\ & (_[1])x^{6}y^{2} + (_[1] + _[2])x^{4}y^{4} + (_[2])x^{2}y^{6} + (-3 _[1] + _[3])x^{4}y^{2} + (-3 _[2] + _[3])x^{2}y^{4} + (-3 _[3])x^{2}y^{2} + (_[1])x^{2} + (_[2])y^{2} + (_[3]) \text{ is SOS}\\ \end{aligned} $$
In [20]:
set_optimizer(m, SCS.Optimizer)
optimize!(m)
@info termination_status(m)

------------------------------------------------------------------
	       SCS v3.2.1 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 48, constraints m: 73
cones: 	  z: primal zero / dual free vars: 27
	  l: linear vars: 1
	  s: psd vars: 45, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 100000, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct-amd-qdldl
	  nnz(A): 105, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 3.41e-01  5.94e-02  1.14e-01  5.71e-02  1.00e-01  2.36e-04 
    25| 2.18e-06  2.97e-07  0.00e+00  0.00e+00  1.00e-01  8.06e-04 
------------------------------------------------------------------
status:  solved
timings: total: 8.09e-04s = setup: 1.32e-04s + solve: 6.77e-04s
	 lin-sys: 6.06e-05s, cones: 4.37e-04s, accel: 3.27e-05s
------------------------------------------------------------------
objective = 0.000000
------------------------------------------------------------------

[ Info: OPTIMAL
In [21]:
denominator = round.(value(m[:denom]), digits=5)
Out[21]:
$$ 0.25621x^{2} + 0.25621y^{2} + 0.54499 $$
In [ ]: