diff --git a/test/1703.09680.jl b/test/1703.09680.jl index cbed108..75e76c1 100644 --- a/test/1703.09680.jl +++ b/test/1703.09680.jl @@ -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 diff --git a/test/1712.07167.jl b/test/1712.07167.jl index 8c9dbd0..96ea904 100644 --- a/test/1712.07167.jl +++ b/test/1712.07167.jl @@ -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 diff --git a/test/1812.03456.jl b/test/1812.03456.jl index bc31c74..c3dbe5c 100644 --- a/test/1812.03456.jl +++ b/test/1812.03456.jl @@ -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 diff --git a/test/actions.jl b/test/actions.jl index 2cee0a7..af7d123 100644 --- a/test/actions.jl +++ b/test/actions.jl @@ -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(Γ)) + + RΓ = SymbolicWedderburn._group_algebra(Γ) + + @time mps, simple = + SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) + @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(Γ)) + + RΓ = SymbolicWedderburn._group_algebra(Γ) + + @time mps, simple = + SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) + @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(Γ)) + + RΓ = SymbolicWedderburn._group_algebra(Γ) + + @time mps, simple = + SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) + @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(Γ)) + + RΓ = SymbolicWedderburn._group_algebra(Γ) + + @time mps, simple = + SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) + @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 diff --git a/test/graded_adj.jl b/test/graded_adj.jl new file mode 100644 index 0000000..117b1dd --- /dev/null +++ b/test/graded_adj.jl @@ -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 + diff --git a/test/optimizers.jl b/test/optimizers.jl new file mode 100644 index 0000000..f415de7 --- /dev/null +++ b/test/optimizers.jl @@ -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 diff --git a/test/runtests.jl b/test/runtests.jl index dfb0413..daf0aa3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -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