Replication Details for 1812.03456
Below are replication details for computations with special linear groups as described in Section 5.1 of paper On property (T) for $\operatorname{Aut}(F_n)$ and $\operatorname{SL}_n$. The content of the accompanying jupyter notebook is reproduced below.
For exact details of computations for $36(\operatorname{Adj}_5 + 2 \operatorname{Op}_5) - 50\Delta_5 \in I\operatorname{SAut}(F_5)$ see the document deposited with the dataset in zenodo data repository.
Table of Contents
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 julia
s 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}).$$
using Pkg
Pkg.activate("..")
using Dates
now()
2019-07-05T22:42:41.473
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
)
N = 4
G = MatrixAlgebra(zz, N)
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
.
S = PropertyT.generating_set(G)
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
.
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]
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.
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.
@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$.
if N == 3
k = 0
elseif N == 4
k = 1
elseif N == 5
k = 1.5
end
elt = AdjN + k*OpN;
elt.coeffs
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{aligned} \text{ maximize : } \quad & \lambda\ \text{under constraints : }\quad & 0 \leqslant \lambda \leqslant \operatorname{UB},\ & x - \lambda \Delta = \sum \xi_i^* \xi_i,\ \end{aligned}$$
# @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)
(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]])
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 accuracymax_iters
sets the number of iterations to run before solver gives upalpha
is a parameter ($\alpha \in (0,2)$) which determines the rate of convergence at the cost of the accuracyacceleration_lookback
: if you experience numerical instability in scs log should be changed to1
(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.
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
============================================================================
λ = 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).$$
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.
function SOS_residual(x::GroupRingElem, Q::Matrix)
RG = parent(x)
@time sos = PropertyT.compute_SOS(RG, Q);
return x - sos
end
SOS_residual (generic function with 1 method)
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
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.
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]
certified_λ = @interval(λ) - 2^2*norm(residual_int,1)
[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:
certified_λ.lo
0.8200499561567799
using Dates
now()
2019-07-05T22:44:09.949
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