reorganize/rewrite tests!

This commit is contained in:
Marek Kaluba 2022-11-07 18:45:12 +01:00
parent 7907137fb5
commit a080ae74c3
No known key found for this signature in database
GPG Key ID: 8BF1A3855328FC15
7 changed files with 871 additions and 387 deletions

View File

@ -1,51 +1,203 @@
function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer)
@time sos_problem =
PropertyT.sos_problem_primal(elt, unit, upper_bound=upper_bound)
status, _ = PropertyT.solve(sos_problem, optimizer)
P = JuMP.value.(sos_problem[:P])
Q = real.(sqrt(P))
certified, λ_cert = PropertyT.certify_solution(
elt,
unit,
JuMP.objective_value(sos_problem),
Q,
halfradius=halfradius,
)
return status, certified, λ_cert
end
@testset "1703.09680 Examples" begin
@testset "SL(2,Z)" begin
N = 2
G = MatrixAlgebra(zz, N)
S = PropertyT.generating_set(G)
G = MatrixGroups.SpecialLinearGroup{N}(Int8)
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
rm("SL($N,Z)", recursive=true, force=true)
sett = PropertyT.Settings("SL($N,Z)", G, S, with_SCS(20000, accel=20); upper_bound=0.1)
Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S)
end
@info sett
elt = Δ^2
unit = Δ
ub = 0.1
λ = PropertyT.spectral_gap(sett)
@test λ < 0.0
@test PropertyT.interpret_results(sett, λ) == false
status, certified, λ = check_positivity(
elt,
unit,
upper_bound=ub,
halfradius=2,
optimizer=scs_optimizer(
eps=1e-10,
max_iters=5_000,
accel=50,
alpha=1.9,
)
)
@test status == JuMP.ALMOST_OPTIMAL
@test !certified
@test λ < 0
end
@testset "SL(3,Z)" begin
@testset "SL(3,F₅)" begin
N = 3
G = MatrixAlgebra(zz, N)
S = PropertyT.generating_set(G)
G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{5})
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
rm("SL($N,Z)", recursive=true, force=true)
sett = PropertyT.Settings("SL($N,Z)", G, S, with_SCS(1000, accel=20); upper_bound=0.1)
Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S)
end
@info sett
elt = Δ^2
unit = Δ
ub = 1.01 # 1.5
λ = PropertyT.spectral_gap(sett)
@test λ > 0.099
@test PropertyT.interpret_results(sett, λ) == true
status, certified, λ = check_positivity(
elt,
unit,
upper_bound=ub,
halfradius=2,
optimizer=scs_optimizer(
eps=1e-10,
max_iters=5_000,
accel=50,
alpha=1.9,
)
)
@test PropertyT.check_property_T(sett) == true #second run should be fast
@test status == JuMP.OPTIMAL
@test certified
@test λ > 1
end
@testset "SAut(F₂)" begin
N = 2
G = SAut(FreeGroup(N))
S = PropertyT.generating_set(G)
G = SpecialAutomorphismGroup(FreeGroup(N))
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
rm("SAut(F$N)", recursive=true, force=true)
sett = PropertyT.Settings("SAut(F$N)", G, S, with_SCS(10000);
upper_bound=0.15)
Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S)
end
@info sett
elt = Δ^2
unit = Δ
ub = 0.1
λ = PropertyT.spectral_gap(sett)
@test λ < 0.0
@test PropertyT.interpret_results(sett, λ) == false
status, certified, λ = check_positivity(
elt,
unit,
upper_bound=0.1,
halfradius=2,
optimizer=scs_optimizer(
eps=1e-10,
max_iters=5_000,
accel=50,
alpha=1.9,
)
)
@test status == JuMP.ALMOST_OPTIMAL
@test λ < 0
@test !certified
end
@testset "SL(3,Z) has (T)" begin
n = 3
SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true)
Δ = RSL(length(S)) - sum(RSL(s) for s in S)
@testset "basic formulation" begin
elt = Δ^2
unit = Δ
ub = 0.1
opt_problem = PropertyT.sos_problem_primal(
elt,
unit,
upper_bound=ub,
augmented=false,
)
status, _ = PropertyT.solve(
opt_problem,
cosmo_optimizer(
eps=1e-10,
max_iters=10_000,
accel=0,
alpha=1.5,
),
)
@test status == JuMP.OPTIMAL
λ = JuMP.value(opt_problem[])
@test λ > 0.09
Q = real.(sqrt(JuMP.value.(opt_problem[:P])))
certified, λ_cert = PropertyT.certify_solution(
elt,
unit,
λ,
Q,
halfradius=2,
augmented=false,
)
@test certified
@test isapprox(λ_cert, λ, rtol=1e-5)
end
@testset "augmented formulation" begin
elt = Δ^2
unit = Δ
ub = 0.20001 # Inf
opt_problem = PropertyT.sos_problem_primal(
elt,
unit,
upper_bound=ub,
augmented=true,
)
status, _ = PropertyT.solve(
opt_problem,
scs_optimizer(
eps=1e-10,
max_iters=10_000,
accel=-10,
alpha=1.5,
),
)
@test status == JuMP.OPTIMAL
λ = JuMP.value(opt_problem[])
Q = real.(sqrt(JuMP.value.(opt_problem[:P])))
certified, λ_cert = PropertyT.certify_solution(
elt,
unit,
λ,
Q,
halfradius=2,
augmented=true,
)
@test certified
@test isapprox(λ_cert, λ, rtol=1e-5)
@test λ_cert > 2 // 10
end
end
end

View File

@ -1,104 +1,225 @@
function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimizer)
@assert aug(elt) == aug(unit) == 0
@time sos_problem, Ps =
PropertyT.sos_problem_primal(elt, unit, wd, upper_bound=upper_bound)
@time status, _ = PropertyT.solve(sos_problem, optimizer)
Q = let Ps = Ps
flPs = [real.(sqrt(JuMP.value.(P))) for P in Ps]
PropertyT.reconstruct(flPs, wd)
end
λ = JuMP.value(sos_problem[])
sos = let RG = parent(elt), Q = Q
z = zeros(eltype(Q), length(basis(RG)))
res = AlgebraElement(z, RG)
cnstrs = PropertyT.constraints(basis(RG), RG.mstructure, augmented=true)
PropertyT._cnstr_sos!(res, Q, cnstrs)
end
residual = elt - λ * unit - sos
λ_fl = PropertyT.sufficient_λ(residual, λ, halfradius=2)
λ_fl < 0 && return status, false, λ_fl
sos = let RG = parent(elt), Q = [PropertyT.IntervalArithmetic.@interval(q) for q in Q]
z = zeros(eltype(Q), length(basis(RG)))
res = AlgebraElement(z, RG)
cnstrs = PropertyT.constraints(basis(RG), RG.mstructure, augmented=true)
PropertyT._cnstr_sos!(res, Q, cnstrs)
end
λ_int = PropertyT.IntervalArithmetic.@interval(λ)
residual_int = elt - λ_int * unit - sos
λ_int = PropertyT.sufficient_λ(residual_int, λ_int, halfradius=2)
return status, λ_int > 0, PropertyT.IntervalArithmetic.inf(λ_int)
end
@testset "1712.07167 Examples" begin
@testset "SL(3,Z)" begin
N = 3
G = MatrixAlgebra(zz, N)
S = PropertyT.generating_set(G)
autS = WreathProduct(SymmetricGroup(2), SymmetricGroup(N))
NAME = "SL($N,Z)_orbit"
rm(NAME, recursive=true, force=true)
sett = PropertyT.Settings(NAME, G, S, autS, with_SCS(1000, accel=20);
upper_bound=0.27, force_compute=false)
@info sett
λ = PropertyT.spectral_gap(sett)
@test λ < 0.0
@test PropertyT.interpret_results(sett, λ) == false
# second run just checks the solution due to force_compute=false above
@test λ == PropertyT.spectral_gap(sett)
@test PropertyT.check_property_T(sett) == false
sett = PropertyT.Settings(NAME, G, S, autS, with_SCS(4000, accel=20);
upper_bound=0.27, force_compute=true)
@info sett
λ = PropertyT.spectral_gap(sett)
@test λ > 0.269999
@test PropertyT.interpret_results(sett, λ) == true
# this should be very fast due to warmstarting:
@test λ PropertyT.spectral_gap(sett) atol=1e-5
@test PropertyT.check_property_T(sett) == true
##########
# Symmetrizing by SymmetricGroup(3):
sett = PropertyT.Settings(NAME, G, S, SymmetricGroup(N), with_SCS(4000, accel=20, warm_start=false);
upper_bound=0.27, force_compute=true)
@info sett
λ = PropertyT.spectral_gap(sett)
@test λ > 0.269999
@test PropertyT.interpret_results(sett, λ) == true
end
@testset "SL(4,Z)" begin
N = 4
G = MatrixAlgebra(zz, N)
S = PropertyT.generating_set(G)
autS = WreathProduct(SymmetricGroup(2), SymmetricGroup(N))
NAME = "SL($N,Z)_orbit"
rm(NAME, recursive=true, force=true)
sett = PropertyT.Settings(NAME, G, S, autS, with_SCS(1000, accel=20);
upper_bound=1.3, force_compute=false)
@info sett
λ = PropertyT.spectral_gap(sett)
@test λ < 0.0
@test PropertyT.interpret_results(sett, λ) == false
# second run just checks the solution due to force_compute=false above
@test λ == PropertyT.spectral_gap(sett)
@test PropertyT.check_property_T(sett) == false
sett = PropertyT.Settings(NAME, G, S, autS, with_SCS(7000, accel=20, warm_start=true);
upper_bound=1.3, force_compute=true)
@info sett
@time λ = PropertyT.spectral_gap(sett)
@test λ > 1.2999
@test PropertyT.interpret_results(sett, λ) == true
# this should be very fast due to warmstarting:
@time @test λ PropertyT.spectral_gap(sett) atol=1e-5
@time @test PropertyT.check_property_T(sett) == true
end
@testset "SAut(F₃)" begin
N = 3
G = SAut(FreeGroup(N))
S = PropertyT.generating_set(G)
autS = WreathProduct(SymmetricGroup(2), SymmetricGroup(N))
G = SpecialAutomorphismGroup(FreeGroup(N))
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
NAME = "SAut(F$N)_orbit"
P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1)))
Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
act = PropertyT.action_by_conjugation(G, Σ)
wd = WedderburnDecomposition(
Float64,
Σ,
act,
basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
)
rm(NAME, recursive=true, force=true)
Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S)
end
sett = PropertyT.Settings(NAME, G, S, autS, with_SCS(1000);
upper_bound=0.15)
elt = Δ^2
unit = Δ
ub = Inf
@info sett
status, certified, λ_cert = check_positivity(
elt,
unit,
wd,
upper_bound=ub,
halfradius=2,
optimizer=cosmo_optimizer(
eps=1e-7,
max_iters=10_000,
accel=50,
alpha=1.9,
),
)
@test PropertyT.check_property_T(sett) == false
@test status == JuMP.OPTIMAL
@test !certified
@test λ_cert < 0
end
@testset "SL(3,Z) has (T)" begin
n = 3
SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true)
Δ = RSL(length(S)) - sum(RSL(s) for s in S)
@testset "Wedderburn formulation" begin
P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
act = PropertyT.action_by_conjugation(SL, Σ)
wd = WedderburnDecomposition(
Rational{Int},
Σ,
act,
basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
)
elt = Δ^2
unit = Δ
ub = 0.2801
@test_throws ErrorException PropertyT.sos_problem_primal(
elt,
unit,
wd,
upper_bound=ub,
augmented=false,
)
wdfl = SymbolicWedderburn.WedderburnDecomposition(
Float64,
Σ,
act,
basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
)
model, varP = PropertyT.sos_problem_primal(
elt,
unit,
wdfl,
upper_bound=ub,
augmented=false,
)
status, warm = PropertyT.solve(
model,
scs_optimizer(
eps=1e-10,
max_iters=20_000,
accel=-20,
alpha=1.2,
),
)
Q = @time let varP = varP
Qs = map(varP) do P
real.(sqrt(JuMP.value.(P)))
end
PropertyT.reconstruct(Qs, wdfl)
end
λ = JuMP.value(model[])
sos = PropertyT.compute_sos(parent(elt), Q; augmented=false)
certified, λ_cert = PropertyT.certify_solution(
elt,
unit,
λ,
Q,
halfradius=2,
augmented=false,
)
@test certified
@test λ_cert >= 28 // 100
end
@testset "augmented Wedderburn formulation" begin
elt = Δ^2
unit = Δ
ub = Inf
P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
act = PropertyT.action_by_conjugation(SL, Σ)
wdfl = SymbolicWedderburn.WedderburnDecomposition(
Float64,
Σ,
act,
basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
)
opt_problem, varP = PropertyT.sos_problem_primal(
elt,
unit,
wdfl,
upper_bound=ub,
# augmented = true # since both elt and unit are augmented
)
status, _ = PropertyT.solve(
opt_problem,
scs_optimizer(
eps=1e-8,
max_iters=20_000,
accel=0,
alpha=1.9,
),
)
@test status == JuMP.OPTIMAL
Q = @time let varP = varP
Qs = map(varP) do P
real.(sqrt(JuMP.value.(P)))
end
PropertyT.reconstruct(Qs, wdfl)
end
certified, λ_cert = PropertyT.certify_solution(
elt,
unit,
JuMP.objective_value(opt_problem),
Q,
halfradius=2,
# augmented = true # since both elt and unit are augmented
)
@test certified
@test λ_cert > 28 // 100
end
end
end

View File

@ -1,3 +1,5 @@
using SparseArrays
@testset "Sq, Adj, Op" begin
function isconstant_on_orbit(v, orb)
isempty(orb) && return true
@ -6,240 +8,237 @@
end
@testset "unit tests" begin
for N in [3,4]
M = MatrixAlgebra(zz, N)
@test PropertyT.EltaryMat(M, 1, 2) isa MatAlgElem
e12 = PropertyT.EltaryMat(M, 1, 2)
@test e12[1,2] == 1
@test inv(e12)[1,2] == -1
S = PropertyT.generating_set(M)
@test e12 S
@test length(PropertyT.generating_set(M)) == 2N*(N-1)
@test all(inv(s) S for s in S)
A = SAut(FreeGroup(N))
@test length(PropertyT.generating_set(A)) == 4N*(N-1)
S = PropertyT.generating_set(A)
@test all(inv(s) S for s in S)
end
@test PropertyT.isopposite(perm"(1,2,3)(4)", perm"(1,4,2)")
@test PropertyT.isadjacent(perm"(1,2,3)", perm"(1,2)(3)")
@test !PropertyT.isopposite(perm"(1,2,3)", perm"(1,2)(3)")
@test !PropertyT.isadjacent(perm"(1,4)", perm"(2,3)(4)")
@test isconstant_on_orbit([1,1,1,2,2], [2,3])
@test !isconstant_on_orbit([1,1,1,2,2], [2,3,4])
@test isconstant_on_orbit([1, 1, 1, 2, 2], [2, 3])
@test !isconstant_on_orbit([1, 1, 1, 2, 2], [2, 3, 4])
end
@testset "Sq, Adj, Op" begin
@testset "Sq, Adj, Op in SL(4,Z)" begin
N = 4
M = MatrixAlgebra(zz, N)
S = PropertyT.generating_set(M)
Δ = PropertyT.Laplacian(S, 2)
RG = parent(Δ)
G = MatrixGroups.SpecialLinearGroup{N}(Int8)
autS = WreathProduct(SymmetricGroup(2), SymmetricGroup(N))
orbits = PropertyT.orbit_decomposition(autS, RG.basis)
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
@test PropertyT.Sq(RG) isa GroupRingElem
sq = PropertyT.Sq(RG)
@test all(isconstant_on_orbit(sq, orb) for orb in orbits)
Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S)
end
@test PropertyT.Adj(RG) isa GroupRingElem
adj = PropertyT.Adj(RG)
@test all(isconstant_on_orbit(adj, orb) for orb in orbits)
P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1)))
Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
act = PropertyT.action_by_conjugation(G, Σ)
@test PropertyT.Op(RG) isa GroupRingElem
op = PropertyT.Op(RG)
@test all(isconstant_on_orbit(op, orb) for orb in orbits)
wd = WedderburnDecomposition(
Float64,
Σ,
act,
basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
)
ivs = SymbolicWedderburn.invariant_vectors(wd)
sq, adj, op = PropertyT.SqAdjOp(RG, N)
@test sq == PropertyT.Sq(RG)
@test adj == PropertyT.Adj(RG)
@test op == PropertyT.Op(RG)
e = one(M)
g = PropertyT.EltaryMat(M, 1,2)
h = PropertyT.EltaryMat(M, 1,3)
k = PropertyT.EltaryMat(M, 3,4)
@test all(
isconstant_on_orbit(sq, SparseArrays.nonzeroinds(iv)) for iv in ivs
)
edges = N*(N-1)÷2
@test sq[e] == 20*edges
@test all(
isconstant_on_orbit(adj, SparseArrays.nonzeroinds(iv)) for iv in ivs
)
@test all(
isconstant_on_orbit(op, SparseArrays.nonzeroinds(iv)) for iv in ivs
)
e = one(G)
g = G([alphabet(G)[MatrixGroups.ElementaryMatrix{N}(1, 2, Int8(1))]])
h = G([alphabet(G)[MatrixGroups.ElementaryMatrix{N}(1, 3, Int8(1))]])
k = G([alphabet(G)[MatrixGroups.ElementaryMatrix{N}(3, 4, Int8(1))]])
@test sq[e] == 120
@test sq[g] == sq[h] == -8
@test sq[g^2] == sq[h^2] == 1
@test sq[g*h] == sq[h*g] == 0
# @test adj[e] == ...
@test adj[g] == adj[h] # == ...
@test adj[e] == 384
@test adj[g] == adj[h] == -32
@test adj[g^2] == adj[h^2] == 0
@test adj[g*h] == adj[h*g] # == ...
@test adj[g*h] == adj[h*g] == 2
@test adj[k*h] == adj[h*k] == 1
# @test op[e] == ...
@test op[g] == op[h] # == ...
@test op[e] == 96
@test op[g] == op[h] == -8
@test op[g^2] == op[h^2] == 0
@test op[g*h] == op[h*g] == 0
@test op[g*k] == op[k*g] # == ...
@test op[g*k] == op[k*g] == 2
@test op[h*k] == op[k*h] == 0
end
@testset "SAut(F₃)" begin
n = 3
G = SpecialAutomorphismGroup(FreeGroup(n))
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
sq, adj, op = PropertyT.SqAdjOp(RG, n)
@test sq(one(G)) == 216
@test all(sq(g) == -16 for g in gens(G))
@test adj(one(G)) == 384
@test all(adj(g) == -32 for g in gens(G))
@test iszero(op)
end
end
@testset "1812.03456 examples" begin
function SOS_residual(x::GroupRingElem, Q::Matrix)
RG = parent(x)
@time sos = PropertyT.compute_SOS(RG, Q);
return x - sos
end
function check_positivity(elt, Δ, orbit_data, upper_bound, warm=nothing; with_solver=with_SCS(20_000, accel=10))
SDP_problem, varP = PropertyT.SOS_problem_primal(elt, Δ, orbit_data; upper_bound=upper_bound)
status, warm = PropertyT.solve(SDP_problem, with_solver, warm);
Base.Libc.flush_cstdio()
@info "Optimization status:" status
λ = value(SDP_problem[])
Ps = [value.(P) for P in varP]
Qs = real.(sqrt.(Ps));
Q = PropertyT.reconstruct(Qs, orbit_data);
b = SOS_residual(elt - λ*Δ, Q)
return b, λ, warm
end
@testset "SL(3,Z)" begin
N = 3
halfradius = 2
M = MatrixAlgebra(zz, N)
S = PropertyT.generating_set(M)
Δ = PropertyT.Laplacian(S, halfradius)
RG = parent(Δ)
orbit_data = PropertyT.BlockDecomposition(RG, WreathProduct(SymmetricGroup(2), SymmetricGroup(N)))
orbit_data = PropertyT.decimate(orbit_data);
n = 3
G = MatrixGroups.SpecialLinearGroup{n}(Int8)
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
Δ = RG(length(S)) - sum(RG(s) for s in S)
P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
act = PropertyT.action_by_conjugation(G, Σ)
wd = SymbolicWedderburn.WedderburnDecomposition(
Float64,
Σ,
act,
basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
)
sq, adj, op = PropertyT.SqAdjOp(RG, n)
@testset "Sq₃ is SOS" begin
elt = PropertyT.Sq(RG)
UB = 0.05 # 0.105?
elt = sq
UB = Inf # λ ≈ 0.1040844
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
@test 2^2*norm(residual, 1) < 2λ/100
status, certified, λ_cert = check_positivity(
elt,
Δ,
wd,
upper_bound=UB,
halfradius=2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
)
@test status == JuMP.OPTIMAL
@test certified
@test λ_cert > 104 // 1000
end
@testset "Adj₃ is SOS" begin
elt = PropertyT.Adj(RG)
UB = 0.1 # 0.157?
elt = adj
UB = Inf # λ ≈ 0.15858018
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) < λ
@test 2^2*norm(residual, 1) < λ/100
status, certified, λ_cert = check_positivity(
elt,
Δ,
wd,
upper_bound=UB,
halfradius=2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
)
@test status == JuMP.OPTIMAL
@test certified
@test λ_cert > 1585 // 10000
end
@testset "Op₃ is empty, so can not be certified" begin
elt = PropertyT.Op(RG)
elt = op
UB = Inf
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) > λ
status, certified, λ_cert = check_positivity(
elt,
Δ,
wd,
upper_bound=UB,
halfradius=2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
)
@test status == JuMP.OPTIMAL
@test !certified
@test λ_cert < 0
end
end
@testset "SL(4,Z)" begin
N = 4
halfradius = 2
M = MatrixAlgebra(zz, N)
S = PropertyT.generating_set(M)
Δ = PropertyT.Laplacian(S, halfradius)
RG = parent(Δ)
orbit_data = PropertyT.BlockDecomposition(RG, WreathProduct(SymmetricGroup(2), SymmetricGroup(N)))
orbit_data = PropertyT.decimate(orbit_data);
n = 4
@testset "Sq₄ is SOS" begin
elt = PropertyT.Sq(RG)
UB = 0.2 # 0.3172
G = MatrixGroups.SpecialLinearGroup{n}(Int8)
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
Δ = RG(length(S)) - sum(RG(s) for s in S)
@test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
@test 2^2*norm(residual, 1) < λ/100
P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
act = PropertyT.action_by_conjugation(G, Σ)
wd = SymbolicWedderburn.WedderburnDecomposition(
Float64,
Σ,
act,
basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
)
sq, adj, op = PropertyT.SqAdjOp(RG, n)
@testset "Sq is SOS" begin
elt = sq
UB = Inf # λ ≈ 0.31670
status, certified, λ_cert = check_positivity(
elt,
Δ,
wd,
upper_bound=UB,
halfradius=2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
)
@test status == JuMP.OPTIMAL
@test certified
@test λ_cert > 316 // 1000
end
@testset "Adj₄ is SOS" begin
elt = PropertyT.Adj(RG)
UB = 0.3 # 0.5459?
@testset "Adj is SOS" begin
elt = adj
UB = 0.541 # λ ≈ 0.545710
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
@test 2^2*norm(residual, 1) < λ/100
status, certified, λ_cert = check_positivity(
elt,
Δ,
wd,
upper_bound=UB,
halfradius=2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
)
@test status == JuMP.OPTIMAL
@test certified
@test λ_cert > 54 // 100
end
@testset "we can't cerify that Op₄ SOS" begin
elt = PropertyT.Op(RG)
UB = 2.0
@testset "Op is a sum of squares, but not an order unit" begin
elt = op
UB = Inf
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB,
with_solver=with_SCS(20_000, accel=10, eps=2e-10))
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) > λ # i.e. we can't certify positivity
end
@testset "Adj₄ + Op₄ is SOS" begin
elt = PropertyT.Adj(RG) + PropertyT.Op(RG)
UB = 0.6 # 0.82005
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
@test 2^2*norm(residual, 1) < λ/100
status, certified, λ_cert = check_positivity(
elt,
Δ,
wd,
upper_bound=UB,
halfradius=2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
)
@test status == JuMP.OPTIMAL
@test !certified
@test -1e-2 < λ_cert < 0
end
end
# @testset "Adj₄ + 100 Op₄ ∈ ISAut(F₄) is SOS" begin
# N = 4
# halfradius = 2
# M = SAut(FreeGroup(N))
# S = PropertyT.generating_set(M)
# Δ = PropertyT.Laplacian(S, halfradius)
# RG = parent(Δ)
# orbit_data = PropertyT.BlockDecomposition(RG, WreathProduct(SymmetricGroup(2), SymmetricGroup(N)))
# orbit_data = PropertyT.decimate(orbit_data);
#
# @time elt = PropertyT.Adj(RG) + 100PropertyT.Op(RG)
# UB = 0.05
#
# warm = nothing
#
# residual, λ, warm = check_positivity(elt, Δ, orbit_data, UB, warm, with_solver=with_SCS(20_000, accel=10))
# @info "obtained λ and residual" λ norm(residual, 1)
#
# @test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
# @test 2^2*norm(residual, 1) < λ/100
# end
end

View File

@ -1,101 +1,210 @@
@testset "actions on Group[Rings]" begin
Eij = PropertyT.EltaryMat
ssgs(M::MatAlgebra, i, j) = (S = [Eij(M, i, j), Eij(M, j, i)];
S = unique([S; inv.(S)]); S)
function ssgs(A::AutGroup, i, j)
rmuls = [Groups.transvection_R(i,j), Groups.transvection_R(j,i)]
lmuls = [Groups.transvection_L(i,j), Groups.transvection_L(j,i)]
gen_set = A.([rmuls; lmuls])
return unique([gen_set; inv.(gen_set)])
end
@testset "actions on SL(3,Z) and its group ring" begin
N = 3
halfradius = 2
M = MatrixAlgebra(zz, N)
S = PropertyT.generating_set(M)
E_R, sizes = Groups.wlmetric_ball(S, one(M), radius=2halfradius);
rdict = GroupRings.reverse_dict(E_R)
pm = GroupRings.create_pm(E_R, rdict, sizes[halfradius]; twisted=false);
RG = GroupRing(M, E_R, rdict, pm)
@testset "correctness of actions" begin
Δ = length(S)*RG(1) - sum(RG(s) for s in S)
@test Δ == PropertyT.spLaplacian(RG, S)
elt = S[5]
x = RG(1) - RG(elt)
elt2 = E_R[rand(sizes[1]:sizes[2])]
y = 2RG(elt2) - RG(elt)
for G in [SymmetricGroup(N), WreathProduct(SymmetricGroup(2), SymmetricGroup(N))]
@test all(one(M)^g == one(M) for g in G)
@test all(rdict[m^g] <= sizes[1] for g in G for m in S)
@test all(m^g*n^g == (m*n)^g for g in G for m in S for n in S)
@test all(Δ^g == Δ for g in G)
@test all(x^g == RG(1) - RG(elt^g) for g in G)
@test all(2RG(elt2^g) - RG(elt^g) == y^g for g in G)
function test_action(basis, group, act)
action = SymbolicWedderburn.action
return @testset "action definition" begin
@test all(basis) do b
e = one(group)
action(act, e, b) == b
end
end
@testset "small Laplacians" begin
for (i,j) in PropertyT.indexing(N)
Sij = ssgs(M, i,j)
Δij= PropertyT.spLaplacian(RG, Sij)
a = let a = rand(basis)
while isone(a)
a = rand(basis)
end
@assert !isone(a)
a
end
@test all(Δij^p == PropertyT.spLaplacian(RG, ssgs(M, p[i], p[j])) for p in SymmetricGroup(N))
g, h = let g_h = rand(group, 2)
while any(isone, g_h)
g_h = rand(group, 2)
end
@assert all(!isone, g_h)
g_h
end
@test all(Δij^g == PropertyT.spLaplacian(RG, ssgs(M, g.p[i], g.p[j])) for g in WreathProduct(SymmetricGroup(2), SymmetricGroup(N)))
action = SymbolicWedderburn.action
@test action(act, g, a) in basis
@test action(act, h, a) in basis
@test action(act, h, action(act, g, a)) == action(act, g * h, a)
@test all([(g, h) for g in group for h in group]) do (g, h)
x = action(act, h, action(act, g, a))
y = action(act, g * h, a)
x == y
end
if act isa SymbolicWedderburn.ByPermutations
@test all(basis) do b
action(act, g, b) basis && action(act, h, b) basis
end
end
end
end
@testset "actions on SAut(F_3) and its group ring" begin
N = 3
halfradius = 2
M = SAut(FreeGroup(N))
S = PropertyT.generating_set(M)
E_R, sizes = Groups.wlmetric_ball(S, one(M), radius=2halfradius);
## Testing
rdict = GroupRings.reverse_dict(E_R)
pm = GroupRings.create_pm(E_R, rdict, sizes[halfradius]; twisted=false);
RG = GroupRing(M, E_R, rdict, pm)
@testset "Actions on SL(3,)" begin
n = 3
@testset "correctness of actions" begin
SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true)
Δ = length(S)*RG(1) - sum(RG(s) for s in S)
@test Δ == PropertyT.spLaplacian(RG, S)
@testset "Permutation action" begin
elt = S[5]
x = RG(1) - RG(elt)
elt2 = E_R[rand(sizes[1]:sizes[2])]
y = 2RG(elt2) - RG(elt)
Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
ΓpA = PropertyT.action_by_conjugation(SL, Γ)
for G in [SymmetricGroup(N), WreathProduct(SymmetricGroup(2), SymmetricGroup(N))]
@test all(one(M)^g == one(M) for g in G)
@test all(rdict[m^g] <= sizes[1] for g in G for m in S)
@test all(m^g*n^g == (m*n)^g for g in G for m in S for n in S)
test_action(basis(RSL), Γ, ΓpA)
@test all(Δ^g == Δ for g in G)
@test all(x^g == RG(1) - RG(elt^g) for g in G)
@testset "mps is successful" begin
@test all(2RG(elt2^g) - RG(elt^g) == y^g for g in G)
charsΓ =
SymbolicWedderburn.Character{
Rational{Int},
}.(SymbolicWedderburn.irreducible_characters(Γ))
= SymbolicWedderburn._group_algebra(Γ)
@time mps, simple =
SymbolicWedderburn.minimal_projection_system(charsΓ, )
@test all(simple)
end
@testset "Wedderburn decomposition" begin
wd = SymbolicWedderburn.WedderburnDecomposition(
Rational{Int},
Γ,
ΓpA,
basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]])
)
@test length(invariant_vectors(wd)) == 918
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [40, 23, 18]
@test all(issimple, direct_summands(wd))
end
end
for (i,j) in PropertyT.indexing(N)
Sij = ssgs(M, i,j)
Δij= PropertyT.spLaplacian(RG, Sij)
@testset "Wreath action" begin
@test all(Δij^p == PropertyT.spLaplacian(RG, ssgs(M, p[i], p[j])) for p in SymmetricGroup(N))
Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
end
@test all(Δij^g == PropertyT.spLaplacian(RG, ssgs(M, g.p[i], g.p[j])) for g in WreathProduct(SymmetricGroup(2), SymmetricGroup(N)))
ΓpA = PropertyT.action_by_conjugation(SL, Γ)
test_action(basis(RSL), Γ, ΓpA)
@testset "mps is successful" begin
charsΓ =
SymbolicWedderburn.Character{
Rational{Int},
}.(SymbolicWedderburn.irreducible_characters(Γ))
= SymbolicWedderburn._group_algebra(Γ)
@time mps, simple =
SymbolicWedderburn.minimal_projection_system(charsΓ, )
@test all(simple)
end
@testset "Wedderburn decomposition" begin
wd = SymbolicWedderburn.WedderburnDecomposition(
Rational{Int},
Γ,
ΓpA,
basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]])
)
@test length(invariant_vectors(wd)) == 247
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [14, 9, 6, 14, 12]
@test all(issimple, direct_summands(wd))
end
end
end
@testset "Actions on SAut(F4)" begin
n = 4
SAutFn = SpecialAutomorphismGroup(FreeGroup(n))
RSAutFn, S, sizes = PropertyT.group_algebra(SAutFn, halfradius=1, twisted=true)
@testset "Permutation action" begin
Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
ΓpA = PropertyT.action_by_conjugation(SAutFn, Γ)
test_action(basis(RSAutFn), Γ, ΓpA)
@testset "mps is successful" begin
charsΓ =
SymbolicWedderburn.Character{
Rational{Int},
}.(SymbolicWedderburn.irreducible_characters(Γ))
= SymbolicWedderburn._group_algebra(Γ)
@time mps, simple =
SymbolicWedderburn.minimal_projection_system(charsΓ, )
@test all(simple)
end
@testset "Wedderburn decomposition" begin
wd = SymbolicWedderburn.WedderburnDecomposition(
Rational{Int},
Γ,
ΓpA,
basis(RSAutFn),
StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]])
)
@test length(invariant_vectors(wd)) == 93
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [4, 8, 5, 4]
@test all(issimple, direct_summands(wd))
end
end
@testset "Wreath action" begin
Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
end
ΓpA = PropertyT.action_by_conjugation(SAutFn, Γ)
test_action(basis(RSAutFn), Γ, ΓpA)
@testset "mps is successful" begin
charsΓ =
SymbolicWedderburn.Character{
Rational{Int},
}.(SymbolicWedderburn.irreducible_characters(Γ))
= SymbolicWedderburn._group_algebra(Γ)
@time mps, simple =
SymbolicWedderburn.minimal_projection_system(charsΓ, )
@test all(simple)
end
@testset "Wedderburn decomposition" begin
wd = SymbolicWedderburn.WedderburnDecomposition(
Rational{Int},
Γ,
ΓpA,
basis(RSAutFn),
StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]])
)
@test length(invariant_vectors(wd)) == 18
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [1, 1, 2, 2, 1, 2, 2, 1]
@test all(issimple, direct_summands(wd))
end
end
end

40
test/graded_adj.jl Normal file
View File

@ -0,0 +1,40 @@
@testset "Adj for SpN via grading" begin
genus = 3
halfradius = 2
SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true)
Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
Δ = RG(length(S)) - sum(RG(s) for s in S)
Δs = PropertyT.laplacians(
RG,
S,
x -> (gx = PropertyT.grading(ψ(x)); Set([gx, -gx])),
)
Δ, Δs
end
@testset "Adj correctness: genus=$genus" begin
all_subtypes = (
:A₁, :C₁, Symbol("A₁×A₁"), Symbol("C₁×C₁"), Symbol("A₁×C₁"), :A₂, :C₂
)
@test PropertyT.Adj(Δs, :A₂)[one(SpN)] == 384
@test iszero(PropertyT.Adj(Δs, Symbol("A₁×A₁")))
@test iszero(PropertyT.Adj(Δs, Symbol("C₁×C₁")))
@testset "divisibility by 16" begin
for subtype in all_subtypes
subtype in (:A₁, :C₁) && continue
@test isinteger(PropertyT.Adj(Δs, subtype)[one(SpN)] / 16)
end
end
@test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) == Δ^2
end
end

54
test/optimizers.jl Normal file
View File

@ -0,0 +1,54 @@
## Optimizers
import JuMP
import SCS
function scs_optimizer(;
accel=10,
alpha=1.5,
eps=1e-9,
max_iters=100_000,
verbose=true,
linear_solver=SCS.DirectSolver
)
return JuMP.optimizer_with_attributes(
SCS.Optimizer,
"acceleration_lookback" => accel,
"acceleration_interval" => 10,
"alpha" => alpha,
"eps_abs" => eps,
"eps_rel" => eps,
"linear_solver" => linear_solver,
"max_iters" => max_iters,
"warm_start" => true,
"verbose" => verbose,
)
end
import COSMO
function cosmo_optimizer(;
accel=15,
alpha=1.6,
eps=1e-9,
max_iters=100_000,
verbose=true,
verbose_timing=verbose,
decompose=false
)
return JuMP.optimizer_with_attributes(
COSMO.Optimizer,
"accelerator" => COSMO.with_options(
COSMO.AndersonAccelerator,
mem=max(accel, 2)
),
"alpha" => alpha,
"decompose" => decompose,
"eps_abs" => eps,
"eps_rel" => eps,
"max_iter" => max_iters,
"verbose" => verbose,
"verbose_timing" => verbose_timing,
"check_termination" => 200,
)
end

View File

@ -1,18 +1,27 @@
using Test
using LinearAlgebra, SparseArrays
using AbstractAlgebra, Groups, GroupRings
using LinearAlgebra
using SparseArrays
BLAS.set_num_threads(1)
ENV["OMP_NUM_THREADS"] = 4
using Groups
using Groups.GroupsCore
import Groups.MatrixGroups
using PropertyT
using JLD
using SymbolicWedderburn
using SymbolicWedderburn.StarAlgebras
using SymbolicWedderburn.PermutationGroups
using JuMP, SCS
include("optimizers.jl")
with_SCS(iters; accel=0, eps=1e-10, warm_start=true) =
with_optimizer(SCS.Optimizer,
linear_solver=SCS.DirectSolver, max_iters=iters,
acceleration_lookback=accel, eps=eps, warm_start=warm_start)
@testset "PropertyT" begin
include("1703.09680.jl")
include("actions.jl")
include("1712.07167.jl")
include("SOS_correctness.jl")
include("1812.03456.jl")
include("actions.jl")
include("1703.09680.jl")
include("1712.07167.jl")
include("1812.03456.jl")
include("graded_adj.jl")
end