kalmar
/
PropertyT.jl
mirror of https://github.com/kalmarek/PropertyT.jl.git
1
0
Fork 0
Code Issues Releases Wiki Activity

Merge pull request #11 from kalmarek/enh/G₂_here_we_come 

G₂ - here we come!
This commit is contained in:
Marek Kaluba 2023-05-09 12:37:56 +02:00 committed by GitHub
parent 0127d05594 d40a0fe117
commit 0908c8ea1c
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
28 changed files with 1437 additions and 577 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
.JuliaFormatter.toml Normal file
View File

@ -0,0 +1,9 @@
# Configuration file for JuliaFormatter.jl
# For more information, see: https://domluna.github.io/JuliaFormatter.jl/stable/config/
always_for_in = true
always_use_return = true
margin = 80
remove_extra_newlines = true
separate_kwargs_with_semicolon = true
short_to_long_function_def = true

6
Project.toml
View File

@ -6,6 +6,7 @@ version = "0.4.0"
[deps] [deps]
Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557" Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572" JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
@ -17,11 +18,12 @@ SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
COSMO = "0.8" COSMO = "0.8"
Groups = "0.7" Groups = "0.7"
IntervalArithmetic = "0.20" IntervalArithmetic = "0.20"
IntervalMatrices = "0.8"
JuMP = "1.3" JuMP = "1.3"
ProgressMeter = "1.7" ProgressMeter = "1.7"
SCS = "1.1.0" SCS = "1.1"
StaticArrays = "1" StaticArrays = "1"
SymbolicWedderburn = "0.3.2" SymbolicWedderburn = "0.3.4"
julia = "1.6" julia = "1.6"
[extras] [extras]

98
scripts/G₂_Adj.jl Normal file
View File

@ -0,0 +1,98 @@
using LinearAlgebra
BLAS.set_num_threads(8)
using MKL_jll
ENV["OMP_NUM_THREADS"] = 4
using Groups
import Groups.MatrixGroups
include(joinpath(@__DIR__, "../test/optimizers.jl"))
using PropertyT
using PropertyT.SymbolicWedderburn
using PropertyT.PermutationGroups
using PropertyT.StarAlgebras
include(joinpath(@__DIR__, "argparse.jl"))
include(joinpath(@__DIR__, "utils.jl"))
# const N = parsed_args["N"]
const HALFRADIUS = parsed_args["halfradius"]
const UPPER_BOUND = parsed_args["upper_bound"]
include(joinpath(@__DIR__, "./G₂_gens.jl"))
G, roots, Weyl = G₂_roots_weyl()
@info "Running Adj² - λ·Δ sum of squares decomposition for G₂"
@info "computing group algebra structure"
RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
@info "computing WedderburnDecomposition"
wd = let Σ = Weyl, RG = RG
act = PropertyT.AlphabetPermutation{eltype(Σ),Int64}(
Dict(g => PermutationGroups.perm(g) for g in Σ),
)
@time SymbolicWedderburn.WedderburnDecomposition(
Float64,
Σ,
act,
basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
semisimple = false,
)
end
@info wd
function desubscriptify(symbol::Symbol)
digits = [
Int(l) - 0x2080 for
l in reverse(string(symbol)) if 0 ≤ Int(l) - 0x2080 ≤ 9
]
res = 0
for (i, d) in enumerate(digits)
res += 10^(i - 1) * d
end
return res
end
function PropertyT.grading(g::MatrixGroups.MatrixElt, roots = roots)
id = desubscriptify(g.id)
return roots[id]
end
Δ = RG(length(S)) - sum(RG(s) for s in S)
Δs = PropertyT.laplacians(
RG,
S,
x -> (gx = PropertyT.grading(x); Set([gx, -gx])),
)
elt = PropertyT.Adj(Δs)
@assert elt == Δ^2 - PropertyT.Sq(Δs)
unit = Δ
@time model, varP = PropertyT.sos_problem_primal(
elt,
unit,
wd;
upper_bound = UPPER_BOUND,
augmented = true,
show_progress = true,
)
solve_in_loop(
model,
wd,
varP;
logdir = "./log/G2/r=$HALFRADIUS/Adj-$(UPPER_BOUND)Δ",
optimizer = scs_optimizer(;
linear_solver = SCS.MKLDirectSolver,
eps = 1e-9,
max_iters = 100_000,
accel = 50,
alpha = 1.95,
),
data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
)

308
scripts/G₂_gens.jl Normal file
View File

@ -0,0 +1,308 @@
#= GAP code to generate matrices
alg := SimpleLieAlgebra("G", 2, Rationals);
root_sys := RootSystem(alg);
pos_gens := PositiveRootVectors(root_sys);
pos_rts := PositiveRoots(root_sys);
neg_gens := NegativeRootVectors(root_sys);
neg_rts := NegativeRoots(root_sys);
alg_gens := ShallowCopy(pos_gens);;
Append(alg_gens, neg_gens);
grading := ShallowCopy(pos_rts);
Append(grading, neg_rts);
mats := List(alg_gens, x->AdjointMatrix(Basis(alg), x));
W := WeylGroup(root_sys);
PW := Action(W, grading, OnRight);
=#
using LinearAlgebra
function matrix_exp(M::AbstractMatrix{<:Integer})
res = zeros(Rational{eltype(M)}, size(M))
res += I
k = 0
expM = one(M)
while !iszero(expM)
k += 1
expM *= M
@. res += 1 // factorial(k) * expM
if k == 20
@warn "matrix exponential did not converge" norm(expM - exp(M))
break
end
end
@debug "matrix_exp converged after $k iterations"
return res
end
const gap_adj_mats = [
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -2],
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1],
[2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0],
[3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 1],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2],
[0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1],
[0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0, 0, 0],
[0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, -3, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0],
[0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -1],
[0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0],
],
]
function G₂_matrices_roots()
adj_mats = map(gap_adj_mats) do m
return hcat(m...)
end
adj_mats = filter!(!isdiag, adj_mats) # remove the ones from center
gens_mats = [convert(Matrix{Int}, matrix_exp(m')) for m in adj_mats]
#=
The roots from
G₂roots_gap = [
[2, -1], # α = e₁ - e₂
[-3, 2], # A = -α + β = -e₁ + 2e₂ - e₃
[-1, 1], # β = e₂ - e₃
[1, 0], # α + β = e₁ - e₃
[3, -1], # B = 2α + β = 2e₁ - e₂ - e₃
[0, 1], # A + B = α + 2β = e₁ + e₂ - 2e₃
[-2, 1], # -α
[3, -2], # -A
[1, -1], # -β
[-1, 0], # -α - β
[-3, 1], # -B
[0, -1], # -A - B
]
G₂roots_gap are the ones from cartan matrix. To obtain the standard
(hexagonal) picture map them by `T` defined as follows:
```julia
cartan = hcat(G₂roots_gap[1:2]...)
rot(α) = [cos(α) -sin(α); sin(α) cos(α)]
c₁ = [√2, 0]
c₂ = rot(5π / 6) * [√2, 0] * √3 # (= 1/2[√6, 1])
T = hcat(c₁, c₂) * inv(cartan)
```
By plotting one against the others (or by blind calculation) one can see
the following assignment. Here `⟨α, β⟩_ℤ = A₂` and `⟨A, B⟩_ℤ ≅ √3/√2 A₂`.
=#
e₁ = PropertyT.Roots.𝕖(3, 1)
e₂ = PropertyT.Roots.𝕖(3, 2)
e₃ = PropertyT.Roots.𝕖(3, 3)
α = e₁ - e₂
β = e₂ - e₃
A = -α + β
B = α + (α + β)
roots = [α, A, β, α + β, B, A + B, -α, -A, -β, -α - β, -B, -A - B]
return gens_mats, roots
end
function G₂_roots_weyl()
(mats, roots) = G₂_matrices_roots()
d = size(first(mats), 1)
G₂ = Groups.MatrixGroup{d}(mats)
m = Groups.gens(G₂)
σ = let w = m[1] * inv(m[7]) * m[1], m = union(m, inv.(m))
PermutationGroups.Perm([findfirst(==(inv(w) * x * w), m) for x in m])
end
τ = let w = m[2] * inv(m[8]) * m[2], m = union(m, inv.(m))
PermutationGroups.Perm([findfirst(==(inv(w) * x * w), m) for x in m])
end
W = PermGroup(σ, τ)
return G₂, roots, W
end

98
scripts/G₂_has_T.jl Normal file
View File

@ -0,0 +1,98 @@
using LinearAlgebra
BLAS.set_num_threads(8)
ENV["OMP_NUM_THREADS"] = 4
using Groups
import Groups.MatrixGroups
include(joinpath(@__DIR__, "../test/optimizers.jl"))
using PropertyT
using PropertyT.SymbolicWedderburn
using PropertyT.PermutationGroups
using PropertyT.StarAlgebras
include(joinpath(@__DIR__, "argparse.jl"))
include(joinpath(@__DIR__, "utils.jl"))
# const N = parsed_args["N"]
const HALFRADIUS = parsed_args["halfradius"]
const UPPER_BOUND = parsed_args["upper_bound"]
include(joinpath(@__DIR__, "./G₂_gens.jl"))
G, roots, Weyl = G₂_roots_weyl()
@info "Running Δ² - λ·Δ sum of squares decomposition for G₂"
@info "computing group algebra structure"
RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
@info "computing WedderburnDecomposition"
wd = let Σ = Weyl, RG = RG
act = PropertyT.AlphabetPermutation{eltype(Σ),Int64}(
Dict(g => PermutationGroups.perm(g) for g in Σ),
)
@time SymbolicWedderburn.WedderburnDecomposition(
Float64,
Σ,
act,
basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
semisimple = false,
)
end
@info wd
Δ = RG(length(S)) - sum(RG(s) for s in S)
elt = Δ^2
unit = Δ
@time model, varP = PropertyT.sos_problem_primal(
elt,
unit,
wd;
upper_bound = UPPER_BOUND,
augmented = true,
show_progress = false,
)
warm = nothing
status = JuMP.OPTIMIZE_NOT_CALLED
certified, λ = false, nothing
while status ≠ JuMP.OPTIMAL
@time status, warm = PropertyT.solve(
model,
scs_optimizer(;
eps = 1e-10,
max_iters = 20_000,
accel = 50,
alpha = 1.95,
),
warm,
)
@info "reconstructing the solution"
Q = @time let wd = wd, Ps = [JuMP.value.(P) for P in varP]
Qs = real.(sqrt.(Ps))
PropertyT.reconstruct(Qs, wd)
end
@info "certifying the solution"
@time certified, λ = PropertyT.certify_solution(
elt,
unit,
JuMP.objective_value(model),
Q;
halfradius = HALFRADIUS,
augmented = true,
)
end
if certified && λ > 0
Κ(λ, S) = round(sqrt(2λ / length(S)), Base.RoundDown; digits = 5)
@info "Certified result: G₂ has property (T):" N λ Κ(λ, S)
else
@info "Could NOT certify the result:" certified λ
end

30
scripts/SpN_Adj.jl
View File

@ -24,8 +24,7 @@ const GENUS = 2N
G = MatrixGroups.SymplecticGroup{GENUS}(Int8) G = MatrixGroups.SymplecticGroup{GENUS}(Int8)
RG, S, sizes = RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
@time PropertyT.group_algebra(G, halfradius=HALFRADIUS, twisted=true)
wd = let RG = RG, N = N wd = let RG = RG, N = N
G = StarAlgebras.object(RG) G = StarAlgebras.object(RG)
@ -42,6 +41,8 @@ wd = let RG = RG, N = N
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]), StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
) )
@info wdfl
wdfl
end end
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -58,23 +59,22 @@ unit = Δ
@time model, varP = PropertyT.sos_problem_primal( @time model, varP = PropertyT.sos_problem_primal(
elt, elt,
unit, unit,
wd, wd;
upper_bound=UPPER_BOUND, upper_bound = UPPER_BOUND,
augmented=true, augmented = true,
show_progress=true show_progress = true,
) )
solve_in_loop( solve_in_loop(
model, model,
wd, wd,
varP, varP;
logdir="./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-InfΔ", logdir = "./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-$(UPPER_BOUND)Δ",
optimizer=cosmo_optimizer( optimizer = cosmo_optimizer(;
eps=1e-10, eps = 1e-10,
max_iters=20_000, max_iters = 20_000,
accel=50, accel = 50,
alpha=1.95, alpha = 1.95,
), ),
data=(elt=elt, unit=unit, halfradius=HALFRADIUS) data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
) )

47
scripts/utils.jl
View File

@ -11,12 +11,19 @@ function get_solution(model)
return solution return solution
end end
function get_solution(model, wd, varP; logdir) function get_solution(model, wd, varP, eps = 1e-10)
λ = JuMP.value(model[:λ]) λ = JuMP.value(model[:λ])
Qs = [real.(sqrt(JuMP.value.(P))) for P in varP] @info "reconstructing the solution"
Q = PropertyT.reconstruct(Qs, wd) Q = @time let wd = wd, Ps = [JuMP.value.(P) for P in varP], eps = eps
PropertyT.__droptol!.(Ps, 100eps)
Qs = real.(sqrt.(Ps))
PropertyT.__droptol!.(Qs, eps)
PropertyT.reconstruct(Qs, wd)
end
solution = Dict(:λ => λ, :Q => Q) solution = Dict(:λ => λ, :Q => Q)
return solution return solution
end end
@ -42,22 +49,23 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
# logstream = current_logger().logger.stream # logstream = current_logger().logger.stream
# v = @ccall setvbuf(logstream.handle::Ptr{Cvoid}, C_NULL::Ptr{Cvoid}, 1::Cint, 0::Cint)::Cint # v = @ccall setvbuf(logstream.handle::Ptr{Cvoid}, C_NULL::Ptr{Cvoid}, 1::Cint, 0::Cint)::Cint
# @warn v # @warn v
status, warm = @time PropertyT.solve(log_file, model, optimizer, warm) status, warm =
@time PropertyT.solve(log_file, model, optimizer, warm)
solution = get_solution(model, args...; logdir=logdir)
solution[:warm] = warm
solution = get_solution(model, args...)
serialize(joinpath(logdir, "solution_$date.sjl"), solution) serialize(joinpath(logdir, "solution_$date.sjl"), solution)
serialize(joinpath(logdir, "solution.sjl"), solution) serialize(joinpath(logdir, "solution.sjl"), solution)
flag, λ_cert = open(log_file, append=true) do io solution[:warm] = warm
flag, λ_cert = open(log_file; append = true) do io
with_logger(SimpleLogger(io)) do with_logger(SimpleLogger(io)) do
PropertyT.certify_solution( return PropertyT.certify_solution(
data.elt, data.elt,
data.unit, data.unit,
solution[:λ], solution[:λ],
solution[:Q], solution[:Q];
halfradius=data.halfradius, halfradius = data.halfradius,
) )
end end
end end
@ -66,19 +74,20 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
end end
if flag == true && certified_λ ≥ 0 if flag == true && certified_λ ≥ 0
@info "Certification done with λ = $certified_λ" @info "Certification done with λ = $certified_λ" certified_λ status
return certified_λ return certified_λ
else else
rel_change = abs(certified_λ - old_lambda) / (abs(certified_λ) + abs(old_lambda)) rel_change =
@info "Certification failed with λ = $λ" certified_λ rel_change abs(certified_λ - old_lambda) /
(abs(certified_λ) + abs(old_lambda))
@info "Certification failed with λ = $λ" certified_λ rel_change status
if rel_change < 1e-9
@info "No progress detected, breaking" certified_λ rel_change status
break
end
end end
old_lambda = certified_λ old_lambda = certified_λ
if rel_change < 1e-9
@info "No progress detected, breaking"
break
end
end end
return status == JuMP.OPTIMAL ? old_lambda : NaN return status == JuMP.OPTIMAL ? old_lambda : NaN

11
src/PropertyT.jl
View File

@ -3,7 +3,6 @@ module PropertyT
using LinearAlgebra using LinearAlgebra
using SparseArrays using SparseArrays
using IntervalArithmetic
using JuMP using JuMP
using Groups using Groups
@ -26,17 +25,17 @@ include("gradings.jl")
include("actions/actions.jl") include("actions/actions.jl")
function group_algebra(G::Groups.Group, S=gens(G); halfradius::Integer, twisted::Bool) function group_algebra(G::Groups.Group, S = gens(G); halfradius::Integer)
S = union!(S, inv.(S)) S = union!(S, inv.(S))
@info "generating wl-metric ball of radius $(2halfradius)" @info "generating wl-metric ball of radius $(2halfradius)"
@time E, sizes = Groups.wlmetric_ball(S, radius=2halfradius) @time E, sizes = Groups.wlmetric_ball(S; radius = 2halfradius)
@info "sizes = $(sizes)" @info "sizes = $(sizes)"
@info "computing the *-algebra structure for G" @info "computing the *-algebra structure for G"
@time RG = StarAlgebras.StarAlgebra{twisted}( @time RG = StarAlgebras.StarAlgebra(
G, G,
StarAlgebras.Basis{UInt32}(E), StarAlgebras.Basis{UInt32}(E),
(sizes[halfradius], sizes[halfradius]), (sizes[halfradius], sizes[halfradius]);
precompute=false, precompute = false,
) )
return RG, S, sizes return RG, S, sizes
end end

3
src/actions/actions.jl
View File

@ -1,5 +1,4 @@
import SymbolicWedderburn.action import SymbolicWedderburn.action
StarAlgebras.star(g::Groups.GroupElement) = inv(g)
include("alphabet_permutation.jl") include("alphabet_permutation.jl")
@ -10,7 +9,7 @@ include("autfn_conjugation.jl")
function SymbolicWedderburn.action( function SymbolicWedderburn.action(
act::SymbolicWedderburn.ByPermutations, act::SymbolicWedderburn.ByPermutations,
g::Groups.GroupElement, g::Groups.GroupElement,
α::StarAlgebras.AlgebraElement α::StarAlgebras.AlgebraElement,
) )
res = StarAlgebras.zero!(similar(α)) res = StarAlgebras.zero!(similar(α))
B = basis(parent(α)) B = basis(parent(α))

78
src/certify.jl
View File

@ -1,3 +1,6 @@
import IntervalArithmetic
import IntervalMatrices
function augment_columns!(Q::AbstractMatrix) function augment_columns!(Q::AbstractMatrix)
for c in eachcol(Q) for c in eachcol(Q)
c .-= sum(c) ./ length(c) c .-= sum(c) ./ length(c)
@ -8,26 +11,25 @@ end
function __sos_via_sqr!( function __sos_via_sqr!(
res::StarAlgebras.AlgebraElement, res::StarAlgebras.AlgebraElement,
P::AbstractMatrix; P::AbstractMatrix;
augmented::Bool augmented::Bool,
id = (b = basis(parent(res)); b[one(first(b))]),
) )
StarAlgebras.zero!(res)
A = parent(res) A = parent(res)
b = basis(A) mstr = A.mstructure
@assert size(A.mstructure) == size(P) @assert size(mstr) == size(P)
e = b[one(b[1])]
for i in axes(A.mstructure, 1) StarAlgebras.zero!(res)
x = StarAlgebras._istwisted(A.mstructure) ? StarAlgebras.star(b[i]) : b[i] for j in axes(mstr, 2)
for j in axes(A.mstructure, 2) for i in axes(mstr, 1)
p = P[i, j] p = P[i, j]
xy = b[A.mstructure[i, j]] x_star_y = mstr[-i, j]
# either result += P[x,y]*(x*y) res[x_star_y] += p
res[xy] += p # either result += P[x,y]*(x'*y)
if augmented if augmented
# or result += P[x,y]*(1-x)*(1-y) == P[x,y]*(2-x-y+xy) # or result += P[x,y]*(1-x)'*(1-y) == P[x,y]*(1-x'-y+x'y)
y = b[j] res[id] += p
res[e] += p x_star, y = mstr[-i, id], j
res[x] -= p res[x_star] -= p
res[y] -= p res[y] -= p
end end
end end
@ -36,7 +38,11 @@ function __sos_via_sqr!(
return res return res
end end
function __sos_via_cnstr!(res::StarAlgebras.AlgebraElement, Q²::AbstractMatrix, cnstrs) function __sos_via_cnstr!(
res::StarAlgebras.AlgebraElement,
Q²::AbstractMatrix,
cnstrs,
)
StarAlgebras.zero!(res) StarAlgebras.zero!(res)
for (g, A_g) in cnstrs for (g, A_g) in cnstrs
res[g] = dot(A_g, Q²) res[g] = dot(A_g, Q²)
@ -44,10 +50,14 @@ function __sos_via_cnstr!(res::StarAlgebras.AlgebraElement, Q²::AbstractMatrix,
return res return res
end end
function compute_sos(A::StarAlgebras.StarAlgebra, Q::AbstractMatrix; augmented::Bool) function compute_sos(
A::StarAlgebras.StarAlgebra,
Q::AbstractMatrix;
augmented::Bool,
)
Q² = Q' * Q Q² = Q' * Q
res = StarAlgebras.AlgebraElement(zeros(eltype(Q²), length(basis(A))), A) res = StarAlgebras.AlgebraElement(zeros(eltype(Q²), length(basis(A))), A)
res = __sos_via_sqr!(res, Q², augmented=augmented) res = __sos_via_sqr!(res, Q²; augmented = augmented)
return res return res
end end
@ -56,7 +66,7 @@ function sufficient_λ(residual::StarAlgebras.AlgebraElement, λ; halfradius)
suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm
eq_sign = let T = eltype(residual) eq_sign = let T = eltype(residual)
if T <: Interval if T <: IntervalArithmetic.Interval
"∈" "∈"
elseif T <: Union{Rational,Integer} elseif T <: Union{Rational,Integer}
"=" "="
@ -81,13 +91,12 @@ function sufficient_λ(
order_unit::StarAlgebras.AlgebraElement, order_unit::StarAlgebras.AlgebraElement,
λ, λ,
sos::StarAlgebras.AlgebraElement; sos::StarAlgebras.AlgebraElement;
halfradius halfradius,
) )
@assert parent(elt) === parent(order_unit) == parent(sos) @assert parent(elt) === parent(order_unit) == parent(sos)
residual = (elt - λ * order_unit) - sos residual = (elt - λ * order_unit) - sos
return sufficient_λ(residual, λ; halfradius=halfradius) return sufficient_λ(residual, λ; halfradius = halfradius)
end end
function certify_solution( function certify_solution(
@ -96,24 +105,27 @@ function certify_solution(
λ, λ,
Q::AbstractMatrix{<:AbstractFloat}; Q::AbstractMatrix{<:AbstractFloat};
halfradius, halfradius,
augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)) augmented = iszero(StarAlgebras.aug(elt)) &&
iszero(StarAlgebras.aug(orderunit)),
) )
should_we_augment =
should_we_augment = !augmented && StarAlgebras.aug(elt) == StarAlgebras.aug(orderunit) == 0 !augmented && StarAlgebras.aug(elt) == StarAlgebras.aug(orderunit) == 0
Q = should_we_augment ? augment_columns!(Q) : Q Q = should_we_augment ? augment_columns!(Q) : Q
@time sos = compute_sos(parent(elt), Q, augmented=augmented) @time sos = compute_sos(parent(elt), Q; augmented = augmented)
@info "Checking in $(eltype(sos)) arithmetic with" λ @info "Checking in $(eltype(sos)) arithmetic with" λ
λ_flpoint = sufficient_λ(elt, orderunit, λ, sos, halfradius=halfradius) λ_flpoint = sufficient_λ(elt, orderunit, λ, sos; halfradius = halfradius)
if λ_flpoint ≤ 0 if λ_flpoint ≤ 0
return false, λ_flpoint return false, λ_flpoint
end end
λ_int = @interval(λ) λ_int = IntervalArithmetic.@interval(λ)
Q_int = [@interval(q) for q in Q] Q_int = IntervalMatrices.IntervalMatrix([
IntervalArithmetic.@interval(q) for q in Q
])
check, sos_int = @time if should_we_augment check, sos_int = @time if should_we_augment
@info("Projecting columns of Q to the augmentation ideal...") @info("Projecting columns of Q to the augmentation ideal...")
@ -123,16 +135,16 @@ function certify_solution(
check_augmented || @error( check_augmented || @error(
"Augmentation failed! The following numbers are not certified!" "Augmentation failed! The following numbers are not certified!"
) )
sos_int = compute_sos(parent(elt), Q_int; augmented=augmented) sos_int = compute_sos(parent(elt), Q_int; augmented = augmented)
check_augmented, sos_int check_augmented, sos_int
else else
true, compute_sos(parent(elt), Q_int, augmented=augmented) true, compute_sos(parent(elt), Q_int; augmented = augmented)
end end
@info "Checking in $(eltype(sos_int)) arithmetic with" λ_int @info "Checking in $(eltype(sos_int)) arithmetic with" λ_int
λ_certified = λ_certified =
sufficient_λ(elt, orderunit, λ_int, sos_int, halfradius=halfradius) sufficient_λ(elt, orderunit, λ_int, sos_int; halfradius = halfradius)
return check && inf(λ_certified) > 0.0, inf(λ_certified) return check && IntervalArithmetic.inf(λ_certified) > 0.0, λ_certified
end end

96
src/constraint_matrix.jl
View File

@ -28,7 +28,12 @@ struct ConstraintMatrix{T,I} <: AbstractMatrix{T}
size::Tuple{Int,Int} size::Tuple{Int,Int}
val::T val::T
function ConstraintMatrix{T}(nzeros::AbstractArray{<:Integer}, n, m, val) where {T} function ConstraintMatrix{T}(
nzeros::AbstractArray{<:Integer},
n,
m,
val,
) where {T}
@assert n ≥ 1 @assert n ≥ 1
@assert m ≥ 1 @assert m ≥ 1
@ -45,15 +50,26 @@ struct ConstraintMatrix{T,I} <: AbstractMatrix{T}
end end
end end
ConstraintMatrix(nzeros::AbstractArray{<:Integer}, n, m, val::T) where {T} = function ConstraintMatrix(
ConstraintMatrix{T}(nzeros, n, m, val) nzeros::AbstractArray{<:Integer},
n,
m,
val::T,
) where {T}
return ConstraintMatrix{T}(nzeros, n, m, val)
end
Base.size(cm::ConstraintMatrix) = cm.size Base.size(cm::ConstraintMatrix) = cm.size
__get_positive(cm::ConstraintMatrix, idx::Integer) = function __get_positive(cm::ConstraintMatrix, idx::Integer)
convert(eltype(cm), cm.val * length(searchsorted(cm.pos, idx))) return convert(eltype(cm), cm.val * length(searchsorted(cm.pos, idx)))
__get_negative(cm::ConstraintMatrix, idx::Integer) = end
convert(eltype(cm), cm.val * length(searchsorted(cm.neg, idx))) function __get_negative(cm::ConstraintMatrix, idx::Integer)
return convert(
eltype(cm),
cm.val * length(searchsorted(cm.neg, idx; rev = true)),
)
end
Base.@propagate_inbounds function Base.getindex( Base.@propagate_inbounds function Base.getindex(
cm::ConstraintMatrix, cm::ConstraintMatrix,
@ -67,26 +83,29 @@ Base.@propagate_inbounds function Base.getindex(
return pos - neg return pos - neg
end end
struct NZPairsIter{T} struct NZPairsIter{T,I}
m::ConstraintMatrix{T} m::ConstraintMatrix{T,I}
end end
Base.eltype(::Type{NZPairsIter{T}}) where {T} = Pair{Int,T} Base.eltype(::Type{NZPairsIter{T}}) where {T} = Pair{Int,T}
Base.IteratorSize(::Type{<:NZPairsIter}) = Base.SizeUnknown() Base.IteratorSize(::Type{<:NZPairsIter}) = Base.SizeUnknown()
# TODO: iterate over (idx=>val) pairs combining vals # TODO: iterate over (idx=>val) pairs combining vals
function Base.iterate(itr::NZPairsIter, state::Tuple{Int,Int}=(1, 1)) function Base.iterate(
itr::NZPairsIter,
state::Tuple{Int,Nothing} = (1, nothing),
)
k = iterate(itr.m.pos, state[1]) k = iterate(itr.m.pos, state[1])
isnothing(k) && return iterate(itr, state[2]) isnothing(k) && return iterate(itr, (nothing, 1))
idx, st = k idx, st = k
return idx => itr.m.val, (st, 1) return idx => itr.m.val, (st, nothing)
end end
function Base.iterate(itr::NZPairsIter, state::Int) function Base.iterate(itr::NZPairsIter, state::Tuple{Nothing,Int})
k = iterate(itr.m.neg, state[1]) k = iterate(itr.m.neg, state[2])
isnothing(k) && return nothing isnothing(k) && return nothing
idx, st = k idx, st = k
return idx => -itr.m.val, st return idx => -itr.m.val, (nothing, st)
end end
""" """
@ -127,3 +146,50 @@ function LinearAlgebra.dot(cm::ConstraintMatrix, m::AbstractMatrix{T}) where {T}
neg = isempty(cm.neg) ? zero(first(m)) : sum(@view m[cm.neg]) neg = isempty(cm.neg) ? zero(first(m)) : sum(@view m[cm.neg])
return convert(eltype(cm), cm.val) * (pos - neg) return convert(eltype(cm), cm.val) * (pos - neg)
end end
function constraints(A::StarAlgebras.StarAlgebra; augmented::Bool)
return constraints(basis(A), A.mstructure; augmented = augmented)
end
function constraints(
basis::StarAlgebras.AbstractBasis,
mstr::StarAlgebras.MultiplicativeStructure;
augmented = false,
)
cnstrs = _constraints(
mstr;
augmented = augmented,
num_constraints = length(basis),
id = basis[one(first(basis))],
)
return Dict(
basis[i] => ConstraintMatrix(c, size(mstr)..., 1) for
(i, c) in pairs(cnstrs)
)
end
function _constraints(
mstr::StarAlgebras.MultiplicativeStructure;
augmented::Bool = false,
num_constraints = maximum(mstr),
id,
)
cnstrs = [signed(eltype(mstr))[] for _ in 1:num_constraints]
LI = LinearIndices(size(mstr))
for ci in CartesianIndices(size(mstr))
k = LI[ci]
i, j = Tuple(ci)
a_star_b = mstr[-i, j]
push!(cnstrs[a_star_b], k)
if augmented
# (1-a)'(1-b) = 1 - a' - b + a'b
push!(cnstrs[id], k)
a_star, b = mstr[-i, id], j
push!(cnstrs[a_star], -k)
push!(cnstrs[b], -k)
end
end
return cnstrs
end

29
src/gradings.jl
View File

@ -1,7 +1,8 @@
## something about roots ## something about roots
Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} = function Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N}
Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j) return Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j)
end
function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N} function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N}
if s.symbol === :A if s.symbol === :A
@ -51,20 +52,28 @@ function Adj(rootsystem::AbstractDict, subtype::Symbol)
+, +,
( (
Δα * Δβ for (α, Δα) in rootsystem for (β, Δβ) in rootsystem if Δα * Δβ for (α, Δα) in rootsystem for (β, Δβ) in rootsystem if
Roots.classify_sub_root_system( Roots.classify_sub_root_system(roots, first(α), first(β)) == subtype
roots, );
first(α), init = zero(first(values(rootsystem))),
first(β), )
) == subtype end
),
init=zero(first(values(rootsystem))), Adj(rootsystem::AbstractDict) = sum(values(rootsystem))^2 - Sq(rootsystem)
function Sq(rootsystem::AbstractDict)
return reduce(
+,
Δα^2 for (_, Δα) in rootsystem;
init = zero(first(values(rootsystem))),
) )
end end
function level(rootsystem, level::Integer) function level(rootsystem, level::Integer)
1 ≤ level ≤ 4 || throw("level is implemented only for i ∈{1,2,3,4}") 1 ≤ level ≤ 4 || throw("level is implemented only for i ∈{1,2,3,4}")
level == 1 && return Adj(rootsystem, :C₁) # always positive level == 1 && return Adj(rootsystem, :C₁) # always positive
level == 2 && return Adj(rootsystem, :A₁) + Adj(rootsystem, Symbol("C₁×C₁")) + Adj(rootsystem, :C₂) # C₂ is not positive level == 2 && return Adj(rootsystem, :A₁) +
Adj(rootsystem, Symbol("C₁×C₁")) +
Adj(rootsystem, :C₂) # C₂ is not positive
level == 3 && return Adj(rootsystem, :A₂) + Adj(rootsystem, Symbol("A₁×C₁")) level == 3 && return Adj(rootsystem, :A₂) + Adj(rootsystem, Symbol("A₁×C₁"))
level == 4 && return Adj(rootsystem, Symbol("A₁×A₁")) # positive level == 4 && return Adj(rootsystem, Symbol("A₁×A₁")) # positive
end end

37
src/reconstruct.jl
View File

@ -7,33 +7,28 @@ end
function reconstruct( function reconstruct(
Ms::AbstractVector{<:AbstractMatrix}, Ms::AbstractVector{<:AbstractMatrix},
wbdec::WedderburnDecomposition; wbdec::WedderburnDecomposition,
atol=eps(real(eltype(wbdec))) * 10__outer_dim(wbdec)
) )
n = __outer_dim(wbdec) n = __outer_dim(wbdec)
res = sum(zip(Ms, SymbolicWedderburn.direct_summands(wbdec))) do (M, ds) res = sum(zip(Ms, SymbolicWedderburn.direct_summands(wbdec))) do (M, ds)
res = similar(M, n, n) res = similar(M, n, n)
reconstruct!(res, M, ds, __group_of(wbdec), wbdec.hom, atol=atol) res = _reconstruct!(res, M, ds)
return res
end end
res = average!(zero(res), res, __group_of(wbdec), wbdec.hom)
return res return res
end end
function reconstruct!( function _reconstruct!(
res::AbstractMatrix, res::AbstractMatrix,
M::AbstractMatrix, M::AbstractMatrix,
ds::SymbolicWedderburn.DirectSummand, ds::SymbolicWedderburn.DirectSummand,
G,
hom;
atol=eps(real(eltype(ds))) * 10max(size(res)...)
) )
res .= zero(eltype(res)) res .= zero(eltype(res))
U = SymbolicWedderburn.image_basis(ds) if !iszero(M)
d = SymbolicWedderburn.degree(ds) U = SymbolicWedderburn.image_basis(ds)
tmp = (U' * M * U) .* d d = SymbolicWedderburn.degree(ds)
res = (U' * M * U) .* d
res = average!(res, tmp, G, hom)
if eltype(res) <: AbstractFloat
__droptol!(res, atol) # TODO: is this really necessary?!
end end
return res return res
end end
@ -52,18 +47,22 @@ function average!(
res::AbstractMatrix, res::AbstractMatrix,
M::AbstractMatrix, M::AbstractMatrix,
G::Groups.Group, G::Groups.Group,
hom::SymbolicWedderburn.InducedActionHomomorphism{<:SymbolicWedderburn.ByPermutations} hom::SymbolicWedderburn.InducedActionHomomorphism{
<:SymbolicWedderburn.ByPermutations,
},
) )
res .= zero(eltype(res))
@assert size(M) == size(res) @assert size(M) == size(res)
o = Groups.order(Int, G)
for g in G for g in G
gext = SymbolicWedderburn.induce(hom, g) p = SymbolicWedderburn.induce(hom, g)
Threads.@threads for c in axes(res, 2) Threads.@threads for c in axes(res, 2)
for r in axes(res, 1) for r in axes(res, 1)
res[r, c] += M[r^gext, c^gext] if !iszero(M[r, c])
res[r^p, c^p] += M[r, c] / o
end
end end
end end
end end
o = Groups.order(Int, G)
res ./= o
return res return res
end end

216
src/roots.jl
View File

@ -5,111 +5,104 @@ using LinearAlgebra
export Root, isproportional, isorthogonal, ~, ⟂ export Root, isproportional, isorthogonal, ~, ⟂
abstract type AbstractRoot{N,T} end abstract type AbstractRoot{N,T} end # <: AbstractVector{T} ?
struct Root{N,T} <: AbstractRoot{N,T} ℓ₂length(r::AbstractRoot) = norm(r, 2)
coord::SVector{N,T} ambient_dim(r::AbstractRoot) = length(r)
end Base.:*(r::AbstractRoot, a::Number) = a * r
Root(a) = Root(SVector(a...))
function Base.:(==)(r::Root{N}, s::Root{M}) where {M,N}
M == N || return false
r.coord == s.coord || return false
return true
end
Base.hash(r::Root, h::UInt) = hash(r.coord, hash(Root, h))
Base.:+(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord + s.coord)
Base.:-(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord - s.coord)
Base.:-(r::Root{N}) where {N} = Root(-r.coord)
Base.:*(a::Number, r::Root) = Root(a * r.coord)
Base.:*(r::Root, a::Number) = a * r
Base.length(r::AbstractRoot) = norm(r, 2)
LinearAlgebra.norm(r::Root, p::Real=2) = norm(r.coord, p)
LinearAlgebra.dot(r::Root, s::Root) = dot(r.coord, s.coord)
cos_angle(a, b) = dot(a, b) / (norm(a) * norm(b)) cos_angle(a, b) = dot(a, b) / (norm(a) * norm(b))
function isproportional(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M} function isproportional(α::AbstractRoot, β::AbstractRoot)
N == M || return false ambient_dim(α) == ambient_dim(β) || return false
val = abs(cos_angle(α, β)) val = abs(cos_angle(α, β))
return isapprox(val, one(val), atol=eps(one(val))) return isapprox(val, one(val); atol = eps(one(val)))
end end
function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M} function isorthogonal(α::AbstractRoot, β::AbstractRoot)
N == M || return false ambient_dim(α) == ambient_dim(β) || return false
val = cos_angle(α, β) val = cos_angle(α, β)
return isapprox(val, zero(val), atol=eps(one(val))) return isapprox(val, zero(val); atol = eps(one(val)))
end end
function _positive_direction(α::Root{N}) where {N} function positive(roots::AbstractVector{<:AbstractRoot})
v = α.coord + 1 / (N * 100) * rand(N) isempty(roots) && return empty(roots)
return Root{N,Float64}(v / norm(v, 2))
end
function positive(roots::AbstractVector{<:Root{N}}) where {N}
# return those roots for which dot(α, Root([½, ¼, …])) > 0.0
pd = _positive_direction(first(roots)) pd = _positive_direction(first(roots))
return filter(α -> dot(α, pd) > 0.0, roots) return filter(α -> dot(α, pd) > 0.0, roots)
end end
function Base.show(io::IO, r::Root) function Base.show(io::IO, r::AbstractRoot)
print(io, "Root$(r.coord)") return print(io, "Root $(r.coord)")
end end
function Base.show(io::IO, ::MIME"text/plain", r::Root{N}) where {N} function Base.show(io::IO, ::MIME"text/plain", r::AbstractRoot)
lngth² = sum(x -> x^2, r.coord) l₂l = ℓ₂length(r)
l = isinteger(sqrt(lngth²)) ? "$(sqrt(lngth²))" : "√$(lngth²)" l = round(Int, l₂l) ≈ l₂l ? "$(round(Int, l₂l))" : "√$(round(Int, l₂l^2))"
print(io, "Root in ℝ^$N of length $l\n", r.coord) return print(io, "Root in ℝ^$(length(r)) of length $l\n", r.coord)
end end
𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N)) function reflection(α::AbstractRoot, β::AbstractRoot)
𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N)) return β - Int(2dot(α, β) // dot(α, α)) * α
end
function cartan(α::AbstractRoot, β::AbstractRoot)
ambient_dim(α) == ambient_dim(β) || throw("incompatible ambient dimensions")
return [
ℓ₂length(reflection(a, b) - b) / ℓ₂length(a) for a in (α, β),
b in (α, β)
]
end
""" """
classify_root_system(α, β) classify_root_system(α, β)
Return the symbol of smallest system generated by roots `α` and `β`. Return the symbol of smallest system generated by roots `α` and `β`.
The classification is based only on roots length and The classification is based only on roots length,
proportionality/orthogonality. proportionality/orthogonality and Cartan matrix.
""" """
function classify_root_system(α::AbstractRoot, β::AbstractRoot) function classify_root_system(
lα, lβ = length(α), length(β) α::AbstractRoot,
β::AbstractRoot,
long::Tuple{Bool,Bool},
)
if isproportional(α, β) if isproportional(α, β)
if lα ≈ lβ ≈ √2 if all(long)
return :A₁
elseif lα ≈ lβ ≈ 2.0
return :C₁ return :C₁
elseif all(.!long) # both short
return :A₁
else else
@error "Proportional roots of different length"
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β") error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
end end
elseif isorthogonal(α, β) elseif isorthogonal(α, β)
if lα ≈ lβ ≈ √2 if all(long)
return Symbol("A₁×A₁")
elseif lα ≈ lβ ≈ 2.0
return Symbol("C₁×C₁") return Symbol("C₁×C₁")
elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2) elseif all(.!long) # both short
return Symbol("A₁×A₁")
elseif any(long)
return Symbol("A₁×C₁") return Symbol("A₁×C₁")
else
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
end end
else # ⟨α, β⟩ is 2-dimensional, but they're not orthogonal else # ⟨α, β⟩ is 2-dimensional, but they're not orthogonal
if lα ≈ lβ ≈ √2 a, b, c, d = abs.(cartan(α, β))
@assert a == d == 2
b, c = b < c ? (b, c) : (c, b)
if b == c == 1
return :A₂ return :A₂
elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2) elseif b == 1 && c == 2
return :C₂ return :C₂
elseif b == 1 && c == 3
@warn ":G₂? really?"
return :G₂
else else
@error a, b, c, d
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β") error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
end end
end end
end end
function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root) function proportional_root_from_system(
Ω::AbstractVector{<:AbstractRoot},
α::AbstractRoot,
)
k = findfirst(v -> isproportional(α, v), Ω) k = findfirst(v -> isproportional(α, v), Ω)
if isnothing(k) if isnothing(k)
error("Line L_α not contained in root system Ω:\n α = $α\n Ω = $Ω") error("Line L_α not contained in root system Ω:\n α = $α\n Ω = $Ω")
@ -117,25 +110,31 @@ function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root)
return Ω[k] return Ω[k]
end end
struct Plane{R<:Root} struct Plane{R<:AbstractRoot}
v1::R v1::R
v2::R v2::R
vectors::Vector{R} vectors::Vector{R}
end end
Plane(α::R, β::R) where {R<:Root} = function Plane(α::AbstractRoot, β::AbstractRoot)
Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3]) return Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
end
function Base.in(r::R, plane::Plane{R}) where {R} function Base.in(r::AbstractRoot, plane::Plane)
return any(isproportional(r, v) for v in plane.vectors) return any(isproportional(r, v) for v in plane.vectors)
end end
function classify_sub_root_system( function _islong(α::AbstractRoot, Ω)
Ω::AbstractVector{<:Root{N}}, lα = ℓ₂length(α)
α::Root{N}, return any(r -> lα - ℓ₂length(r) > eps(lα), Ω)
β::Root{N}, end
) where {N}
function classify_sub_root_system(
Ω::AbstractVector{<:AbstractRoot{N}},
α::AbstractRoot{N},
β::AbstractRoot{N},
) where {N}
@assert 1 ≤ length(unique(ℓ₂length, Ω)) ≤ 2
v = proportional_root_from_system(Ω, α) v = proportional_root_from_system(Ω, α)
w = proportional_root_from_system(Ω, β) w = proportional_root_from_system(Ω, β)
@ -146,32 +145,75 @@ function classify_sub_root_system(
l = length(subsystem) l = length(subsystem)
if l == 1 if l == 1
x = first(subsystem) x = first(subsystem)
return classify_root_system(x, x) long = _islong(x, Ω)
return classify_root_system(x, -x, (long, long))
elseif l == 2 elseif l == 2
return classify_root_system(subsystem...) x, y = subsystem
return classify_root_system(x, y, (_islong(x, Ω), _islong(y, Ω)))
elseif l == 3 elseif l == 3
a = classify_root_system(subsystem[1], subsystem[2]) x, y, z = subsystem
b = classify_root_system(subsystem[2], subsystem[3]) l1, l2, l3 = _islong(x, Ω), _islong(y, Ω), _islong(z, Ω)
c = classify_root_system(subsystem[1], subsystem[3]) a = classify_root_system(x, y, (l1, l2))
b = classify_root_system(y, z, (l2, l3))
c = classify_root_system(x, z, (l1, l3))
if a == b == c # it's only A₂ if :A₂ == a == b == c # it's only A₂
return a return a
end end
C = (:C₂, Symbol("C₁×C₁")) throw("Unknown subroot system! $((x,y,z))")
if (a ∈ C && b ∈ C && c ∈ C) && (:C₂ ∈ (a, b, c))
return :C₂
end
elseif l == 4 elseif l == 4
for i = 1:l subtypes = [
for j = (i+1):l classify_root_system(x, y, (_islong(x, Ω), _islong(y, Ω))) for
T = classify_root_system(subsystem[i], subsystem[j]) x in subsystem for y in subsystem if x ≠ y
T == :C₂ && return :C₂ ]
end if :C₂ in subtypes
return :C₂
end end
end end
@error "Unknown root subsystem generated by" α β @error "Unknown root subsystem generated by" α β
throw("Unknown root system: $subsystem") throw("Unknown root system: $subsystem")
end end
## concrete implementation:
struct Root{N,T} <: AbstractRoot{N,T}
coord::SVector{N,T}
end
Root(a) = Root(SVector(a...))
# convienience constructors
𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
function Base.:(==)(r::Root{N}, s::Root{M}) where {M,N}
M == N || return false
r.coord == s.coord || return false
return true
end
Base.hash(r::Root, h::UInt) = hash(r.coord, hash(Root, h))
function Base.:+(r::Root, s::Root)
ambient_dim(r) == ambient_dim(s) || throw("incompatible ambient dimensions")
return Root(r.coord + s.coord)
end
function Base.:-(r::Root, s::Root)
ambient_dim(r) == ambient_dim(s) || throw("incompatible ambient dimensions")
return Root(r.coord - s.coord)
end
Base.:-(r::Root) = Root(-r.coord)
Base.:*(a::Number, r::Root) = Root(a * r.coord)
Base.length(r::Root) = length(r.coord)
LinearAlgebra.norm(r::Root, p::Real = 2) = norm(r.coord, p)
LinearAlgebra.dot(r::Root, s::Root) = dot(r.coord, s.coord)
function _positive_direction(α::Root{N}) where {N}
v = α.coord + 1 / (N * 100) * rand(N)
return Root{N,Float64}(v / norm(v, 2))
end
end # of module Roots end # of module Roots

23
src/solve.jl
View File

@ -7,12 +7,15 @@ function setwarmstart!(model::JuMP.Model, warmstart)
ct => JuMP.all_constraints(model, ct...) for ct => JuMP.all_constraints(model, ct...) for
ct in JuMP.list_of_constraint_types(model) ct in JuMP.list_of_constraint_types(model)
) )
try
JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal)
JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal) for (ct, idx) in pairs(constraint_map)
JuMP.set_start_value.(idx, warmstart.slack[ct])
for (ct, idx) in pairs(constraint_map) JuMP.set_dual_start_value.(idx, warmstart.dual[ct])
JuMP.set_start_value.(idx, warmstart.slack[ct]) end
JuMP.set_dual_start_value.(idx, warmstart.dual[ct]) catch e
@warn "Warmstarting was not succesfull!" e
end end
return model return model
end end
@ -28,11 +31,10 @@ function getwarmstart(model::JuMP.Model)
slack = Dict(k => value.(v) for (k, v) in constraint_map) slack = Dict(k => value.(v) for (k, v) in constraint_map)
duals = Dict(k => JuMP.dual.(v) for (k, v) in constraint_map) duals = Dict(k => JuMP.dual.(v) for (k, v) in constraint_map)
return (primal=primal, dual=duals, slack=slack) return (primal = primal, dual = duals, slack = slack)
end end
function solve(m::JuMP.Model, optimizer, warmstart=nothing) function solve(m::JuMP.Model, optimizer, warmstart = nothing)
JuMP.set_optimizer(m, optimizer) JuMP.set_optimizer(m, optimizer)
MOIU.attach_optimizer(m) MOIU.attach_optimizer(m)
@ -45,8 +47,7 @@ function solve(m::JuMP.Model, optimizer, warmstart=nothing)
return status, getwarmstart(m) return status, getwarmstart(m)
end end
function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart=nothing) function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart = nothing)
isdir(dirname(solverlog)) || mkpath(dirname(solverlog)) isdir(dirname(solverlog)) || mkpath(dirname(solverlog))
Base.flush(Base.stdout) Base.flush(Base.stdout)
@ -55,7 +56,7 @@ function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart=nothing)
redirect_stdout(logfile) do redirect_stdout(logfile) do
status, warmstart = solve(m, optimizer, warmstart) status, warmstart = solve(m, optimizer, warmstart)
Base.Libc.flush_cstdio() Base.Libc.flush_cstdio()
status, warmstart return status, warmstart
end end
end end

205
src/sos_sdps.jl
View File

@ -6,15 +6,13 @@ See also [sos_problem_primal](@ref).
""" """
function sos_problem_dual( function sos_problem_dual(
elt::StarAlgebras.AlgebraElement, elt::StarAlgebras.AlgebraElement,
order_unit::StarAlgebras.AlgebraElement=zero(elt); order_unit::StarAlgebras.AlgebraElement = zero(elt);
lower_bound=-Inf lower_bound = -Inf,
) )
@assert parent(elt) == parent(order_unit) @assert parent(elt) == parent(order_unit)
algebra = parent(elt) algebra = parent(elt)
mstructure = if StarAlgebras._istwisted(algebra.mstructure) moment_matrix = let m = algebra.mstructure
algebra.mstructure [m[-i, j] for i in axes(m, 1), j in axes(m, 2)]
else
StarAlgebras.MTable{true}(basis(algebra), table_size=size(algebra.mstructure))
end end
# 1 variable for every primal constraint # 1 variable for every primal constraint
@ -22,60 +20,21 @@ function sos_problem_dual(
# Symmetrized: # Symmetrized:
# 1 dual variable for every orbit of G acting on basis # 1 dual variable for every orbit of G acting on basis
model = Model() model = Model()
@variable model y[1:length(basis(algebra))] JuMP.@variable(model, y[1:length(basis(algebra))])
@constraint model λ_dual dot(order_unit, y) == 1 JuMP.@constraint(model, λ_dual, dot(order_unit, y) == 1)
@constraint(model, psd, y[mstructure] in PSDCone()) JuMP.@constraint(model, psd, y[moment_matrix] in PSDCone())
if !isinf(lower_bound) if !isinf(lower_bound)
throw("Not Implemented yet") throw("Not Implemented yet")
@variable model λ_ub_dual JuMP.@variable(model, λ_ub_dual)
@objective model Min dot(elt, y) + lower_bound * λ_ub_dual JuMP.@objective(model, Min, dot(elt, y) + lower_bound * λ_ub_dual)
else else
@objective model Min dot(elt, y) JuMP.@objective(model, Min, dot(elt, y))
end end
return model return model
end end
function constraints(
basis::StarAlgebras.AbstractBasis,
mstr::AbstractMatrix{<:Integer};
augmented::Bool=false,
table_size=size(mstr)
)
cnstrs = [signed(eltype(mstr))[] for _ in basis]
LI = LinearIndices(table_size)
for ci in CartesianIndices(table_size)
k = LI[ci]
a_star_b = basis[mstr[k]]
push!(cnstrs[basis[a_star_b]], k)
if augmented
# (1-a_star)(1-b) = 1 - a_star - b + a_star_b
i, j = Tuple(ci)
a, b = basis[i], basis[j]
push!(cnstrs[basis[one(a)]], k)
push!(cnstrs[basis[StarAlgebras.star(a)]], -k)
push!(cnstrs[basis[b]], -k)
end
end
return Dict(
basis[i] => ConstraintMatrix(c, table_size..., 1) for (i, c) in pairs(cnstrs)
)
end
function constraints(A::StarAlgebras.StarAlgebra; augmented::Bool, twisted::Bool)
mstructure = if StarAlgebras._istwisted(A.mstructure) == twisted
A.mstructure
else
StarAlgebras.MTable{twisted}(basis(A), table_size=size(A.mstructure))
end
return constraints(basis(A), mstructure, augmented=augmented)
end
""" """
sos_problem_primal(X, [u = zero(X); upper_bound=Inf]) sos_problem_primal(X, [u = zero(X); upper_bound=Inf])
Formulate sum of squares decomposition problem for `X - λ·u`. Formulate sum of squares decomposition problem for `X - λ·u`.
@ -96,9 +55,10 @@ The default `u = zero(X)` formulates a simple feasibility problem.
""" """
function sos_problem_primal( function sos_problem_primal(
elt::StarAlgebras.AlgebraElement, elt::StarAlgebras.AlgebraElement,
order_unit::StarAlgebras.AlgebraElement=zero(elt); order_unit::StarAlgebras.AlgebraElement = zero(elt);
upper_bound=Inf, upper_bound = Inf,
augmented::Bool=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(order_unit)) augmented::Bool = iszero(StarAlgebras.aug(elt)) &&
iszero(StarAlgebras.aug(order_unit)),
) )
@assert parent(elt) === parent(order_unit) @assert parent(elt) === parent(order_unit)
@ -111,7 +71,7 @@ function sos_problem_primal(
@warn "Setting `upper_bound` together with zero `order_unit` has no effect" @warn "Setting `upper_bound` together with zero `order_unit` has no effect"
end end
A = constraints(parent(elt), augmented=augmented, twisted=true) A = constraints(parent(elt); augmented = augmented)
if !iszero(order_unit) if !iszero(order_unit)
λ = JuMP.@variable(model, λ) λ = JuMP.@variable(model, λ)
@ -135,11 +95,11 @@ end
function invariant_constraint!( function invariant_constraint!(
result::AbstractMatrix, result::AbstractMatrix,
basis::StarAlgebras.AbstractBasis, basis::StarAlgebras.AbstractBasis,
cnstrs::AbstractDict{K,CM}, cnstrs::AbstractDict{K,<:ConstraintMatrix},
invariant_vec::SparseVector, invariant_vec::SparseVector,
) where {K,CM<:ConstraintMatrix} ) where {K}
result .= zero(eltype(result)) result .= zero(eltype(result))
for i in SparseArrays.nonzeroinds(invariant_vec) @inbounds for i in SparseArrays.nonzeroinds(invariant_vec)
g = basis[i] g = basis[i]
A = cnstrs[g] A = cnstrs[g]
for (idx, v) in nzpairs(A) for (idx, v) in nzpairs(A)
@ -149,25 +109,47 @@ function invariant_constraint!(
return result return result
end end
function invariant_constraint(basis, cnstrs, invariant_vec)
I = UInt32[]
J = UInt32[]
V = Float64[]
_M = first(values(cnstrs))
CI = CartesianIndices(_M)
@inbounds for i in SparseArrays.nonzeroinds(invariant_vec)
g = basis[i]
A = cnstrs[g]
for (idx, v) in nzpairs(A)
ci = CI[idx]
push!(I, ci[1])
push!(J, ci[2])
push!(V, invariant_vec[i] * v)
end
end
return sparse(I, J, V, size(_M)...)
end
function isorth_projection(ds::SymbolicWedderburn.DirectSummand) function isorth_projection(ds::SymbolicWedderburn.DirectSummand)
U = SymbolicWedderburn.image_basis(ds) U = SymbolicWedderburn.image_basis(ds)
return isapprox(U * U', I) return isapprox(U * U', I)
end end
sos_problem_primal( function sos_problem_primal(
elt::StarAlgebras.AlgebraElement, elt::StarAlgebras.AlgebraElement,
wedderburn::WedderburnDecomposition; wedderburn::WedderburnDecomposition;
kwargs... kwargs...,
) = sos_problem_primal(elt, zero(elt), wedderburn; kwargs...) )
return sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)
end
function __fast_recursive_dot!( function __fast_recursive_dot!(
res::JuMP.AffExpr, res::JuMP.AffExpr,
Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}}, Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}},
Ms::AbstractArray{<:AbstractSparseMatrix}; Ms::AbstractArray{<:AbstractSparseMatrix},
weights,
) )
@assert length(Ps) == length(Ms) @assert length(Ps) == length(Ms)
for (A, P) in zip(Ms, Ps) for (w, A, P) in zip(weights, Ms, Ps)
iszero(Ms) && continue iszero(Ms) && continue
rows = rowvals(A) rows = rowvals(A)
vals = nonzeros(A) vals = nonzeros(A)
@ -175,13 +157,21 @@ function __fast_recursive_dot!(
for i in nzrange(A, cidx) for i in nzrange(A, cidx)
ridx = rows[i] ridx = rows[i]
val = vals[i] val = vals[i]
JuMP.add_to_expression!(res, P[ridx, cidx], val) JuMP.add_to_expression!(res, P[ridx, cidx], w * val)
end end
end end
end end
return res return res
end end
function _dot(
Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}},
Ms::AbstractArray{<:AbstractMatrix{T}},
weights = Iterators.repeated(one(T), length(Ms)),
) where {T}
return __fast_recursive_dot!(JuMP.AffExpr(), Ps, Ms, weights)
end
import ProgressMeter import ProgressMeter
__show_itrs(n, total) = () -> [(Symbol("constraint"), "$n/$total")] __show_itrs(n, total) = () -> [(Symbol("constraint"), "$n/$total")]
@ -189,21 +179,38 @@ function sos_problem_primal(
elt::StarAlgebras.AlgebraElement, elt::StarAlgebras.AlgebraElement,
orderunit::StarAlgebras.AlgebraElement, orderunit::StarAlgebras.AlgebraElement,
wedderburn::WedderburnDecomposition; wedderburn::WedderburnDecomposition;
upper_bound=Inf, upper_bound = Inf,
augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)), augmented = iszero(StarAlgebras.aug(elt)) &&
check_orthogonality=true, iszero(StarAlgebras.aug(orderunit)),
show_progress=false check_orthogonality = true,
show_progress = false,
) )
@assert parent(elt) === parent(orderunit) @assert parent(elt) === parent(orderunit)
if check_orthogonality if check_orthogonality
if any(!isorth_projection, direct_summands(wedderburn)) if any(!isorth_projection, direct_summands(wedderburn))
error("Wedderburn decomposition contains a non-orthogonal projection") error(
"Wedderburn decomposition contains a non-orthogonal projection",
)
end end
end end
id_one = findfirst(invariant_vectors(wedderburn)) do v
b = basis(parent(elt))
return sparsevec([b[one(first(b))]], [1 // 1], length(v)) == v
end
prog = ProgressMeter.Progress(
length(invariant_vectors(wedderburn));
dt = 1,
desc = "Adding constraints: ",
enabled = show_progress,
barlen = 60,
showspeed = true,
)
feasibility_problem = iszero(orderunit) feasibility_problem = iszero(orderunit)
# problem creation starts here
model = JuMP.Model() model = JuMP.Model()
if !feasibility_problem # add λ or not? if !feasibility_problem # add λ or not?
λ = JuMP.@variable(model, λ) λ = JuMP.@variable(model, λ)
@ -217,58 +224,52 @@ function sos_problem_primal(
end end
end end
P = map(direct_summands(wedderburn)) do ds # semidefinite constraints as described by wedderburn
Ps = map(direct_summands(wedderburn)) do ds
dim = size(ds, 1) dim = size(ds, 1)
P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric) P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric)
JuMP.@constraint(model, P in PSDCone()) JuMP.@constraint(model, P in PSDCone())
P return P
end end
begin # preallocating begin # Ms are preallocated for the constraints loop
T = eltype(wedderburn) T = eltype(wedderburn)
Ms = [spzeros.(T, size(p)...) for p in P] Ms = [spzeros.(T, size(p)...) for p in Ps]
M_orb = zeros(T, size(parent(elt).mstructure)...) _eps = 10 * eps(T) * max(size(parent(elt).mstructure)...)
end end
X = convert(Vector{T}, StarAlgebras.coeffs(elt)) X = StarAlgebras.coeffs(elt)
U = convert(Vector{T}, StarAlgebras.coeffs(orderunit)) U = StarAlgebras.coeffs(orderunit)
# defining constraints based on the multiplicative structure # defining constraints based on the multiplicative structure
cnstrs = constraints(parent(elt), augmented=augmented, twisted=true) cnstrs = constraints(parent(elt); augmented = augmented)
prog = ProgressMeter.Progress(
length(invariant_vectors(wedderburn)),
dt=1,
desc="Adding constraints... ",
enabled=show_progress,
barlen=60,
showspeed=true
)
# adding linear constraints: one per orbit
for (i, iv) in enumerate(invariant_vectors(wedderburn)) for (i, iv) in enumerate(invariant_vectors(wedderburn))
ProgressMeter.next!(prog, showvalues=__show_itrs(i, prog.n)) ProgressMeter.next!(prog; showvalues = __show_itrs(i, prog.n))
augmented && i == id_one && continue
# i == 500 && break
x = dot(X, iv) x = dot(X, iv)
u = dot(U, iv) u = dot(U, iv)
M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv) spM_orb = invariant_constraint(basis(parent(elt)), cnstrs, iv)
Ms = SymbolicWedderburn.diagonalize!(Ms, M_orb, wedderburn)
SparseArrays.droptol!.(Ms, 10 * eps(T) * max(size(M_orb)...))
# @info [nnz(m) / length(m) for m in Ms]
Ms = SymbolicWedderburn.diagonalize!(
Ms,
spM_orb,
wedderburn;
trace_preserving = true,
)
for M in Ms
SparseArrays.droptol!(M, _eps)
end
if feasibility_problem if feasibility_problem
JuMP.@constraint( JuMP.@constraint(model, x == _dot(Ps, Ms))
model,
x == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
)
else else
JuMP.@constraint( JuMP.@constraint(model, x - λ * u == _dot(Ps, Ms))
model,
x - λ * u == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
)
end end
end end
ProgressMeter.finish!(prog) ProgressMeter.finish!(prog)
return model, P return model, Ps
end end

155
test/1703.09680.jl
View File

@ -1,9 +1,8 @@
@testset "1703.09680 Examples" begin @testset "1703.09680 Examples" begin
@testset "SL(2,Z)" begin @testset "SL(2,Z)" begin
N = 2 N = 2
G = MatrixGroups.SpecialLinearGroup{N}(Int8) G = MatrixGroups.SpecialLinearGroup{N}(Int8)
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -15,15 +14,15 @@
status, certified, λ = check_positivity( status, certified, λ = check_positivity(
elt, elt,
unit, unit;
upper_bound=ub, upper_bound = ub,
halfradius=2, halfradius = 2,
optimizer=scs_optimizer( optimizer = scs_optimizer(;
eps=1e-10, eps = 1e-10,
max_iters=5_000, max_iters = 5_000,
accel=50, accel = 50,
alpha=1.9, alpha = 1.9,
) ),
) )
@test status == JuMP.ALMOST_OPTIMAL @test status == JuMP.ALMOST_OPTIMAL
@ -33,8 +32,10 @@
@testset "SL(3,F₅)" begin @testset "SL(3,F₅)" begin
N = 3 N = 3
G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{5}) G = MatrixGroups.SpecialLinearGroup{N}(
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) SymbolicWedderburn.Characters.FiniteFields.GF{5},
)
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -46,15 +47,15 @@
status, certified, λ = check_positivity( status, certified, λ = check_positivity(
elt, elt,
unit, unit;
upper_bound=ub, upper_bound = ub,
halfradius=2, halfradius = 2,
optimizer=scs_optimizer( optimizer = scs_optimizer(;
eps=1e-10, eps = 1e-10,
max_iters=5_000, max_iters = 5_000,
accel=50, accel = 50,
alpha=1.9, alpha = 1.9,
) ),
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@ -62,12 +63,15 @@
@test λ > 1 @test λ > 1
m = PropertyT.sos_problem_dual(elt, unit) m = PropertyT.sos_problem_dual(elt, unit)
PropertyT.solve(m, cosmo_optimizer( PropertyT.solve(
eps=1e-6, m,
max_iters=5_000, cosmo_optimizer(;
accel=50, eps = 1e-6,
alpha=1.9, max_iters = 5_000,
)) accel = 50,
alpha = 1.9,
),
)
@test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL) @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
@test JuMP.objective_value(m) ≈ 1.5 atol = 1e-2 @test JuMP.objective_value(m) ≈ 1.5 atol = 1e-2
@ -76,7 +80,7 @@
@testset "SAut(F₂)" begin @testset "SAut(F₂)" begin
N = 2 N = 2
G = SpecialAutomorphismGroup(FreeGroup(N)) G = SpecialAutomorphismGroup(FreeGroup(N))
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -88,53 +92,46 @@
status, certified, λ = check_positivity( status, certified, λ = check_positivity(
elt, elt,
unit, unit;
upper_bound=ub, upper_bound = ub,
halfradius=2, halfradius = 2,
optimizer=scs_optimizer( optimizer = scs_optimizer(;
eps=1e-10, eps = 1e-10,
max_iters=5_000, max_iters = 5_000,
accel=50, accel = 50,
alpha=1.9, alpha = 1.9,
) ),
) )
@test status == JuMP.ALMOST_OPTIMAL @test status == JuMP.ALMOST_OPTIMAL
@test λ < 0 @test λ < 0
@test !certified @test !certified
@time sos_problem = @time sos_problem = PropertyT.sos_problem_primal(elt; upper_bound = ub)
PropertyT.sos_problem_primal(elt, upper_bound=ub)
status, _ = PropertyT.solve( status, _ = PropertyT.solve(
sos_problem, sos_problem,
cosmo_optimizer( cosmo_optimizer(;
eps=1e-7, eps = 1e-7,
max_iters=10_000, max_iters = 10_000,
accel=0, accel = 0,
alpha=1.9, alpha = 1.9,
) ),
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
P = JuMP.value.(sos_problem[:P]) P = JuMP.value.(sos_problem[:P])
Q = real.(sqrt(P)) Q = real.(sqrt(P))
certified, λ_cert = PropertyT.certify_solution( certified, λ_cert =
elt, PropertyT.certify_solution(elt, zero(elt), 0.0, Q; halfradius = 2)
zero(elt),
0.0,
Q,
halfradius=2,
)
@test !certified @test !certified
@test λ_cert < 0 @test λ_cert < 0
end end
@testset "SL(3,Z) has (T)" begin @testset "SL(3,Z) has (T)" begin
n = 3 n = 3
SL = MatrixGroups.SpecialLinearGroup{n}(Int8) SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true) RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2)
Δ = RSL(length(S)) - sum(RSL(s) for s in S) Δ = RSL(length(S)) - sum(RSL(s) for s in S)
@ -145,18 +142,18 @@
opt_problem = PropertyT.sos_problem_primal( opt_problem = PropertyT.sos_problem_primal(
elt, elt,
unit, unit;
upper_bound=ub, upper_bound = ub,
augmented=false, augmented = false,
) )
status, _ = PropertyT.solve( status, _ = PropertyT.solve(
opt_problem, opt_problem,
cosmo_optimizer( cosmo_optimizer(;
eps=1e-10, eps = 1e-10,
max_iters=10_000, max_iters = 10_000,
accel=0, accel = 0,
alpha=1.5, alpha = 1.5,
), ),
) )
@ -170,13 +167,13 @@
elt, elt,
unit, unit,
λ, λ,
Q, Q;
halfradius=2, halfradius = 2,
augmented=false, augmented = false,
) )
@test certified @test certified
@test isapprox(λ_cert, λ, rtol=1e-5) @test isapprox(λ_cert, λ, rtol = 1e-5)
end end
@testset "augmented formulation" begin @testset "augmented formulation" begin
@ -186,18 +183,18 @@
opt_problem = PropertyT.sos_problem_primal( opt_problem = PropertyT.sos_problem_primal(
elt, elt,
unit, unit;
upper_bound=ub, upper_bound = ub,
augmented=true, augmented = true,
) )
status, _ = PropertyT.solve( status, _ = PropertyT.solve(
opt_problem, opt_problem,
scs_optimizer( scs_optimizer(;
eps=1e-10, eps = 1e-10,
max_iters=10_000, max_iters = 10_000,
accel=-10, accel = -10,
alpha=1.5, alpha = 1.5,
), ),
) )
@ -210,13 +207,13 @@
elt, elt,
unit, unit,
λ, λ,
Q, Q;
halfradius=2, halfradius = 2,
augmented=true, augmented = true,
) )
@test certified @test certified
@test isapprox(λ_cert, λ, rtol=1e-5) @test isapprox(λ_cert, λ, rtol = 1e-5)
@test λ_cert > 2 // 10 @test λ_cert > 2 // 10
end end
end end

88
test/1712.07167.jl
View File

@ -1,9 +1,9 @@
@testset "1712.07167 Examples" begin @testset "1712.07167 Examples" begin
@testset "SAut(F₃)" begin @testset "SAut(F₃)" begin
N = 3 N = 3
G = SpecialAutomorphismGroup(FreeGroup(N)) G = SpecialAutomorphismGroup(FreeGroup(N))
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) @info "running tests for" G
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1))) P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1)))
Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
@ -15,6 +15,7 @@
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
) )
@info wd
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -27,14 +28,14 @@
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
unit, unit,
wd, wd;
upper_bound=ub, upper_bound = ub,
halfradius=2, halfradius = 2,
optimizer=cosmo_optimizer( optimizer = cosmo_optimizer(;
eps=1e-7, eps = 1e-7,
max_iters=10_000, max_iters = 10_000,
accel=50, accel = 50,
alpha=1.9, alpha = 1.9,
), ),
) )
@ -47,7 +48,8 @@
n = 3 n = 3
SL = MatrixGroups.SpecialLinearGroup{n}(Int8) SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true) @info "running tests for" SL
RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2)
Δ = RSL(length(S)) - sum(RSL(s) for s in S) Δ = RSL(length(S)) - sum(RSL(s) for s in S)
@ -62,6 +64,7 @@
basis(RSL), basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
) )
@info wd
elt = Δ^2 elt = Δ^2
unit = Δ unit = Δ
@ -71,8 +74,8 @@
elt, elt,
unit, unit,
wd, wd,
upper_bound=ub, upper_bound = ub,
augmented=false, augmented = false,
) )
wdfl = SymbolicWedderburn.WedderburnDecomposition( wdfl = SymbolicWedderburn.WedderburnDecomposition(
@ -86,18 +89,18 @@
model, varP = PropertyT.sos_problem_primal( model, varP = PropertyT.sos_problem_primal(
elt, elt,
unit, unit,
wdfl, wdfl;
upper_bound=ub, upper_bound = ub,
augmented=false, augmented = false,
) )
status, warm = PropertyT.solve( status, warm = PropertyT.solve(
model, model,
cosmo_optimizer( cosmo_optimizer(;
eps=1e-10, eps = 1e-10,
max_iters=20_000, max_iters = 20_000,
accel=50, accel = 50,
alpha=1.9, alpha = 1.9,
), ),
) )
@ -105,34 +108,34 @@
status, _ = PropertyT.solve( status, _ = PropertyT.solve(
model, model,
scs_optimizer( scs_optimizer(;
eps=1e-10, eps = 1e-10,
max_iters=100, max_iters = 100,
accel=-20, accel = -20,
alpha=1.2, alpha = 1.2,
), ),
warm warm,
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
Q = @time let varP = varP Q = @time let varP = varP
Qs = map(varP) do P Qs = map(varP) do P
real.(sqrt(JuMP.value.(P))) return real.(sqrt(JuMP.value.(P)))
end end
PropertyT.reconstruct(Qs, wdfl) PropertyT.reconstruct(Qs, wdfl)
end end
λ = JuMP.value(model[:λ]) λ = JuMP.value(model[:λ])
sos = PropertyT.compute_sos(parent(elt), Q; augmented=false) sos = PropertyT.compute_sos(parent(elt), Q; augmented = false)
certified, λ_cert = PropertyT.certify_solution( certified, λ_cert = PropertyT.certify_solution(
elt, elt,
unit, unit,
λ, λ,
Q, Q;
halfradius=2, halfradius = 2,
augmented=false, augmented = false,
) )
@test certified @test certified
@ -155,22 +158,23 @@
basis(RSL), basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
) )
@info wdfl
opt_problem, varP = PropertyT.sos_problem_primal( opt_problem, varP = PropertyT.sos_problem_primal(
elt, elt,
unit, unit,
wdfl, wdfl;
upper_bound=ub, upper_bound = ub,
# augmented = true # since both elt and unit are augmented # augmented = true # since both elt and unit are augmented
) )
status, _ = PropertyT.solve( status, _ = PropertyT.solve(
opt_problem, opt_problem,
scs_optimizer( scs_optimizer(;
eps=1e-8, eps = 1e-8,
max_iters=20_000, max_iters = 20_000,
accel=0, accel = 0,
alpha=1.9, alpha = 1.9,
), ),
) )
@ -178,7 +182,7 @@
Q = @time let varP = varP Q = @time let varP = varP
Qs = map(varP) do P Qs = map(varP) do P
real.(sqrt(JuMP.value.(P))) return real.(sqrt(JuMP.value.(P)))
end end
PropertyT.reconstruct(Qs, wdfl) PropertyT.reconstruct(Qs, wdfl)
end end
@ -187,8 +191,8 @@
elt, elt,
unit, unit,
JuMP.objective_value(opt_problem), JuMP.objective_value(opt_problem),
Q, Q;
halfradius=2, halfradius = 2,
# augmented = true # since both elt and unit are augmented # augmented = true # since both elt and unit are augmented
) )

67
test/1812.03456.jl
View File

@ -21,8 +21,9 @@ using SparseArrays
@testset "Sq, Adj, Op in SL(4,Z)" begin @testset "Sq, Adj, Op in SL(4,Z)" begin
N = 4 N = 4
G = MatrixGroups.SpecialLinearGroup{N}(Int8) G = MatrixGroups.SpecialLinearGroup{N}(Int8)
@info "running tests for" G
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -39,6 +40,7 @@ using SparseArrays
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
) )
@info wd
ivs = SymbolicWedderburn.invariant_vectors(wd) ivs = SymbolicWedderburn.invariant_vectors(wd)
sq, adj, op = PropertyT.SqAdjOp(RG, N) sq, adj, op = PropertyT.SqAdjOp(RG, N)
@ -82,7 +84,7 @@ using SparseArrays
@testset "SAut(F₃)" begin @testset "SAut(F₃)" begin
n = 3 n = 3
G = SpecialAutomorphismGroup(FreeGroup(n)) G = SpecialAutomorphismGroup(FreeGroup(n))
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
sq, adj, op = PropertyT.SqAdjOp(RG, n) sq, adj, op = PropertyT.SqAdjOp(RG, n)
@test sq(one(G)) == 216 @test sq(one(G)) == 216
@ -98,7 +100,8 @@ end
n = 3 n = 3
G = MatrixGroups.SpecialLinearGroup{n}(Int8) G = MatrixGroups.SpecialLinearGroup{n}(Int8)
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) @info "running tests for" G
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -113,6 +116,7 @@ end
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
) )
@info wd
sq, adj, op = PropertyT.SqAdjOp(RG, n) sq, adj, op = PropertyT.SqAdjOp(RG, n)
@ -123,10 +127,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd, wd;
upper_bound=UB, upper_bound = UB,
halfradius=2, halfradius = 2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9) optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test certified @test certified
@ -140,10 +144,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd, wd;
upper_bound=UB, upper_bound = UB,
halfradius=2, halfradius = 2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9) optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test certified @test certified
@ -152,10 +156,11 @@ end
m, _ = PropertyT.sos_problem_primal(elt, wd) m, _ = PropertyT.sos_problem_primal(elt, wd)
PropertyT.solve( PropertyT.solve(
m, m,
scs_optimizer(max_iters=5000, accel=50, alpha=1.9) scs_optimizer(; max_iters = 5000, accel = 50, alpha = 1.9),
) )
@test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL, JuMP.ITERATION_LIMIT) @test JuMP.termination_status(m) in
(JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL, JuMP.ITERATION_LIMIT)
@test abs(JuMP.objective_value(m)) < 1e-3 @test abs(JuMP.objective_value(m)) < 1e-3
end end
@ -168,10 +173,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd, wd;
upper_bound=UB, upper_bound = UB,
halfradius=2, halfradius = 2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9) optimizer = scs_optimizer(; accel = 50, alpha = 1.9),
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test !certified @test !certified
@ -183,7 +188,8 @@ end
n = 4 n = 4
G = MatrixGroups.SpecialLinearGroup{n}(Int8) G = MatrixGroups.SpecialLinearGroup{n}(Int8)
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) @info "running tests for" G
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -198,6 +204,7 @@ end
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
) )
@info wd
sq, adj, op = PropertyT.SqAdjOp(RG, n) sq, adj, op = PropertyT.SqAdjOp(RG, n)
@ -208,10 +215,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd, wd;
upper_bound=UB, upper_bound = UB,
halfradius=2, halfradius = 2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9) optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test certified @test certified
@ -225,10 +232,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd, wd;
upper_bound=UB, upper_bound = UB,
halfradius=2, halfradius = 2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9) optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test certified @test certified
@ -242,10 +249,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd, wd;
upper_bound=UB, upper_bound = UB,
halfradius=2, halfradius = 2,
optimizer=cosmo_optimizer(accel=50, alpha=1.9) optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test !certified @test !certified

168
test/Chevalley.jl Normal file
View File

@ -0,0 +1,168 @@
countmap(v) = countmap(identity, v)
function countmap(f, v)
counts = Dict{eltype(f(first(v))),Int}()
for x in v
fx = f(x)
counts[fx] = get!(counts, fx, 0) + 1
end
return counts
end
@testset "classify_root_system" begin
α = PropertyT.Roots.Root([1, -1, 0])
β = PropertyT.Roots.Root([0, 1, -1])
γ = PropertyT.Roots.Root([2, 0, 0])
@test PropertyT.Roots.classify_root_system(α, β, (false, false)) == :A₂
@test PropertyT.Roots.classify_root_system(α, γ, (false, true)) == :C₂
@test PropertyT.Roots.classify_root_system(β, γ, (false, true)) ==
Symbol("A₁×C₁")
end
@testset "Exceptional root systems" begin
@testset "F4" begin
F4 = let Σ = PermutationGroups.PermGroup(perm"(1,2,3,4)", perm"(1,2)")
long = let x = (1, 1, 0, 0) .// 1
PropertyT.Roots.Root.(
union(
(x^g for g in Σ),
((x .* (-1, 1, 1, 1))^g for g in Σ),
((-1 .* x)^g for g in Σ),
),
)
end
short = let x = (1, 0, 0, 0) .// 1
PropertyT.Roots.Root.(
union((x^g for g in Σ), ((-1 .* x)^g for g in Σ))
)
end
signs = collect(Iterators.product(fill([-1, +1], 4)...))
halfs = let x = (1, 1, 1, 1) .// 2
PropertyT.Roots.Root.(union(x .* sgn for sgn in signs))
end
union(long, short, halfs)
end
@test length(F4) == 48
a = F4[1]
@test isapprox(PropertyT.Roots.ℓ₂length(a), sqrt(2))
b = F4[6]
@test isapprox(PropertyT.Roots.ℓ₂length(b), sqrt(2))
c = a + b
@test isapprox(PropertyT.Roots.ℓ₂length(c), 2.0)
@test PropertyT.Roots.classify_root_system(b, c, (false, true)) == :C₂
long = F4[findfirst(r -> PropertyT.Roots.ℓ₂length(r) == sqrt(2), F4)]
short = F4[findfirst(r -> PropertyT.Roots.ℓ₂length(r) == 1.0, F4)]
subtypes = Set([:C₂, :A₂, Symbol("A₁×C₁")])
let Ω = F4, α = long
counts = countmap([
PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
γ in Ω if !PropertyT.Roots.isproportional(α, γ)
])
@test Set(keys(counts)) == subtypes
d, r = divrem(counts[:C₂], 6)
@test r == 0 && d == 3
d, r = divrem(counts[:A₂], 4)
@test r == 0 && d == 4
end
let Ω = F4, α = short
counts = countmap([
PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
γ in Ω if !PropertyT.Roots.isproportional(α, γ)
])
@test Set(keys(counts)) == subtypes
d, r = divrem(counts[:C₂], 6)
@test r == 0 && d == 3
d, r = divrem(counts[:A₂], 4)
@test r == 0 && d == 4
end
end
@testset "E6-7-8 exceptional root systems" begin
E8 =
let Σ = PermutationGroups.PermGroup(
perm"(1,2,3,4,5,6,7,8)",
perm"(1,2)",
)
long = let x = (1, 1, 0, 0, 0, 0, 0, 0) .// 1
PropertyT.Roots.Root.(
union(
(x^g for g in Σ),
((x .* (-1, 1, 1, 1, 1, 1, 1, 1))^g for g in Σ),
((-1 .* x)^g for g in Σ),
),
)
end
signs = collect(
p for p in Iterators.product(fill([-1, +1], 8)...) if
iseven(count(==(-1), p))
)
halfs = let x = (1, 1, 1, 1, 1, 1, 1, 1) .// 2
rts = unique(PropertyT.Roots.Root(x .* sgn) for sgn in signs)
end
union(long, halfs)
end
subtypes = Set([:A₂, Symbol("A₁×A₁")])
@testset "E8" begin
@test length(E8) == 240
@test all(r -> PropertyT.Roots.ℓ₂length(r) ≈ sqrt(2), E8)
let Ω = E8, α = first(Ω)
counts = countmap([
PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
γ in Ω if !PropertyT.Roots.isproportional(α, γ)
])
@test Set(keys(counts)) == subtypes
d, r = divrem(counts[:A₂], 4)
@test r == 0 && d == 28
end
end
@testset "E7" begin
E7 = filter(r -> iszero(sum(r.coord)), E8)
@test length(E7) == 126
let Ω = E7, α = first(Ω)
counts = countmap([
PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
γ in Ω if !PropertyT.Roots.isproportional(α, γ)
])
@test Set(keys(counts)) == subtypes
d, r = divrem(counts[:A₂], 4)
@test r == 0 && d == 16
end
end
@testset "E6" begin
E6 = filter(
r -> r.coord[end] == r.coord[end-1] == r.coord[end-2],
E8,
)
@test length(E6) == 72
let Ω = E6, α = first(Ω)
counts = countmap([
PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
γ in Ω if !PropertyT.Roots.isproportional(α, γ)
])
@test Set(keys(counts)) == subtypes
d, r = divrem(counts[:A₂], 4)
@info d, r
@test r == 0 && d == 10
end
end
end
end

56
test/actions.jl
View File

@ -3,7 +3,7 @@ function test_action(basis, group, act)
return @testset "action definition" begin return @testset "action definition" begin
@test all(basis) do b @test all(basis) do b
e = one(group) e = one(group)
action(act, e, b) == b return action(act, e, b) == b
end end
a = let a = rand(basis) a = let a = rand(basis)
@ -30,12 +30,12 @@ function test_action(basis, group, act)
@test all([(g, h) for g in group for h in group]) do (g, h) @test all([(g, h) for g in group for h in group]) do (g, h)
x = action(act, h, action(act, g, a)) x = action(act, h, action(act, g, a))
y = action(act, g * h, a) y = action(act, g * h, a)
x == y return x == y
end end
if act isa SymbolicWedderburn.ByPermutations if act isa SymbolicWedderburn.ByPermutations
@test all(basis) do b @test all(basis) do b
action(act, g, b) ∈ basis && action(act, h, b) ∈ basis return action(act, g, b) ∈ basis && action(act, h, b) ∈ basis
end end
end end
end end
@ -44,21 +44,18 @@ end
## Testing ## Testing
@testset "Actions on SL(3,ℤ)" begin @testset "Actions on SL(3,ℤ)" begin
n = 3 n = 3
SL = MatrixGroups.SpecialLinearGroup{n}(Int8) SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true) RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2)
@testset "Permutation action" begin @testset "Permutation action" begin
Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
ΓpA = PropertyT.action_by_conjugation(SL, Γ) ΓpA = PropertyT.action_by_conjugation(SL, Γ)
test_action(basis(RSL), Γ, ΓpA) test_action(basis(RSL), Γ, ΓpA)
@testset "mps is successful" begin @testset "mps is successful" begin
charsΓ = charsΓ =
SymbolicWedderburn.Character{ SymbolicWedderburn.Character{
Rational{Int}, Rational{Int},
@ -66,9 +63,9 @@ end
RΓ = SymbolicWedderburn._group_algebra(Γ) RΓ = SymbolicWedderburn._group_algebra(Γ)
@time mps, simple = @time mps, ranks =
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
@test all(simple) @test all(isone, ranks)
end end
@testset "Wedderburn decomposition" begin @testset "Wedderburn decomposition" begin
@ -77,17 +74,17 @@ end
Γ, Γ,
ΓpA, ΓpA,
basis(RSL), basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]) StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
) )
@test length(invariant_vectors(wd)) == 918 @test length(invariant_vectors(wd)) == 918
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [40, 23, 18] @test SymbolicWedderburn.size.(direct_summands(wd), 1) ==
[40, 23, 18]
@test all(issimple, direct_summands(wd)) @test all(issimple, direct_summands(wd))
end end
end end
@testset "Wreath action" begin @testset "Wreath action" begin
Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
end end
@ -97,7 +94,6 @@ end
test_action(basis(RSL), Γ, ΓpA) test_action(basis(RSL), Γ, ΓpA)
@testset "mps is successful" begin @testset "mps is successful" begin
charsΓ = charsΓ =
SymbolicWedderburn.Character{ SymbolicWedderburn.Character{
Rational{Int}, Rational{Int},
@ -105,9 +101,9 @@ end
RΓ = SymbolicWedderburn._group_algebra(Γ) RΓ = SymbolicWedderburn._group_algebra(Γ)
@time mps, simple = @time mps, ranks =
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
@test all(simple) @test all(isone, ranks)
end end
@testset "Wedderburn decomposition" begin @testset "Wedderburn decomposition" begin
@ -116,32 +112,30 @@ end
Γ, Γ,
ΓpA, ΓpA,
basis(RSL), basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]) StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
) )
@test length(invariant_vectors(wd)) == 247 @test length(invariant_vectors(wd)) == 247
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [14, 9, 6, 14, 12] @test SymbolicWedderburn.size.(direct_summands(wd), 1) ==
[14, 9, 6, 14, 12]
@test all(issimple, direct_summands(wd)) @test all(issimple, direct_summands(wd))
end end
end end
end end
@testset "Actions on SAut(F4)" begin @testset "Actions on SAut(F4)" begin
n = 4 n = 4
SAutFn = SpecialAutomorphismGroup(FreeGroup(n)) SAutFn = SpecialAutomorphismGroup(FreeGroup(n))
RSAutFn, S, sizes = PropertyT.group_algebra(SAutFn, halfradius=1, twisted=true) RSAutFn, S, sizes = PropertyT.group_algebra(SAutFn; halfradius = 1)
@testset "Permutation action" begin @testset "Permutation action" begin
Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
ΓpA = PropertyT.action_by_conjugation(SAutFn, Γ) ΓpA = PropertyT.action_by_conjugation(SAutFn, Γ)
test_action(basis(RSAutFn), Γ, ΓpA) test_action(basis(RSAutFn), Γ, ΓpA)
@testset "mps is successful" begin @testset "mps is successful" begin
charsΓ = charsΓ =
SymbolicWedderburn.Character{ SymbolicWedderburn.Character{
Rational{Int}, Rational{Int},
@ -149,9 +143,9 @@ end
RΓ = SymbolicWedderburn._group_algebra(Γ) RΓ = SymbolicWedderburn._group_algebra(Γ)
@time mps, simple = @time mps, ranks =
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
@test all(simple) @test all(isone, ranks)
end end
@testset "Wedderburn decomposition" begin @testset "Wedderburn decomposition" begin
@ -160,17 +154,17 @@ end
Γ, Γ,
ΓpA, ΓpA,
basis(RSAutFn), basis(RSAutFn),
StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]]) StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]]),
) )
@test length(invariant_vectors(wd)) == 93 @test length(invariant_vectors(wd)) == 93
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [4, 8, 5, 4] @test SymbolicWedderburn.size.(direct_summands(wd), 1) ==
[4, 8, 5, 4]
@test all(issimple, direct_summands(wd)) @test all(issimple, direct_summands(wd))
end end
end end
@testset "Wreath action" begin @testset "Wreath action" begin
Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
end end
@ -180,7 +174,6 @@ end
test_action(basis(RSAutFn), Γ, ΓpA) test_action(basis(RSAutFn), Γ, ΓpA)
@testset "mps is successful" begin @testset "mps is successful" begin
charsΓ = charsΓ =
SymbolicWedderburn.Character{ SymbolicWedderburn.Character{
Rational{Int}, Rational{Int},
@ -188,9 +181,9 @@ end
RΓ = SymbolicWedderburn._group_algebra(Γ) RΓ = SymbolicWedderburn._group_algebra(Γ)
@time mps, simple = @time mps, ranks =
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
@test all(simple) @test all(isone, ranks)
end end
@testset "Wedderburn decomposition" begin @testset "Wedderburn decomposition" begin
@ -199,11 +192,12 @@ end
Γ, Γ,
ΓpA, ΓpA,
basis(RSAutFn), basis(RSAutFn),
StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]]) StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]]),
) )
@test length(invariant_vectors(wd)) == 18 @test length(invariant_vectors(wd)) == 18
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [1, 1, 2, 2, 1, 2, 2, 1] @test SymbolicWedderburn.size.(direct_summands(wd), 1) ==
[1, 1, 2, 2, 1, 2, 2, 1]
@test all(issimple, direct_summands(wd)) @test all(issimple, direct_summands(wd))
end end
end end

35
test/check_positivity.jl
View File

@ -1,6 +1,12 @@
function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer) function check_positivity(
elt,
unit;
upper_bound = Inf,
halfradius = 2,
optimizer,
)
@time sos_problem = @time sos_problem =
PropertyT.sos_problem_primal(elt, unit, upper_bound=upper_bound) PropertyT.sos_problem_primal(elt, unit; upper_bound = upper_bound)
status, _ = PropertyT.solve(sos_problem, optimizer) status, _ = PropertyT.solve(sos_problem, optimizer)
P = JuMP.value.(sos_problem[:P]) P = JuMP.value.(sos_problem[:P])
@ -9,16 +15,23 @@ function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer)
elt, elt,
unit, unit,
JuMP.objective_value(sos_problem), JuMP.objective_value(sos_problem),
Q, Q;
halfradius=halfradius, halfradius = halfradius,
) )
return status, certified, λ_cert return status, certified, λ_cert
end end
function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimizer) function check_positivity(
elt,
unit,
wd;
upper_bound = Inf,
halfradius = 2,
optimizer,
)
@assert aug(elt) == aug(unit) == 0 @assert aug(elt) == aug(unit) == 0
@time sos_problem, Ps = @time sos_problem, Ps =
PropertyT.sos_problem_primal(elt, unit, wd, upper_bound=upper_bound) PropertyT.sos_problem_primal(elt, unit, wd; upper_bound = upper_bound)
@time status, _ = PropertyT.solve(sos_problem, optimizer) @time status, _ = PropertyT.solve(sos_problem, optimizer)
@ -29,13 +42,7 @@ function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimize
λ = JuMP.value(sos_problem[:λ]) λ = JuMP.value(sos_problem[:λ])
certified, λ_cert = PropertyT.certify_solution( certified, λ_cert =
elt, PropertyT.certify_solution(elt, unit, λ, Q; halfradius = halfradius)
unit,
λ,
Q,
halfradius=halfradius
)
return status, certified, λ_cert return status, certified, λ_cert
end end

14
test/constratint_matrices.jl
View File

@ -1,12 +1,17 @@
@testset "ConstraintMatrix" begin @testset "ConstraintMatrix" begin
@test PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π) isa AbstractMatrix @test PropertyT.ConstraintMatrix{Float64}(
[-1, 2, -1, 1, 4, 2, 6],
3,
2,
π,
) isa AbstractMatrix
cm = PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π) cm = PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π)
@test cm == Float64[ @test cm == Float64[
-π π -π π
2π 0.0 2π 0
0.0 π 0 π
] ]
@test collect(PropertyT.nzpairs(cm)) == [ @test collect(PropertyT.nzpairs(cm)) == [
@ -18,4 +23,7 @@
1 => -3.141592653589793 1 => -3.141592653589793
1 => -3.141592653589793 1 => -3.141592653589793
] ]
@test PropertyT.ConstraintMatrix{Float64}([-9:-1; 1:9], 3, 3, 1.0) ==
zeros(3, 3)
end end

26
test/graded_adj.jl
View File

@ -1,10 +1,9 @@
@testset "Adj via grading" begin @testset "Adj via grading" begin
@testset "SL(n,Z) & Aut(F₄)" begin @testset "SL(n,Z) & Aut(F₄)" begin
n = 4 n = 4
halfradius = 1 halfradius = 1
SL = MatrixGroups.SpecialLinearGroup{n}(Int8) SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=halfradius, twisted=true) RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = halfradius)
Δ = RSL(length(S)) - sum(RSL(s) for s in S) Δ = RSL(length(S)) - sum(RSL(s) for s in S)
@ -22,10 +21,9 @@
@test PropertyT.Adj(Δs, :A₂) == adj @test PropertyT.Adj(Δs, :A₂) == adj
@test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op @test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
halfradius = 1 halfradius = 1
G = SpecialAutomorphismGroup(FreeGroup(n)) G = SpecialAutomorphismGroup(FreeGroup(n))
RG, S, sizes = PropertyT.group_algebra(G, halfradius=halfradius, twisted=true) RG, S, sizes = PropertyT.group_algebra(G; halfradius = halfradius)
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -44,7 +42,6 @@
@test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op @test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
end end
@testset "Symplectic group" begin @testset "Symplectic group" begin
@testset "Sp2(ℤ)" begin @testset "Sp2(ℤ)" begin
genus = 2 genus = 2
@ -52,7 +49,8 @@
SpN = MatrixGroups.SymplecticGroup{2genus}(Int8) SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true) RSpN, S_sp, sizes_sp =
PropertyT.group_algebra(SpN; halfradius = halfradius)
Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -73,7 +71,8 @@
SpN = MatrixGroups.SymplecticGroup{2genus}(Int8) SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true) RSpN, S_sp, sizes_sp =
PropertyT.group_algebra(SpN; halfradius = halfradius)
Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -86,9 +85,14 @@
end end
@testset "Adj numerics for genus=$genus" begin @testset "Adj numerics for genus=$genus" begin
all_subtypes = ( all_subtypes = (
:A₁, :C₁, Symbol("A₁×A₁"), Symbol("C₁×C₁"), Symbol("A₁×C₁"), :A₂, :C₂ :A₁,
:C₁,
Symbol("A₁×A₁"),
Symbol("C₁×C₁"),
Symbol("A₁×C₁"),
:A₂,
:C₂,
) )
@test PropertyT.Adj(Δs, :A₂)[one(SpN)] == 384 @test PropertyT.Adj(Δs, :A₂)[one(SpN)] == 384
@ -101,8 +105,8 @@
@test isinteger(PropertyT.Adj(Δs, subtype)[one(SpN)] / 16) @test isinteger(PropertyT.Adj(Δs, subtype)[one(SpN)] / 16)
end end
end end
@test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) == Δ^2 @test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) ==
Δ^2
end end
end end
end end

36
test/optimizers.jl
View File

@ -4,12 +4,12 @@ import JuMP
import SCS import SCS
function scs_optimizer(; function scs_optimizer(;
accel=10, accel = 10,
alpha=1.5, alpha = 1.5,
eps=1e-9, eps = 1e-9,
max_iters=100_000, max_iters = 100_000,
verbose=true, verbose = true,
linear_solver=SCS.DirectSolver linear_solver = SCS.DirectSolver,
) )
return JuMP.optimizer_with_attributes( return JuMP.optimizer_with_attributes(
SCS.Optimizer, SCS.Optimizer,
@ -28,27 +28,27 @@ end
import COSMO import COSMO
function cosmo_optimizer(; function cosmo_optimizer(;
accel=15, accel = 15,
alpha=1.6, alpha = 1.6,
eps=1e-9, eps = 1e-9,
max_iters=100_000, max_iters = 100_000,
verbose=true, verbose = true,
verbose_timing=verbose, verbose_timing = verbose,
decompose=false decompose = false,
) )
return JuMP.optimizer_with_attributes( return JuMP.optimizer_with_attributes(
COSMO.Optimizer, COSMO.Optimizer,
"accelerator" => COSMO.with_options( "accelerator" =>
COSMO.AndersonAccelerator, COSMO.with_options(COSMO.AndersonAccelerator; mem = max(accel, 2)),
mem=max(accel, 2)
),
"alpha" => alpha, "alpha" => alpha,
"decompose" => decompose, "decompose" => decompose,
"eps_abs" => eps, "eps_abs" => eps,
"eps_rel" => eps, "eps_rel" => eps,
"eps_prim_inf" => eps,
"eps_dual_inf" => eps,
"max_iter" => max_iters, "max_iter" => max_iters,
"verbose" => verbose, "verbose" => verbose,
"verbose_timing" => verbose_timing, "verbose_timing" => verbose_timing,
"check_termination" => 200, "check_termination" => 250,
) )
end end

63
test/quick_tests.jl
View File

@ -1,11 +1,12 @@
@testset "Quick tests" begin @testset "Quick tests" begin
@testset "SL(2,F₇)" begin @testset "SL(2,F₇)" begin
N = 2 N = 2
p = 7 p = 7
halfradius = 3 halfradius = 3
G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{p}) G = MatrixGroups.SpecialLinearGroup{N}(
RG, S, sizes = PropertyT.group_algebra(G, halfradius=3, twisted=true) SymbolicWedderburn.Characters.FiniteFields.GF{p},
)
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 3)
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -18,15 +19,15 @@
@testset "standard formulation" begin @testset "standard formulation" begin
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
unit, unit;
upper_bound=ub, upper_bound = ub,
halfradius=2, halfradius = 2,
optimizer=cosmo_optimizer( optimizer = cosmo_optimizer(;
eps=1e-7, eps = 1e-7,
max_iters=5_000, max_iters = 5_000,
accel=50, accel = 50,
alpha=1.95, alpha = 1.95,
) ),
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@ -34,14 +35,18 @@
@test λ_cert > 5857 // 10000 @test λ_cert > 5857 // 10000
m = PropertyT.sos_problem_dual(elt, unit) m = PropertyT.sos_problem_dual(elt, unit)
PropertyT.solve(m, cosmo_optimizer( PropertyT.solve(
eps=1e-7, m,
max_iters=10_000, cosmo_optimizer(;
accel=50, eps = 1e-7,
alpha=1.95, max_iters = 10_000,
)) accel = 50,
alpha = 1.95,
),
)
@test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL) @test JuMP.termination_status(m) in
(JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
@test JuMP.objective_value(m) ≈ λ_cert atol = 1e-2 @test JuMP.objective_value(m) ≈ λ_cert atol = 1e-2
end end
@ -55,20 +60,22 @@
Σ, Σ,
act, act,
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[halfradius]]), StarAlgebras.Basis{UInt16}(
@view basis(RG)[1:sizes[halfradius]]
),
) )
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
unit, unit,
wd, wd;
upper_bound=ub, upper_bound = ub,
halfradius=2, halfradius = 2,
optimizer=cosmo_optimizer( optimizer = cosmo_optimizer(;
eps=1e-7, eps = 1e-7,
max_iters=10_000, max_iters = 10_000,
accel=50, accel = 50,
alpha=1.9, alpha = 1.9,
), ),
) )

10
test/roots.jl Normal file
View File

@ -0,0 +1,10 @@
using PropertyT.Roots
@testset "Roots" begin
@test Roots.Root{3,Int}([1, 2, 3]) isa Roots.AbstractRoot{}
@test Roots.Root([1, 2, 3]) isa Roots.AbstractRoot{3,Int}
# io
r = Roots.Root{3,Int}([1, 2, 3])
@test contains(sprint(show, MIME"text/plain"(), r), "of length √14\n")
r = Roots.Root{3,Int}([1, 2, 2])
@test contains(sprint(show, MIME"text/plain"(), r), "of length 3\n")
end

2
test/runtests.jl
View File

@ -24,6 +24,8 @@ if haskey(ENV, "FULL_TEST") || haskey(ENV, "CI")
include("1712.07167.jl") include("1712.07167.jl")
include("1812.03456.jl") include("1812.03456.jl")
include("roots.jl")
include("graded_adj.jl") include("graded_adj.jl")
include("Chevalley.jl")
end end
end end