Notebook

Installation¶

The following instructions were prepared using julia-1.1.1.

Before exploring the notebook you need to clone the main repository:

git clone https://github.com/kalmarek/1812.03456.git

This notebook should be located in 1812.03456/notebooks directory.

In the main directory (1812.03456) you should run the following code in julias REPL console to instantiate the environment for computations:

using Pkg
Pkg.activate(".")
Pkg.instantiate()

(this needs to be done once per installation).

Instantiation should install (among others): the SCS solver, JuMP package for mathematical programming and IntervalArithmetic.jl package from ValidatedNumerics.jl.

The environment uses Groups.jl, GroupRings.jl (which are built on the framework of AbstractAlgebra.jl) and PropertyT.jl packages.

The computation¶

The following programme certifies that

$\operatorname{Adj}_4 + \operatorname{Op}_4 - 0.82\Delta_4 =\Sigma_i \xi_i^*\xi_i \in \Sigma^2_2\mathbb{R}\operatorname{SL}(4,\mathbb{Z}).$

With small changes (which we will indicate) it also certifies that

$\operatorname{Adj}_3 - 0.157999\Delta_3 \in \Sigma^2_2\mathbb{R}\operatorname{SL}(3,\mathbb{Z})$ and that

$\operatorname{Adj}_5 +1.5 \mathrm{Op}_5 - 1.5\Delta_5 \in \Sigma^2_2\mathbb{R}\operatorname{SL}(5,\mathbb{Z}).$

In [1]:

using Pkg
Pkg.activate("..")
using Dates
now()

Out[1]:

2019-07-05T22:42:41.473

In [2]:

using LinearAlgebra
using AbstractAlgebra
using Groups
using GroupRings
using PropertyT

So far we only made the needed packages available in the notebook. In the next cell we define G to be the set of all $4\times 4$ matrices over $\mathbb Z$ . (For the second computation, set N=3 below; for the third, set N=5)

In [3]:

N = 4
G = MatrixAlgebra(zz, N)

Out[3]:

Matrix Algebra of degree 4 over Integers

Generating set¶

Now we create the elementary matrices $E_{i,j}$ . The set of all such matrices and their inverses is denoted by S.

In [4]:

S = PropertyT.generating_set(G)

Out[4]:

24-element Array{AbstractAlgebra.Generic.MatAlgElem{Int64},1}:
 [1 1 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1] 
 [1 0 1 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1] 
 [1 0 0 1]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1] 
 [1 0 0 0]
[1 1 0 0]
[0 0 1 0]
[0 0 0 1] 
 [1 0 0 0]
[0 1 1 0]
[0 0 1 0]
[0 0 0 1] 
 [1 0 0 0]
[0 1 0 1]
[0 0 1 0]
[0 0 0 1] 
 [1 0 0 0]
[0 1 0 0]
[1 0 1 0]
[0 0 0 1] 
 [1 0 0 0]
[0 1 0 0]
[0 1 1 0]
[0 0 0 1] 
 [1 0 0 0]
[0 1 0 0]
[0 0 1 1]
[0 0 0 1] 
 [1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[1 0 0 1] 
 [1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 1 0 1] 
 [1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 1 1] 
 [1 -1 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1]
 [1 0 -1 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1]
 [1 0 0 -1]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1]
 [1 0 0 0]
[-1 1 0 0]
[0 0 1 0]
[0 0 0 1]
 [1 0 0 0]
[0 1 -1 0]
[0 0 1 0]
[0 0 0 1]
 [1 0 0 0]
[0 1 0 -1]
[0 0 1 0]
[0 0 0 1]
 [1 0 0 0]
[0 1 0 0]
[-1 0 1 0]
[0 0 0 1]
 [1 0 0 0]
[0 1 0 0]
[0 -1 1 0]
[0 0 0 1]
 [1 0 0 0]
[0 1 0 0]
[0 0 1 -1]
[0 0 0 1]
 [1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[-1 0 0 1]
 [1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 -1 0 1]
 [1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 -1 1]

Group Ring and Laplacians¶

Now we will generate the ball E_R of radius $R=4$ in $\operatorname{SL}(N,\mathbb{Z})$ and use this as a (partial) basis in a group ring (denoted by RG below). Such group ring also needs a multiplication table (pm, which is actually a division table) which is created as follows: when $x,y$ reside at positions i-th and j-th in E_R, then pm[i,j] = k, where k is the position of $x^{-1}y$ in E_R.

In [5]:

halfradius = 2
E_R, sizes = Groups.generate_balls(S, radius=2*halfradius);
E_rdict = GroupRings.reverse_dict(E_R)
pm = GroupRings.create_pm(E_R, E_rdict, sizes[halfradius]; twisted=true);
RG = GroupRing(G, E_R, E_rdict, pm)
@show sizes;
Δ = length(S)*one(RG) - sum(RG(s) for s in S)

sizes = [25, 433, 6149, 75197]

Out[5]:

24[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1] - 1[1 1 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1] - 1[1 0 1 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1] - 1[1 0 0 1]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1] - 1[1 0 0 0]
[1 1 0 0]
[0 0 1 0]
[0 0 0 1] - 1[1 0 0 0]
[0 1 1 0]
[0 0 1 0]
[0 0 0 1] - 1[1 0 0 0]
[0 1 0 1]
[0 0 1 0]
[0 0 0 1] - 1[1 0 0 0]
[0 1 0 0]
[1 0 1 0]
[0 0 0 1] - 1[1 0 0 0]
[0 1 0 0]
[0 1 1 0]
[0 0 0 1] - 1[1 0 0 0]
[0 1 0 0]
[0 0 1 1]
[0 0 0 1] - 1[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[1 0 0 1] - 1[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 1 0 1] - 1[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 1 1] - 1[1 -1 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1] - 1[1 0 -1 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1] - 1[1 0 0 -1]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1] - 1[1 0 0 0]
[-1 1 0 0]
[0 0 1 0]
[0 0 0 1] - 1[1 0 0 0]
[0 1 -1 0]
[0 0 1 0]
[0 0 0 1] - 1[1 0 0 0]
[0 1 0 -1]
[0 0 1 0]
[0 0 0 1] - 1[1 0 0 0]
[0 1 0 0]
[-1 0 1 0]
[0 0 0 1] - 1[1 0 0 0]
[0 1 0 0]
[0 -1 1 0]
[0 0 0 1] - 1[1 0 0 0]
[0 1 0 0]
[0 0 1 -1]
[0 0 0 1] - 1[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[-1 0 0 1] - 1[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 -1 0 1] - 1[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 -1 1]

Orbit Decomposition¶

Now something happens: in the next cell we split the subspace of $\mathbb{R} \operatorname{SL}(N, \mathbb{Z})$ supported on E_R into irreducible representations of the wreath product $\mathbb Z / 2 \mathbb Z \wr \operatorname{Sym}_N$ . The action of wreath product on the elements of the matrix space is by conjugation, i.e. permutation of rows and columns. We also compute projections on the invariant subspaces to later speed up the optimisation step.

In [6]:

od = PropertyT.OrbitData(RG, WreathProduct(PermGroup(2), PermGroup(N)))
orbit_data = PropertyT.decimate(od);

┌ Info: Decomposing basis of RG into orbits of
│   autS = Wreath Product of Permutation group over 2 elements by Permutation group over 4 elements
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:15

  0.381708 seconds (1.18 M allocations: 82.083 MiB, 6.45% gc time)

┌ Info: The action has 558 orbits
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:18
┌ Info: Finding projections in the Group Ring of
│   autS = Wreath Product of Permutation group over 2 elements by Permutation group over 4 elements
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:20

  6.154452 seconds (9.74 M allocations: 487.622 MiB, 4.69% gc time)

┌ Info: Finding AutS-action matrix representation
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:23

  0.294136 seconds (963.71 k allocations: 64.184 MiB, 6.98% gc time)
  0.440521 seconds (717.96 k allocations: 49.094 MiB, 4.87% gc time)

┌ Info: Computing the projection matrices Uπs
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:27

  1.778654 seconds (1.75 M allocations: 841.323 MiB, 5.12% gc time)

┌ Info: 
│ multiplicities  =   3  13  19  12  10   0   0   0   9  11  13  15   0   0   0   1   1   1   2   1
│     dimensions  =   1   3   3   2   1   4   8   4   6   6   6   6   4   8   4   1   3   3   2   1
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:37
┌ Info: Sparsifying (433, 3)-matrix... 
│  0.7013086989992302   → 0.18475750577367206 ), (240 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 13)-matrix... 
│  0.5821637946349263   → 0.3084029134837449  ), (1736 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 19)-matrix... 
│  0.6031360155585268   → 0.4934970220007293  ), (4060 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 12)-matrix... 
│  0.6123941493456505   → 0.2817551963048499  ), (1464 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 10)-matrix... 
│  0.48221709006928404  → 0.1                 ), (433 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 9)-matrix... 
│  0.24044136515268155  → 0.03695150115473441 ), (144 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 11)-matrix... 
│  0.19462523619567498  → 0.033592273777031285), (160 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 13)-matrix... 
│  0.16237342334340024  → 0.02842423165748801 ), (160 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 15)-matrix... 
│  0.16874518860662047  → 0.027097767513471902), (176 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 1)-matrix... 
│  0.11085450346420324  → 0.11085450346420324 ), (48 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 1)-matrix... 
│  0.07390300230946882  → 0.07390300230946882 ), (32 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 1)-matrix... 
│  0.07390300230946882  → 0.07390300230946882 ), (32 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 2)-matrix... 
│  0.11200923787528869  → 0.09237875288683603 ), (80 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114
┌ Info: Sparsifying (433, 1)-matrix... 
│  0.11085450346420324  → 0.11085450346420324 ), (48 non-zeros)
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/orbitdata.jl:114

Elements Adj and Op¶

Now we define the elements $\operatorname{Adj}_N$ and $\operatorname{Op}_N$ . The functions Sq, Adj, Op returning the appropriate elements are defined in the src/sqadjop.jl source file.

In [7]:

@time AdjN = PropertyT.Adj(RG, N)
@time OpN = PropertyT.Op(RG, N);

  1.719218 seconds (2.70 M allocations: 134.807 MiB, 3.24% gc time)
  0.284995 seconds (367.56 k allocations: 18.374 MiB, 2.55% gc time)

Finally we compute the element elt of our interest:

if N=3: $\operatorname{elt} = \operatorname{Adj}_3$
if N=4: $\operatorname{elt} = \operatorname{Adj}_4 + \operatorname{Op}_4$
if N=5: $\operatorname{elt} = \operatorname{Adj}_5 + 1.5\operatorname{Op}_5.$

In [8]:

if N == 3
    k = 0
elseif N == 4
    k = 1
elseif N == 5
    k = 1.5
end
elt = AdjN + k*OpN;
elt.coeffs

Out[8]:

75197-element SparseArrays.SparseVector{Int64,Int64} with 361 stored entries:
  [1    ]  =  480
  [2    ]  =  -40
  [3    ]  =  -40
  [4    ]  =  -40
  [5    ]  =  -40
  [6    ]  =  -40
  [7    ]  =  -40
  [8    ]  =  -40
  [9    ]  =  -40
  [10   ]  =  -40
           ⋮
  [418  ]  =  1
  [420  ]  =  1
  [422  ]  =  2
  [423  ]  =  1
  [424  ]  =  1
  [425  ]  =  1
  [426  ]  =  1
  [428  ]  =  1
  [429  ]  =  1
  [430  ]  =  1
  [431  ]  =  1

Optimization Problem¶

We are ready to define the optimisation problem. Function

PropertyT.SOS_problem(x, Δ, orbit_data; upper_bound=UB)

defines the optimisation problem equivalent to the one of the form

$\begin{align} \text{ maximize : } \quad & \lambda\\ \text{under constraints : }\quad & 0 \leqslant \lambda \leqslant \operatorname{UB},\\ & x - \lambda \Delta = \sum \xi_i^* \xi_i,\\ & \text{each $\xi_i$ is invariant under $\mathbb{Z}/2\mathbb{Z} \wr \operatorname{Sym}_N$}. \end{align}$

In [9]:

# @time SDP_problem, varλ, varP = PropertyT.SOS_problem(elt, Δ, orbit_data)
if N == 3
    UB = 0.158
elseif N == 4
    UB = 0.82005
elseif N == 5
    UB = 1.5005
end
SDP_problem, varP = PropertyT.SOS_problem(elt, Δ, orbit_data; upper_bound=UB)

┌ Info: Adding 558 constraints...
└ @ PropertyT /home/kalmar/.julia/packages/PropertyT/yHEMH/src/sos_sdps.jl:92

  1.916622 seconds (3.18 M allocations: 408.777 MiB, 12.54% gc time)

Out[9]:

(A JuMP Model
Maximization problem with:
Variables: 1388
Objective function type: JuMP.VariableRef
`JuMP.VariableRef`-in-`MathOptInterface.LessThan{Float64}`: 1 constraint
`JuMP.GenericAffExpr{Float64,JuMP.VariableRef}`-in-`MathOptInterface.EqualTo{Float64}`: 558 constraints
`Array{JuMP.VariableRef,1}`-in-`MathOptInterface.PositiveSemidefiniteConeSquare`: 14 constraints
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
Names registered in the model: λ, Array{JuMP.VariableRef,2}[[noname noname noname; noname noname noname; noname noname noname], [noname noname … noname noname; noname noname … noname noname; … ; noname noname … noname noname; noname noname … noname noname], [noname noname … noname noname; noname noname … noname noname; … ; noname noname … noname noname; noname noname … noname noname], [noname noname … noname noname; noname noname … noname noname; … ; noname noname … noname noname; noname noname … noname noname], [noname noname … noname noname; noname noname … noname noname; … ; noname noname … noname noname; noname noname … noname noname], [noname noname … noname noname; noname noname … noname noname; … ; noname noname … noname noname; noname noname … noname noname], [noname noname … noname noname; noname noname … noname noname; … ; noname noname … noname noname; noname noname … noname noname], [noname noname … noname noname; noname noname … noname noname; … ; noname noname … noname noname; noname noname … noname noname], [noname noname … noname noname; noname noname … noname noname; … ; noname noname … noname noname; noname noname … noname noname], [noname], [noname], [noname], [noname noname; noname noname], [noname]])

In [10]:

using JuMP
using SCS
λ = Ps = warm = nothing

Solving the problem¶

Depending on the actual problem one may need to tweak the parameters given to the solver:

eps sets the requested accuracy
max_iters sets the number of iterations to run before solver gives up
alpha is a parameter ( $\alpha \in (0,2)$ ) which determines the rate of convergence at the cost of the accuracy
acceleration_lookback: if you experience numerical instability in scs log should be changed to 1 (at the cost of rate of convergence).

The parameters below should be enough to obtain a decent solution for $\operatorname{SL}(4, \mathbb{Z}), \operatorname{SL}(5, \mathbb{Z})$ .
For $\operatorname{SL}(3, \mathbb{Z})$ approximately 1_000_000 of iterations is required; in this case by changing UB to $0.15$ (above) a much faster convergence can be observed.

In [11]:

with_SCS = with_optimizer(SCS.Optimizer, 
    linear_solver=SCS.Direct, 
    eps=3e-13,
    max_iters=10000,
    alpha=1.5,
    acceleration_lookback=10,
    warm_start=true)

status, warm = PropertyT.solve(SDP_problem, with_SCS, warm);

----------------------------------------------------------------------------
	SCS v2.0.2 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012-2017
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 130382
eps = 3.00e-13, alpha = 1.50, max_iters = 10000, normalize = 1, scale = 1.00
acceleration_lookback = 10, rho_x = 1.00e-03
Variables n = 1388, constraints m = 1946
Cones:	primal zero / dual free vars: 1196
	linear vars: 1
	sd vars: 749, sd blks: 14
Setup time: 8.46e-02s
SCS using variable warm-starting
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 1.02e+00  1.23e+00  9.75e-01 -3.44e+01  4.76e+00  0.00e+00  2.96e-03 
   100| 1.68e-04  1.10e-03  1.18e-03 -8.19e-01 -8.16e-01  2.95e-16  1.65e-01 
   200| 1.86e-04  4.32e-04  1.38e-04 -8.19e-01 -8.20e-01  5.89e-16  3.40e-01 
   300| 6.04e-05  5.47e-04  1.30e-04 -8.19e-01 -8.19e-01  2.39e-17  5.22e-01 
   400| 1.35e-04  4.46e-04  2.68e-04 -8.19e-01 -8.20e-01  2.06e-15  6.93e-01 
   500| 1.35e-04  4.21e-04  3.66e-04 -8.19e-01 -8.20e-01  2.36e-15  8.59e-01 
   600| 4.90e-05  2.34e-04  2.32e-04 -8.20e-01 -8.19e-01  3.42e-16  1.02e+00 
   700| 2.12e-05  1.74e-04  1.49e-04 -8.20e-01 -8.19e-01  1.51e-15  1.19e+00 
   800| 6.12e-05  2.23e-04  2.15e-04 -8.20e-01 -8.19e-01  4.95e-16  1.37e+00 
   900| 3.99e-05  1.79e-04  3.29e-05 -8.20e-01 -8.20e-01  5.59e-16  1.54e+00 
  1000| 1.30e-05  1.93e-04  1.04e-04 -8.20e-01 -8.20e-01  1.37e-15  1.72e+00 
  1100| 2.31e-05  2.01e-04  9.26e-05 -8.20e-01 -8.19e-01  2.58e-16  1.89e+00 
  1200| 2.15e-05  1.68e-04  2.12e-04 -8.20e-01 -8.19e-01  1.20e-15  2.07e+00 
  1300| 2.32e-05  1.61e-04  1.64e-04 -8.20e-01 -8.19e-01  5.90e-17  2.24e+00 
  1400| 1.19e-05  1.13e-04  1.75e-04 -8.20e-01 -8.19e-01  5.17e-16  2.42e+00 
  1500| 8.01e-06  1.33e-04  8.49e-05 -8.20e-01 -8.20e-01  3.33e-16  2.59e+00 
  1600| 3.49e-05  1.50e-04  1.04e-04 -8.20e-01 -8.19e-01  3.73e-16  2.76e+00 
  1700| 2.28e-05  8.89e-05  5.56e-05 -8.20e-01 -8.20e-01  2.86e-16  2.92e+00 
  1800| 1.54e-05  9.31e-05  1.06e-04 -8.20e-01 -8.20e-01  4.60e-16  3.09e+00 
  1900| 9.36e-06  1.25e-04  4.95e-06 -8.20e-01 -8.20e-01  3.36e-15  3.26e+00 
  2000| 3.71e-05  8.99e-05  2.00e-05 -8.20e-01 -8.20e-01  3.53e-16  3.46e+00 
  2100| 7.38e-06  7.97e-05  5.16e-05 -8.20e-01 -8.20e-01  1.92e-15  3.64e+00 
  2200| 1.08e-05  7.21e-05  3.00e-05 -8.20e-01 -8.20e-01  8.14e-16  3.81e+00 
  2300| 6.20e-06  9.69e-05  5.00e-05 -8.20e-01 -8.20e-01  1.36e-15  3.98e+00 
  2400| 8.47e-06  8.09e-05  6.05e-05 -8.20e-01 -8.20e-01  2.48e-15  4.15e+00 
  2500| 2.39e-05  8.09e-05  3.22e-05 -8.20e-01 -8.20e-01  3.40e-15  4.33e+00 
  2600| 1.56e-05  5.56e-05  1.86e-05 -8.20e-01 -8.20e-01  1.13e-16  4.50e+00 
  2700| 6.45e-06  6.44e-05  7.46e-07 -8.20e-01 -8.20e-01  1.15e-16  4.67e+00 
  2800| 2.76e-06  4.58e-05  5.44e-05 -8.20e-01 -8.20e-01  3.67e-15  4.84e+00 
  2900| 3.68e-06  5.38e-05  3.21e-05 -8.20e-01 -8.20e-01  1.04e-15  5.02e+00 
  3000| 6.13e-06  5.96e-05  5.55e-05 -8.20e-01 -8.20e-01  2.67e-16  5.19e+00 
  3100| 1.15e-05  4.97e-05  2.32e-06 -8.20e-01 -8.20e-01  4.40e-16  5.37e+00 
  3200| 1.28e-05  3.79e-05  4.61e-05 -8.20e-01 -8.20e-01  2.80e-15  5.54e+00 
  3300| 1.01e-05  8.71e-05  3.11e-05 -8.20e-01 -8.20e-01  1.90e-15  5.73e+00 
  3400| 5.20e-06  4.60e-05  1.14e-05 -8.20e-01 -8.20e-01  1.77e-16  5.90e+00 
  3500| 2.60e-06  2.80e-05  1.28e-05 -8.20e-01 -8.20e-01  7.29e-16  6.06e+00 
  3600| 6.14e-06  4.07e-05  1.17e-05 -8.20e-01 -8.20e-01  1.41e-15  6.23e+00 
  3700| 2.59e-06  2.02e-05  7.04e-06 -8.20e-01 -8.20e-01  7.81e-16  6.40e+00 
  3800| 5.14e-06  3.43e-05  2.91e-05 -8.20e-01 -8.20e-01  3.08e-16  6.59e+00 
  3900| 1.77e-06  2.43e-05  3.56e-07 -8.20e-01 -8.20e-01  1.42e-15  6.75e+00 
  4000| 6.59e-06  2.22e-05  1.03e-05 -8.20e-01 -8.20e-01  2.65e-16  6.92e+00 
  4100| 2.81e-06  2.11e-05  1.35e-05 -8.20e-01 -8.20e-01  1.65e-15  7.09e+00 
  4200| 1.90e-06  2.91e-05  2.26e-05 -8.20e-01 -8.20e-01  1.54e-15  7.26e+00 
  4300| 4.63e-06  3.09e-05  1.42e-06 -8.20e-01 -8.20e-01  1.15e-15  7.43e+00 
  4400| 2.02e-06  1.69e-05  7.11e-06 -8.20e-01 -8.20e-01  8.27e-16  7.61e+00 
  4500| 9.67e-06  2.00e-05  5.91e-06 -8.20e-01 -8.20e-01  5.69e-16  7.77e+00 
  4600| 8.40e-07  4.23e-06  3.47e-06 -8.20e-01 -8.20e-01  2.54e-15  7.93e+00 
  4700| 1.33e-06  9.78e-06  4.42e-08 -8.20e-01 -8.20e-01  8.54e-16  8.09e+00 
  4800| 8.40e-07  2.94e-06  4.89e-06 -8.20e-01 -8.20e-01  2.71e-15  8.25e+00 
  4900| 5.06e-07  1.84e-06  1.09e-06 -8.20e-01 -8.20e-01  2.43e-15  8.40e+00 
  5000| 1.26e-06  3.09e-06  4.54e-06 -8.20e-01 -8.20e-01  5.48e-17  8.55e+00 
  5100| 1.15e-07  2.30e-07  2.20e-07 -8.20e-01 -8.20e-01  4.83e-15  8.71e+00 
  5200| 2.57e-08  2.21e-07  2.86e-09 -8.20e-01 -8.20e-01  3.43e-15  8.85e+00 
  5300| 1.07e-04  2.94e-04  2.86e-04 -8.20e-01 -8.21e-01  4.81e-12  9.00e+00 
  5400| 8.39e-09  1.30e-07  8.58e-08 -8.20e-01 -8.20e-01  1.72e-15  9.14e+00 
  5500| 1.12e-08  1.68e-07  9.26e-09 -8.20e-01 -8.20e-01  2.03e-15  9.30e+00 
  5600| 2.88e-06  1.85e-05  1.04e-05 -8.20e-01 -8.20e-01  6.67e-15  9.45e+00 
  5700| 2.10e-05  2.21e-04  7.80e-05 -8.20e-01 -8.20e-01  8.08e-15  9.65e+00 
  5800| 4.36e-07  1.32e-06  8.45e-07 -8.20e-01 -8.20e-01  1.36e-15  9.80e+00 
  5900| 3.68e-06  3.39e-05  1.88e-06 -8.20e-01 -8.20e-01  1.39e-15  9.95e+00 
  6000| 2.28e-08  2.05e-07  9.89e-08 -8.20e-01 -8.20e-01  1.26e-15  1.01e+01 
  6100| 4.73e-08  7.77e-08  3.29e-08 -8.20e-01 -8.20e-01  1.86e-15  1.03e+01 
  6200| 5.33e-09  1.65e-08  1.07e-08 -8.20e-01 -8.20e-01  2.61e-16  1.04e+01 
  6300| 5.08e-08  3.23e-07  2.44e-07 -8.20e-01 -8.20e-01  1.71e-14  1.05e+01 
  6400| 5.01e-14  8.48e-13  4.12e-13 -8.20e-01 -8.20e-01  4.17e-16  1.07e+01 
  6402| 4.42e-14  1.73e-13  6.97e-14 -8.20e-01 -8.20e-01  7.23e-16  1.07e+01 
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 1.07e+01s
	Lin-sys: nnz in L factor: 288182, avg solve time: 1.00e-03s
	Cones: avg projection time: 4.08e-04s
	Acceleration: avg step time: 1.74e-04s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 8.4711e-10, dist(y, K*) = 2.4476e-09, s'y/|s||y| = -3.2607e-15
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 4.4176e-14
dual res:   |A'y + c|_2 / (1 + |c|_2) = 1.7338e-13
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 6.9723e-14
----------------------------------------------------------------------------
c'x = -0.8201, -b'y = -0.8201
============================================================================

In [12]:

λ = value(SDP_problem[:λ])
Ps = [value.(P) for P in varP]
@show(status, λ);

status = OPTIMAL::TerminationStatusCode = 1
λ = 0.8200500000003405

Checking the solution¶

Now we reconstruct the solution to the original problem over $\mathbb{R} \operatorname{SL}(N,\mathbb{Z})$ , which essentially boils down to averaging the obtained solution over the orbits of wreath product action:

$Q=\frac{1}{|\Sigma|}\sum_{\sigma\in\Sigma}\sum_{\pi\in \widehat{\Sigma}} \dim{\pi}\cdot\sigma\left(U_{\pi}^T \sqrt{P_{\pi}} U_{\pi}\right).$

In [13]:

Qs = real.(sqrt.(Ps));
Q = PropertyT.reconstruct(Qs, orbit_data);

As explained in the paper the columns of the square-root of the solution matrix provide the coefficients for $\xi_i$ 's in basis E_R of the group ring. Below we compute the residual

$b = \left(x - \lambda\Delta\right) - \sum \xi_i^*\xi_i.$ As we do it in floating-point arithmetic, the result can't be taken seriously.

In [14]:

function SOS_residual(x::GroupRingElem, Q::Matrix)
    RG = parent(x)
    @time sos = PropertyT.compute_SOS(RG, Q);
    return x - sos
end

Out[14]:

SOS_residual (generic function with 1 method)

In [15]:

residual = SOS_residual(elt - λ*Δ, Q)
@show norm(residual, 1);

┌ Warning: Scalar and coeffs are in different rings! Promoting result to Float64
└ @ GroupRings /home/kalmar/.julia/packages/GroupRings/UwACc/src/GroupRings.jl:303
┌ Warning: Adding elements with different coefficient rings, Promoting result to Float64
└ @ GroupRings /home/kalmar/.julia/packages/GroupRings/UwACc/src/GroupRings.jl:341

  0.017873 seconds (977 allocations: 2.065 MiB)
norm(residual, 1) = 8.53743030863061e-9

Checking in interval arithmetic¶

In [16]:

using IntervalArithmetic
IntervalArithmetic.setrounding(Interval, :tight)
IntervalArithmetic.setformat(sigfigs=12);

Here we resort to interval arithmetic to provide certified upper and lower bounds on the norm of the residual.

We first change entries of Q to narrow intervals
We project columns of Q so that $0$ is in the sum of coefficients of each column (i.e. $\xi_i \in I \operatorname{SL}(N,\mathbb{Z})$ )
We compute the sum of squares and the $\ell_1$ -norm of the residual in the interval arithmetic.

The returned check_columns_augmentation is a boolean flag to detect if the projection was successful, i.e. if we can guarantee that each column of Q_aug can be represented by an element from the augmentation ideal. (If it were not successful, one may project Q = PropertyT.augIdproj(Q) in the floating point arithmetic prior to the cell below).

The resulting norm of the residual is guaranteed to be contained in the resulting interval. E.g. if each entry of Q were changed into an honest rational number and all the computations were carried out in rational arithmetic, the rational $\ell_1$ -norm will be contained in the interval $\ell_1$ -norm.

In [17]:

Q_aug, check_columns_augmentation = PropertyT.augIdproj(Interval, Q);
@assert check_columns_augmentation
elt_int = elt - @interval(λ)*Δ;
residual_int = SOS_residual(elt_int, Q_aug)
@show norm(residual_int, 1);

┌ Warning: Scalar and coeffs are in different rings! Promoting result to Interval{Float64}
└ @ GroupRings /home/kalmar/.julia/packages/GroupRings/UwACc/src/GroupRings.jl:303
┌ Warning: Adding elements with different coefficient rings, Promoting result to Interval{Float64}
└ @ GroupRings /home/kalmar/.julia/packages/GroupRings/UwACc/src/GroupRings.jl:341

  3.296487 seconds (976 allocations: 4.070 MiB)
norm(residual_int, 1) = [1.06387679475e-08, 1.09608901207e-08]

In [18]:

certified_λ = @interval(λ) - 2^2*norm(residual_int,1)

Out[18]:

[0.820049956156, 0.820049957446]

So $\operatorname{elt} - \lambda_0 \Delta \in \Sigma^2 I\operatorname{SL}(N, \mathbb{Z})$ , where as $\lambda_0$ we could take the left end of the above interval:

In [19]:

certified_λ.lo

Out[19]:

0.8200499561567799

In [20]:

using Dates
now()

Out[20]:

2019-07-05T22:44:09.949

In [21]:

versioninfo()

Julia Version 1.1.1
Commit 55e36cc308 (2019-05-16 04:10 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-6200U CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 2

In [ ]: