Compare commits

..

8 Commits

8 changed files with 301 additions and 83 deletions

2
.gitignore vendored
View File

@ -11,3 +11,5 @@ SL*_*
*.gws
.*
tests*
*.py
*.pyc

View File

@ -1,4 +1,10 @@
This repository contains code for computations in [Certifying Numerical Estimates of Spectral Gaps](https://arxiv.org/abs/1703.09680).
# DEPRECATED!
This repository has not been updated for a while!
If You are interested in replicating results for [1712.07167](https://arxiv.org/abs/1712.07167) please check [these instruction](https://kalmar.faculty.wmi.amu.edu.pl/post/1712.07176/)
Also [this notebook](https://nbviewer.jupyter.org/gist/kalmarek/03510181bc1e7c98615e86e1ec580b2a) could be of some help. If everything else fails the [zenodo dataset](https://zenodo.org/record/1133440) should contain the last-resort instructions.
This repository contains some legacy code for computations in [Certifying Numerical Estimates of Spectral Gaps](https://arxiv.org/abs/1703.09680).
# Installing
To run the code You need `julia-v0.5` (should work on `v0.6`, but with warnings).

View File

@ -4,6 +4,7 @@ using PropertyT
using AbstractAlgebra
using Nemo
using Groups
using GroupRings
export PropertyTGroup, SymmetrizedGroup, GAPGroup,
SpecialLinearGroup,
@ -50,5 +51,6 @@ end
include("mappingclassgroup.jl")
include("higman.jl")
include("caprace.jl")
include("actions.jl")
end # of module PropertyTGroups

92
groups/actions.jl Normal file
View File

@ -0,0 +1,92 @@
function (p::perm)(A::GroupRingElem)
RG = parent(A)
result = zero(RG, eltype(A.coeffs))
for (idx, c) in enumerate(A.coeffs)
if c!= zero(eltype(A.coeffs))
result[p(RG.basis[idx])] = c
end
end
return result
end
###############################################################################
#
# Action of WreathProductElems on Nemo.MatElem
#
###############################################################################
function matrix_emb(n::DirectProductGroupElem, p::perm)
Id = parent(n.elts[1])()
elt = diagm([(-1)^(el == Id ? 0 : 1) for el in n.elts])
return elt[:, p.d]
end
function (g::WreathProductElem)(A::MatElem)
g_inv = inv(g)
G = matrix_emb(g.n, g_inv.p)
G_inv = matrix_emb(g_inv.n, g.p)
M = parent(A)
return M(G)*A*M(G_inv)
end
import Base.*
doc"""
*(x::AbstractAlgebra.MatElem, P::Generic.perm)
> Apply the pemutation $P$ to the rows of the matrix $x$ and return the result.
"""
function *(x::AbstractAlgebra.MatElem, P::Generic.perm)
z = similar(x)
m = rows(x)
n = cols(x)
for i = 1:m
for j = 1:n
z[i, j] = x[i,P[j]]
end
end
return z
end
function (p::perm)(A::MatElem)
length(p.d) == A.r == A.c || throw("Can't act via $p on matrix of size ($(A.r), $(A.c))")
return p*A*inv(p)
end
###############################################################################
#
# Action of WreathProductElems on AutGroupElem
#
###############################################################################
function AutFG_emb(A::AutGroup, g::WreathProductElem)
isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
parent(g).P.n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(g)) into $A")
elt = A()
Id = parent(g.n.elts[1])()
flips = Groups.AutSymbol[Groups.flip_autsymbol(i) for i in 1:length(g.p.d) if g.n.elts[i] != Id]
Groups.r_multiply!(elt, flips, reduced=false)
Groups.r_multiply!(elt, [Groups.perm_autsymbol(g.p)])
return elt
end
function AutFG_emb(A::AutGroup, p::perm)
isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
parent(p).n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(p)) into $A")
return A(Groups.perm_autsymbol(p))
end
function (g::WreathProductElem)(a::Groups.Automorphism)
A = parent(a)
g = AutFG_emb(A,g)
res = A()
Groups.r_multiply!(res, g.symbols, reduced=false)
Groups.r_multiply!(res, a.symbols, reduced=false)
Groups.r_multiply!(res, [inv(s) for s in reverse!(g.symbols)])
return res
end
function (p::perm)(a::Groups.Automorphism)
g = AutFG_emb(parent(a),p)
return g*a*inv(g)
end

View File

@ -19,41 +19,3 @@ end
function autS(G::SpecialAutomorphismGroup{N}) where N
return WreathProduct(PermutationGroup(2), PermutationGroup(N))
end
###############################################################################
#
# Action of WreathProductElems on AutGroupElem
#
###############################################################################
function AutFG_emb(A::AutGroup, g::WreathProductElem)
isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
parent(g).P.n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(g)) into $A")
elt = A()
Id = parent(g.n.elts[1])()
flips = Groups.AutSymbol[Groups.flip_autsymbol(i) for i in 1:length(g.p.d) if g.n.elts[i] != Id]
Groups.r_multiply!(elt, flips, reduced=false)
Groups.r_multiply!(elt, [Groups.perm_autsymbol(g.p)])
return elt
end
function AutFG_emb(A::AutGroup, p::perm)
isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
parent(p).n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(g)) into $A")
return A(Groups.perm_autsymbol(p))
end
function (g::WreathProductElem)(a::Groups.Automorphism)
A = parent(a)
g = AutFG_emb(A,g)
res = A()
Groups.r_multiply!(res, g.symbols, reduced=false)
Groups.r_multiply!(res, a.symbols, reduced=false)
Groups.r_multiply!(res, [inv(s) for s in reverse!(g.symbols)])
return res
end
function (p::perm)(a::Groups.Automorphism)
g = AutFG_emb(parent(a),p)
return g*a*inv(g)
end

View File

@ -60,46 +60,3 @@ end
function autS(G::SpecialLinearGroup{N}) where N
return WreathProduct(PermutationGroup(2), PermutationGroup(N))
end
###############################################################################
#
# Action of WreathProductElems on Nemo.MatElem
#
###############################################################################
function matrix_emb(n::DirectProductGroupElem, p::perm)
Id = parent(n.elts[1])()
elt = diagm([(-1)^(el == Id ? 0 : 1) for el in n.elts])
return elt[:, p.d]
end
function (g::WreathProductElem)(A::MatElem)
g_inv = inv(g)
G = matrix_emb(g.n, g_inv.p)
G_inv = matrix_emb(g_inv.n, g.p)
M = parent(A)
return M(G)*A*M(G_inv)
end
import Base.*
doc"""
*(x::AbstractAlgebra.MatElem, P::Generic.perm)
> Apply the pemutation $P$ to the rows of the matrix $x$ and return the result.
"""
function *(x::AbstractAlgebra.MatElem, P::Generic.perm)
z = similar(x)
m = rows(x)
n = cols(x)
for i = 1:m
for j = 1:n
z[i, j] = x[i,P[j]]
end
end
return z
end
function (p::perm)(A::MatElem)
length(p.d) == A.r == A.c || throw("Can't act via $p on matrix of size ($(A.r), $(A.c))")
return p*A*inv(p)
end

View File

@ -27,7 +27,7 @@ function Settings(Gr::PropertyTGroup, args, solver)
S = PropertyTGroups.generatingset(Gr)
sol = solver
ub = get(args,"upper_bound", Inf)
ub = get(args,"upper-bound", Inf)
tol = get(args,"tol", 1e-10)
ws = get(args, "warmstart", false)

View File

@ -0,0 +1,197 @@
using AbstractAlgebra
using Groups
using GroupRings
using PropertyT
using SCS
solver(tol, iterations) =
SCSSolver(linearsolver=SCS.Direct,
eps=tol, max_iters=iterations,
alpha=1.95, acceleration_lookback=1)
include("../main.jl")
using PropertyTGroups
args = Dict("SAut" => 5, "upper-bound" => 50.0, "radius" => 2, "nosymmetry"=>false, "tol"=>1e-12, "iterations" =>200000, "warmstart" => true)
Gr = PropertyTGroups.PropertyTGroup(args)
sett = PropertyT.Settings(Gr, args,
solver(args["tol"], args["iterations"]))
@show sett
fullpath = PropertyT.fullpath(sett)
isdir(fullpath) || mkpath(fullpath)
# setup_logging(PropertyT.filename(fullpath, :fulllog), :fulllog)
function small_generating_set(RG::GroupRing{AutGroup{N}}, n) where N
indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
rmuls = [Groups.rmul_autsymbol(i,j) for (i,j) in indexing]
lmuls = [Groups.lmul_autsymbol(i,j) for (i,j) in indexing]
gen_set = RG.group.([rmuls; lmuls])
return [gen_set; inv.(gen_set)]
end
function computeX(RG::GroupRing{AutGroup{N}}) where N
Tn = small_generating_set(RG, N-1)
= Int64
Δn = length(Tn)*one(RG, ) - RG(Tn, );
Alt_N = [g for g in elements(PermutationGroup(N)) if parity(g) == 0]
@time X = sum(σ(Δn)*sum(τ(Δn) for τ Alt_N if τ σ) for σ in Alt_N);
return X
end
function Sq(RG::GroupRing{AutGroup{N}}) where N
T2 = small_generating_set(RG, 2)
= Int64
Δ₂ = length(T2)*one(RG, ) - RG(T2, );
Alt_N = [g for g in elements(PermutationGroup(N)) if parity(g) == 0]
elt = sum(σ(Δ₂)^2 for σ in Alt_N)
return elt
end
function Adj(RG::GroupRing{AutGroup{N}}) where N
T2 = small_generating_set(RG, 2)
= Int64
Δ₂ = length(T2)*one(RG, ) - RG(T2, );
Alt_N = [g for g in elements(PermutationGroup(N)) if parity(g) == 0]
adj(σ::perm, τ::perm, i=1, j=2) = Set([σ[i], σ[j]]) Set([τ[i], τ[j]])
adj(σ::perm) = [τ for τ in Alt_N if length(adj(σ, τ)) == 1]
@time elt = sum(σ(Δ₂)*sum(τ(Δ₂) for τ in adj(σ)) for σ in Alt_N);
# return RG(elt.coeffs÷factorial(N-2)^2)
return elt
end
function Op(RG::GroupRing{AutGroup{N}}) where N
T2 = small_generating_set(RG, 2)
= Int64
Δ₂ = length(T2)*one(RG, ) - RG(T2, );
Alt_N = [g for g in elements(PermutationGroup(N)) if parity(g) == 0]
adj(σ::perm, τ::perm, i=1, j=2) = Set([σ[i], σ[j]]) Set([τ[i], τ[j]])
adj(σ::perm) = [τ for τ in Alt_N if length(adj(σ, τ)) == 0]
@time elt = sum(σ(Δ₂)*sum(τ(Δ₂) for τ in adj(σ)) for σ in Alt_N);
# return RG(elt.coeffs÷factorial(N-2)^2)
return elt
end
const ELT_FILE = joinpath(dirname(PropertyT.filename(sett, )), "SqAdjOp_coeffs.jld")
const WARMSTART_FILE = PropertyT.filename(sett, :warmstart)
if isfile(PropertyT.filename(sett,)) && isfile(ELT_FILE) &&
isfile(PropertyT.filename(sett, :OrbitData))
# cached
Δ = PropertyT.loadGRElem(PropertyT.filename(sett,), sett.G)
RG = parent(Δ)
orbit_data = load(PropertyT.filename(sett, :OrbitData), "OrbitData")
sq_c, adj_c, op_c = load(ELT_FILE, "Sq", "Adj", "Op")
# elt = ELT_FILE, sett.G)
sq = GroupRingElem(sq_c, RG)
adj = GroupRingElem(adj_c, RG)
op = GroupRingElem(op_c, RG);
else
info("Compute Laplacian")
Δ = PropertyT.Laplacian(sett.S, sett.radius)
RG = parent(Δ)
info("Compute Sq, Adj, Op")
@time sq, adj, op = Sq(RG), Adj(RG), Op(RG)
PropertyT.saveGRElem(PropertyT.filename(sett, ), Δ)
save(ELT_FILE, "Sq", sq.coeffs, "Adj", adj.coeffs, "Op", op.coeffs)
info("Compute OrbitData")
if !isfile(PropertyT.filename(sett, :OrbitData))
orbit_data = PropertyT.OrbitData(parent(Y), sett.autS)
save(PropertyT.filename(sett, :OrbitData), "OrbitData", orbit_data)
else
orbit_data = load(PropertyT.filename(sett, :OrbitData), "OrbitData")
end
end;
orbit_data = PropertyT.decimate(orbit_data);
elt = adj+2op;
const SOLUTION_FILE = PropertyT.filename(sett, :solution)
if !isfile(SOLUTION_FILE)
SDP_problem, varλ, varP = PropertyT.SOS_problem(elt, Δ, orbit_data; upper_bound=sett.upper_bound)
begin
using SCS
scs_solver = SCS.SCSSolver(linearsolver=SCS.Direct,
eps=sett.tol,
max_iters=args["iterations"],
alpha=1.95,
acceleration_lookback=1)
JuMP.setsolver(SDP_problem, scs_solver)
end
λ = Ps = nothing
ws = PropertyT.warmstart(sett)
# using ProgressMeter
# @showprogress "Running SDP optimization step... " for i in 1:args["repetitions"]
while true
λ, Ps, ws = PropertyT.solve(PropertyT.filename(sett, :solverlog),
SDP_problem, varλ, varP, ws);
if all((!isnan).(ws[1]))
save(WARMSTART_FILE, "warmstart", ws, "λ", λ, "Ps", Ps)
save(WARMSTART_FILE[1:end-4]*"_$(now()).jld", "warmstart", ws, "λ", λ, "Ps", Ps)
else
warn("No valid solution was saved!")
end
end
info("Reconstructing P...")
@time P = PropertyT.reconstruct(Ps, orbit_data);
save(SOLUTION_FILE, "λ", λ, "P", P)
end
P, λ = load(SOLUTION_FILE, "P", "λ")
@show λ;
@time const Q = real(sqrtm(P));
function SOS_residual(eoi::GroupRingElem, Q::Matrix)
RG = parent(eoi)
@time sos = PropertyT.compute_SOS(RG, Q);
return eoi - sos
end
info("Floating Point arithmetic:")
EOI = elt - λ*Δ
b = SOS_residual(EOI, Q)
@show norm(b, 1);
info("Interval arithmetic:")
using IntervalArithmetic
Qint = PropertyT.augIdproj(Q);
@assert all([zero(eltype(Q)) in sum(view(Qint, :, i)) for i in 1:size(Qint, 2)])
EOI_int = elt - @interval(λ)*Δ;
Q_int = PropertyT.augIdproj(Q);
@assert all([zero(eltype(Q)) in sum(view(Q_int, :, i)) for i in 1:size(Q_int, 2)])
b_int = SOS_residual(EOI_int, Q_int)
@show norm(b_int, 1);
info("λ is certified to be > ", (@interval(λ) - 2^2*norm(b_int,1)).lo)