kalmar/PropertyT.jl
1
0
Fork 0
mirror of https://github.com/kalmarek/PropertyT.jl.git synced 2024-11-19 15:25:29 +01:00
Code Issues Releases Wiki Activity

Merge pull request #9 from kalmarek/mk/redesign

Browse Source 
Redesign + graded Adj
This commit is contained in:
Marek Kaluba 2022-11-08 16:09:03 +01:00 committed by GitHub
commit 2b20a5bb34
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
30 changed files with 2138 additions and 1737 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
.github/workflows/CompatHelper.yml vendored Normal file
View File

@ -0,0 +1,16 @@
name: CompatHelper
on:
schedule:
- cron: '00 00 * * *'
workflow_dispatch:
jobs:
CompatHelper:
runs-on: ubuntu-latest
steps:
- name: Pkg.add("CompatHelper")
run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
- name: CompatHelper.main()
env:
GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
COMPATHELPER_PRIV: ${{ secrets.COMPATHELPER_PRIV }} # optional
run: julia -e 'using CompatHelper; CompatHelper.main()'

10
.github/workflows/TagBot.yml vendored
View File

@ -1,11 +1,17 @@
name: TagBot name: TagBot
on: on:
schedule: issue_comment: # THIS BIT IS NEW
- cron: 0 * * * * types:
- created
workflow_dispatch:
jobs: jobs:
TagBot: TagBot:
# THIS 'if' LINE IS NEW
if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
# NOTHING BELOW HAS CHANGED
runs-on: ubuntu-latest runs-on: ubuntu-latest
steps: steps:
- uses: JuliaRegistries/TagBot@v1 - uses: JuliaRegistries/TagBot@v1
with: with:
token: ${{ secrets.GITHUB_TOKEN }} token: ${{ secrets.GITHUB_TOKEN }}
# ssh: ${{ secrets.DOCUMENTER_KEY }}

26
.github/workflows/runtests.yml → .github/workflows/ci.yml vendored
View File

@ -1,11 +1,6 @@
name: CI name: CI
on: on:
push: - pull_request
branches:
- master
pull_request:
branches:
- master
jobs: jobs:
test: test:
name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
@ -13,22 +8,32 @@ jobs:
strategy: strategy:
matrix: matrix:
version: version:
- '1.3' - '1.6'
- '1.5' - '1'
- 'nightly' - 'nightly'
os: os:
- ubuntu-latest - ubuntu-latest
- macOS-latest
arch: arch:
- x64 - x64
allow_failures:
- julia: nightly
fail-fast: false fail-fast: false
steps: steps:
- uses: actions/checkout@v2 - uses: actions/checkout@v2
- uses: julia-actions/setup-julia@v1 - uses: julia-actions/setup-julia@v1
with: with:
version: ${{ matrix.version }} version: ${{ matrix.version }}
arch: ${{ matrix.arch }} arch: ${{ matrix.arch }}
- uses: actions/cache@v1
env:
cache-name: cache-artifacts
with:
path: ~/.julia/artifacts
key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
restore-keys: |
${{ runner.os }}-test-${{ env.cache-name }}-
${{ runner.os }}-test-
${{ runner.os }}-
- uses: julia-actions/julia-buildpkg@latest - uses: julia-actions/julia-buildpkg@latest
- uses: julia-actions/julia-runtest@latest - uses: julia-actions/julia-runtest@latest
- uses: julia-actions/julia-processcoverage@v1 - uses: julia-actions/julia-processcoverage@v1
@ -38,4 +43,3 @@ jobs:
flags: unittests flags: unittests
name: codecov-umbrella name: codecov-umbrella
fail_ci_if_error: false fail_ci_if_error: false
token: ${{ secrets.CODECOV_TOKEN }}

32
Project.toml
View File

@ -1,33 +1,35 @@
name = "PropertyT" name = "PropertyT"
uuid = "03b72c93-0167-51e2-8a1e-eb4ff1fb940d" uuid = "03b72c93-0167-51e2-8a1e-eb4ff1fb940d"
authors = ["Marek Kaluba <kalmar@amu.edu.pl>"] authors = ["Marek Kaluba <kalmar@amu.edu.pl>"]
version = "0.3.2" version = "0.4.0"
[deps] [deps]
AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
GroupRings = "0befed6a-bd73-11e8-1e41-a1190947c2f5"
Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557" Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572" JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
[compat] [compat]
AbstractAlgebra = "^0.10.0" COSMO = "0.8"
GroupRings = "^0.3.2" Groups = "0.7"
Groups = "^0.5.0" IntervalArithmetic = "0.20"
IntervalArithmetic = "^0.16.0" JuMP = "1.3"
JLD = "^0.9.0" SCS = "1.1.0"
JuMP = "^0.20.0" StaticArrays = "1"
SCS = "^0.7.0" SymbolicWedderburn = "0.3.1"
julia = "^1.3.0" julia = "1.6"
[extras] [extras]
COSMO = "1e616198-aa4e-51ec-90a2-23f7fbd31d8d"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[targets] [targets]
test = ["Test", "SCS"] scripts = ["Dates", "Logging", "Serialization", "SCS", "COSMO"]
test = ["Test", "SCS", "COSMO"]

329
src/1712.07167.jl
View File

@ -1,329 +0,0 @@
###############################################################################
# Settings and filenames
abstract type Settings end
struct Naive{El} <: Settings
name::String
G::Union{Group, NCRing}
S::Vector{El}
halfradius::Int
upper_bound::Float64
solver::JuMP.OptimizerFactory
force_compute::Bool
end
struct Symmetrized{El} <: Settings
name::String
G::Union{Group, NCRing}
S::Vector{El}
autS::Group
halfradius::Int
upper_bound::Float64
solver::JuMP.OptimizerFactory
force_compute::Bool
end
function Settings(name::String,
G::Union{Group, NCRing}, S::AbstractVector{El}, solver::JuMP.OptimizerFactory;
halfradius=2, upper_bound=1.0, force_compute=false) where El <: Union{GroupElem, NCRingElem}
return Naive(name, G, S, halfradius, upper_bound, solver, force_compute)
end
function Settings(name::String,
G::Union{Group, NCRing}, S::AbstractVector{El}, autS::Group, solver::JuMP.OptimizerFactory;
halfradius=2, upper_bound=1.0, force_compute=false) where El <: Union{GroupElem, NCRingElem}
return Symmetrized(name, G, S, autS, halfradius, upper_bound, solver, force_compute)
end
suffix(s::Settings) = "$(s.upper_bound)"
prepath(s::Settings) = s.name
fullpath(s::Settings) = joinpath(prepath(s), suffix(s))
filename(sett::Settings, s::Symbol; kwargs...) = filename(sett, Val{s}; kwargs...)
filename(sett::Settings, ::Type{Val{:fulllog}}; kwargs...) =
filename(fullpath(sett), "full", "log", suffix=Dates.now(); kwargs...)
filename(sett::Settings, ::Type{Val{:solverlog}}; kwargs...) =
filename(fullpath(sett), "solver", "log", suffix=Dates.now(); kwargs...)
filename(sett::Settings, ::Type{Val{:Δ}}; kwargs...) =
filename(prepath(sett), "delta", "jld"; kwargs...)
filename(sett::Settings, ::Type{Val{:BlockDecomposition}}; kwargs...) =
filename(prepath(sett), "BlockDecomposition", "jld"; kwargs...)
filename(sett::Settings, ::Type{Val{:solution}}; kwargs...) =
filename(fullpath(sett), "solution", "jld"; kwargs...)
function filename(sett::Settings, ::Type{Val{:warmstart}}; kwargs...)
filename(fullpath(sett), "warmstart", "jld"; kwargs...)
end
function filename(path::String, name, extension; prefix=nothing, suffix=nothing)
pre = isnothing(prefix) ? "" : "$(prefix)_"
suf = isnothing(suffix) ? "" : "_$(suffix)"
return joinpath(path, "$pre$name$suf.$extension")
end
###############################################################################
# Approximation by SOS (logged & warmstarted)
function warmstart(sett::Settings)
warmstart_fname = filename(sett, :warmstart)
try
ws = load(warmstart_fname, "warmstart")
@info "Loaded $warmstart_fname."
return ws
catch ex
@warn "$(ex.msg). Could not provide a warmstart to the solver."
return nothing
end
end
function approximate_by_SOS(sett::Naive,
elt::GroupRingElem, orderunit::GroupRingElem;
solverlog=tempname()*".log")
isdir(fullpath(sett)) || mkpath(fullpath(sett))
@info "Creating SDP problem..."
SDP_problem = SOS_problem_primal(elt, orderunit, upper_bound=sett.upper_bound)
@info Base.repr(SDP_problem)
@info "Logging solver's progress into $solverlog"
ws = warmstart(sett)
@time status, ws = PropertyT.solve(solverlog, SDP_problem, sett.solver, ws)
@info "Optimization finished:" status
P = value.(SDP_problem[:P])
λ = value(SDP_problem[:λ])
if any(isnan, P)
@warn "The solution seems to contain NaNs. Not overriding warmstart.jld"
else
save(filename(sett, :warmstart),
"warmstart", (ws.primal, ws.dual, ws.slack),
"P", P,
"λ", λ)
end
save(filename(sett, :warmstart, suffix=Dates.now()),
"warmstart", (ws.primal, ws.dual, ws.slack),
"P", P,
"λ", λ)
return λ, P
end
function approximate_by_SOS(sett::Symmetrized,
elt::GroupRingElem, orderunit::GroupRingElem;
solverlog=tempname()*".log")
isdir(fullpath(sett)) || mkpath(fullpath(sett))
orbit_data = try
orbit_data = load(filename(sett, :BlockDecomposition), "BlockDecomposition")
@info "Loaded orbit data."
orbit_data
catch ex
@warn ex.msg
GroupRings.hasbasis(parent(orderunit)) ||
throw("You need to define basis of Group Ring to compute orbit decomposition!")
@info "Computing orbit and Wedderburn decomposition..."
orbit_data = BlockDecomposition(parent(orderunit), sett.autS)
save(filename(sett, :BlockDecomposition), "BlockDecomposition", orbit_data)
orbit_data
end
orbit_data = decimate(orbit_data)
@info "Creating SDP problem..."
SDP_problem, varP = SOS_problem_primal(elt, orderunit, orbit_data, upper_bound=sett.upper_bound)
@info Base.repr(SDP_problem)
@info "Logging solver's progress into $solverlog"
ws = warmstart(sett)
@time status, ws = PropertyT.solve(solverlog, SDP_problem, sett.solver, ws)
@info "Optimization finished:" status
λ = value(SDP_problem[:λ])
Ps = [value.(P) for P in varP]
if any(any(isnan, P) for P in Ps)
@warn "The solution seems to contain NaNs. Not overriding warmstart.jld"
else
save(filename(sett, :warmstart),
"warmstart", (ws.primal, ws.dual, ws.slack),
"Ps", Ps,
"λ", λ)
end
save(filename(sett, :warmstart, suffix=Dates.now()),
"warmstart", (ws.primal, ws.dual, ws.slack),
"Ps", Ps,
"λ", λ)
@info "Reconstructing P..."
@time P = reconstruct(Ps, orbit_data)
return λ, P
end
###############################################################################
# Checking solution
function certify_SOS_decomposition(elt::GroupRingElem, orderunit::GroupRingElem,
λ::Number, Q::AbstractMatrix; R::Int=2)
separator = "-"^76
@info "$separator\nChecking in floating-point arithmetic..." λ
eoi = elt - λ*orderunit
@info("Computing sum of squares decomposition...")
@time residual = eoi - compute_SOS(parent(eoi), augIdproj(Q))
L1_norm = norm(residual,1)
floatingpoint_λ = λ - 2.0^(2ceil(log2(R)))*L1_norm
info_strs = ["Numerical metrics of the obtained SOS:",
"ɛ(elt - λu - ∑ξᵢ*ξᵢ) ≈ $(aug(residual))",
"‖elt - λu - ∑ξᵢ*ξᵢ‖₁ ≈ $(L1_norm)",
"Floating point (NOT certified) λ ≈"]
@info join(info_strs, "\n") floatingpoint_λ
if floatingpoint_λ ≤ 0
return floatingpoint_λ
end
λ = @interval(λ)
info_strs = [separator,
"Checking in interval arithmetic...",
"λ ∈ $λ"]
@info(join(info_strs, "\n"))
eoi = elt - λ*orderunit
@info("Projecting columns of Q to the augmentation ideal...")
@time Q, check = augIdproj(Interval, Q)
@info "Checking that sum of every column contains 0.0..." check_augmented=check
check || @error("The following numbers are meaningless!")
@info("Computing sum of squares decomposition...")
@time residual = eoi - compute_SOS(parent(eoi), Q)
L1_norm = norm(residual,1)
certified_λ = λ - 2.0^(2ceil(log2(R)))*L1_norm
info_strs = ["Numerical metrics of the obtained SOS:",
"ɛ(elt - λu - ∑ξᵢ*ξᵢ) ∈ $(aug(residual))",
"‖elt - λu - ∑ξᵢ*ξᵢ‖₁ ∈ $(L1_norm)",
"Interval aritmetic (certified) λ ∈"]
@info join(info_strs, "\n") certified_λ
return certified_λ
end
function spectral_gap(Δ::GroupRingElem, λ::Number, Q::AbstractMatrix; R::Int=2)
@info "elt = Δ², u = Δ"
return certify_SOS_decomposition(Δ^2, Δ, λ, Q, R=R)
end
###############################################################################
# Interpreting the numerical results
Kazhdan_constant(λ::Number, N::Integer) = sqrt(2*λ/N)
Kazhdan_constant(λ::Interval, N::Integer) = IntervalArithmetic.inf(sqrt(2*λ/N))
function check_property_T(sett::Settings)
@info sett
certified_sgap = spectral_gap(sett)
return interpret_results(sett, certified_sgap)
end
function Base.show(io::IO, sett::Settings)
info_strs = ["PropertyT Settings:",
"Group: $(sett.name)",
"Upper bound for λ: $(sett.upper_bound), on halfradius $(sett.halfradius).",
"Force computations: $(sett.force_compute);",
"Results will be stored in ./$(PropertyT.prepath(sett));",
"Solver: $(typeof(sett.solver()))",
"Solvers options: "]
append!(info_strs, [rpad(" $k", 30)* "→ \t$v" for (k,v) in sett.solver().options])
push!(info_strs, "="^76)
print(io, join(info_strs, "\n"))
end
function interpret_results(name::String, sgap::Number, N::Integer)
if sgap > 0
κ = Kazhdan_constant(sgap, N)
@info "κ($name, S) ≥ $κ: Group HAS property (T)!"
return true
end
info_strs = [
"The certified lower bound on the spectral gap is negative:",
"λ($name, S) ≥ 0.0 > $sgap",
"This tells us nothing about property (T)",
]
@info join(info_strs, "\n")
return false
end
interpret_results(sett::Settings, sgap::Number) =
interpret_results(sett.name, sgap, length(sett.S))
function spectral_gap(sett::Settings)
fp = PropertyT.fullpath(sett)
isdir(fp) || mkpath(fp)
Δ = try
Δ = loadGRElem(filename(sett,:Δ), sett.G)
@info "Loaded precomputed Δ."
Δ
catch ex
@warn ex.msg
@info "Computing Δ..."
Δ = Laplacian(sett.S, sett.halfradius)
saveGRElem(filename(sett, :Δ), Δ)
Δ
end
function compute(sett, Δ)
@info "Computing λ and P..."
λ, P = approximate_by_SOS(sett, Δ^2, Δ;
solverlog=filename(sett, :solverlog))
save(filename(sett, :solution), "λ", λ, "P", P)
λ < 0 && @warn "Solver did not produce a valid solution!"
return λ, P
end
if sett.force_compute
λ, P = compute(sett, Δ)
else
λ, P =try
λ, P = load(filename(sett, :solution), "λ", "P")
@info "Loaded existing λ and P."
λ, P
catch ex
@warn ex.msg
compute(sett, Δ)
end
end
info_strs = ["Numerical metrics of matrix solution:",
"sum(P) = $(sum(P))",
"maximum(P) = $(maximum(P))",
"minimum(P) = $(minimum(P))"]
@info join(info_strs, "\n")
isapprox(eigvals(P), abs.(eigvals(P))) ||
@warn "The solution matrix doesn't seem to be positive definite!"
@time Q = real(sqrt(Symmetric( (P.+ P')./2 )))
certified_sgap = spectral_gap(Δ, λ, Q, R=sett.halfradius)
return certified_sgap
end

26
src/1812.03456.jl
View File

@ -1,26 +0,0 @@
indexing(n) = [(i,j) for i in 1:n for j in 1:n if i≠j]
function generating_set(G::AutGroup{N}, n=N) where N
rmuls = [Groups.transvection_R(i,j) for (i,j) in indexing(n)]
lmuls = [Groups.transvection_L(i,j) for (i,j) in indexing(n)]
gen_set = G.([rmuls; lmuls])
return [gen_set; inv.(gen_set)]
end
function EltaryMat(M::MatAlgebra, i::Integer, j::Integer, val=1)
@assert i ≠ j
@assert 1 ≤ i ≤ nrows(M)
@assert 1 ≤ j ≤ ncols(M)
m = one(M)
m[i,j] = val
return m
end
function generating_set(M::MatAlgebra, n=nrows(M))
elts = [EltaryMat(M, i,j) for (i,j) in indexing(n)]
return elem_type(M)[elts; inv.(elts)]
end
include("sqadjop.jl")

43
src/PropertyT.jl
View File

@ -1,26 +1,43 @@
__precompile__() __precompile__()
module PropertyT module PropertyT
using AbstractAlgebra
using LinearAlgebra using LinearAlgebra
using SparseArrays using SparseArrays
using Dates
using Groups using IntervalArithmetic
using GroupRings
using JLD
using JuMP using JuMP
import AbstractAlgebra: Group, NCRing using Groups
import Groups.GroupsCore
using SymbolicWedderburn
import SymbolicWedderburn.StarAlgebras
import SymbolicWedderburn.PermutationGroups
include("laplacians.jl") include("constraint_matrix.jl")
include("RGprojections.jl")
include("blockdecomposition.jl")
include("sos_sdps.jl") include("sos_sdps.jl")
include("checksolution.jl") include("certify.jl")
include("1712.07167.jl") include("sqadjop.jl")
include("1812.03456.jl")
include("roots.jl")
import .Roots
include("gradings.jl")
include("actions/actions.jl")
function group_algebra(G::Groups.Group, S=gens(G); halfradius::Integer, twisted::Bool)
S = union!(S, inv.(S))
@info "generating wl-metric ball of radius $(2halfradius)"
@time E, sizes = Groups.wlmetric_ball_serial(S, radius=2halfradius)
@info "sizes = $(sizes)"
@info "computing the *-algebra structure for G"
@time RG = StarAlgebras.StarAlgebra{twisted}(
G,
StarAlgebras.Basis{UInt32}(E),
(sizes[halfradius], sizes[halfradius]),
precompute=false,
)
return RG, S, sizes
end
end # module Property(T) end # module Property(T)

251
src/RGprojections.jl
View File

@ -1,251 +0,0 @@
module Projections
using AbstractAlgebra
using Groups
using GroupRings
export PermCharacter, DirectProdCharacter, rankOne_projections
###############################################################################
#
# Characters of Symmetric Group and DirectProduct
#
###############################################################################
abstract type AbstractCharacter end
struct PermCharacter <: AbstractCharacter
p::Generic.Partition
end
struct DirectProdCharacter{N, T<:AbstractCharacter} <: AbstractCharacter
chars::NTuple{N, T}
end
function (chi::DirectProdCharacter)(g::DirectPowerGroupElem)
res = 1
for (χ, elt) in zip(chi.chars, g.elts)
res *= χ(elt)
end
return res
end
function (chi::PermCharacter)(g::Generic.Perm)
R = AbstractAlgebra.partitionseq(chi.p)
p = Partition(Generic.permtype(g))
return Int(Generic.MN1inner(R, p, 1, Generic._charvalsTable))
end
AbstractAlgebra.dim(χ::PermCharacter) = dim(YoungTableau(χ.p))
for T in [PermCharacter, DirectProdCharacter]
@eval begin
function (chi::$T)(X::GroupRingElem)
RG = parent(X)
z = zero(eltype(X))
result = z
for i in 1:length(X.coeffs)
if X.coeffs[i] != z
result += chi(RG.basis[i])*X.coeffs[i]
end
end
return result
end
end
end
characters(G::Generic.SymmetricGroup) = (PermCharacter(p) for p in AllParts(G.n))
function characters(G::DirectPowerGroup{N}) where N
nfold_chars = Iterators.repeated(characters(G.group), N)
return (DirectProdCharacter(idx) for idx in Iterators.product(nfold_chars...))
end
###############################################################################
#
# Projections
#
###############################################################################
function central_projection(RG::GroupRing, chi::AbstractCharacter, T::Type=Rational{Int})
result = RG(zeros(T, length(RG.basis)))
dim = chi(one(RG.group))
ord = Int(order(RG.group))
for g in RG.basis
result[g] = convert(T, (dim//ord)*chi(g))
end
return result
end
function alternating_emb(RG::GroupRing{Gr,T}, V::Vector{T}, S::Type=Rational{Int}) where {Gr<:AbstractAlgebra.AbstractPermutationGroup, T<:GroupElem}
res = RG(S)
for g in V
res[g] += sign(g)
end
return res
end
function idempotents(RG::GroupRing{Generic.SymmetricGroup{S}}, T::Type=Rational{Int}) where S<:Integer
if RG.group.n == 1
return GroupRingElem{T}[one(RG,T)]
elseif RG.group.n == 2
Id = one(RG,T)
transp = RG(perm"(1,2)", T)
return GroupRingElem{T}[1//2*(Id + transp), 1//2*(Id - transp)]
end
projs = Vector{Vector{Generic.Perm{S}}}()
for l in 2:RG.group.n
u = RG.group([circshift([i for i in 1:l], -1); [i for i in l+1:RG.group.n]])
i = 0
while (l-1)*i <= RG.group.n
v = RG.group(circshift(collect(1:RG.group.n), i))
k = inv(v)*u*v
push!(projs, generateGroup([k], RG.group.n))
i += 1
end
end
idems = Vector{GroupRingElem{T}}()
for p in projs
append!(idems, [RG(p, T)//length(p), alternating_emb(RG, p, T)//length(p)])
end
return unique(idems)
end
function rankOne_projection(chi::PermCharacter, idems::Vector{T}) where {T<:GroupRingElem}
RG = parent(first(idems))
S = eltype(first(idems))
ids = [one(RG, S); idems]
zzz = zero(S)
for (i,j,k) in Base.product(ids, ids, ids)
if chi(i) == zzz || chi(j) == zzz || chi(k) == zzz
continue
else
elt = i*j*k
if elt^2 != elt
continue
elseif chi(elt) == one(S)
return elt
# return (i,j,k)
end
end
end
throw("Couldn't find rank-one projection for $chi")
end
function rankOne_projections(RG::GroupRing{<:Generic.SymmetricGroup}, T::Type=Rational{Int})
if RG.group.n == 1
return [GroupRingElem([one(T)], RG)]
end
RGidems = idempotents(RG, T)
min_projs = [central_projection(RG,chi)*rankOne_projection(chi,RGidems) for chi in characters(RG.group)]
return min_projs
end
function ifelsetuple(a,b, k, n)
x = [repeat([a], k); repeat([b], n-k)]
return tuple(x...)
end
function orbit_selector(n::Integer, k::Integer,
chi::AbstractCharacter, psi::AbstractCharacter)
return Projections.DirectProdCharacter(ifelsetuple(chi, psi, k, n))
end
function rankOne_projections(RBn::GroupRing{G}, T::Type=Rational{Int}) where {G<:WreathProduct}
Bn = RBn.group
N = Bn.P.n
# projections as elements of the group rings RSₙ
Sn_rankOnePr = [rankOne_projections(
GroupRing(SymmetricGroup(i), collect(SymmetricGroup(i))))
for i in typeof(N)(1):N]
# embedding into group ring of BN
RN = GroupRing(Bn.N, collect(Bn.N))
sign, id = collect(characters(Bn.N.group))
# Bn.N = (Z/2Z)ⁿ characters corresponding to the first k coordinates:
BnN_orbits = Dict(i => orbit_selector(N, i, sign, id) for i in 0:N)
Q = Dict(i => RBn(g -> Bn(g), central_projection(RN, BnN_orbits[i], T)) for i in 0:N)
Q = Dict(key => GroupRings.dense(val) for (key, val) in Q)
all_projs = [Q[0]*RBn(g->Bn(g), p) for p in Sn_rankOnePr[N]]
r = collect(1:N)
for i in 1:N-1
first_emb = g->Bn(Generic.emb!(one(Bn.P), g, view(r, 1:i)))
last_emb = g->Bn(Generic.emb!(one(Bn.P), g, view(r, (i+1):N)))
Sk_first = (RBn(first_emb, p) for p in Sn_rankOnePr[i])
Sk_last = (RBn(last_emb, p) for p in Sn_rankOnePr[N-i])
append!(all_projs,
[Q[i]*p1*p2 for (p1,p2) in Base.product(Sk_first,Sk_last)])
end
append!(all_projs, [Q[N]*RBn(g->Bn(g), p) for p in Sn_rankOnePr[N]])
return all_projs
end
##############################################################################
#
# General Groups Misc
#
##############################################################################
"""
products(X::Vector{GroupElem}, Y::Vector{GroupElem}[, op=*])
Return a vector of all possible products (or `op(x,y)`), where `x ∈ X` and
`y ∈ Y`. You may change which operation is used by specifying `op` argument.
"""
function products(X::AbstractVector{T}, Y::AbstractVector{T}, op=*) where {T<:GroupElem}
result = Vector{T}()
seen = Set{T}()
sizehint!(result, length(X)*length(Y))
sizehint!(seen, length(X)*length(Y))
for x in X
for y in Y
z = op(x,y)
if !in(z, seen)
push!(seen, z)
push!(result, z)
end
end
end
return result
end
"""
generateGroup(gens::Vector{GroupElem}, r, [op=*])
Produce all elements of a group generated by elements in `gens` in ball of
radius `r` (word-length metric induced by `gens`). The binary operation can
be optionally specified.
"""
function generateGroup(gens::Vector{T}, r, op=*) where {T<:GroupElem}
n = 0
R = 1
elts = [one(first(gens)); gens]
while n ≠ length(elts) && R < r
# @show elts
R += 1
n = length(elts)
elts = products(gens, elts, op)
end
return elts
end
end # of module Projections

29
src/actions/actions.jl Normal file
View File

@ -0,0 +1,29 @@
import SymbolicWedderburn.action
StarAlgebras.star(g::Groups.GroupElement) = inv(g)
include("alphabet_permutation.jl")
include("sln_conjugation.jl")
include("autfn_conjugation.jl")
function SymbolicWedderburn.action(
act::SymbolicWedderburn.ByPermutations,
g::Groups.GroupElement,
α::StarAlgebras.AlgebraElement
)
res = StarAlgebras.zero!(similar(α))
B = basis(parent(α))
for (idx, val) in StarAlgebras._nzpairs(StarAlgebras.coeffs(α))
a = B[idx]
a_g = SymbolicWedderburn.action(act, g, a)
res[a_g] += val
end
return res
end
function Base.:^(
w::W,
p::PermutationGroups.AbstractPerm,
) where {W<:Groups.AbstractWord}
return W([l^p for l in w])
end

44
src/actions/alphabet_permutation.jl Normal file
View File

@ -0,0 +1,44 @@
## action induced from permuting letters of an alphabet
import Groups: Constructions
struct AlphabetPermutation{GEl,I} <: SymbolicWedderburn.ByPermutations
perms::Dict{GEl,PermutationGroups.Perm{I}}
end
function AlphabetPermutation(
A::Alphabet,
Γ::PermutationGroups.AbstractPermutationGroup,
op,
)
return AlphabetPermutation(
Dict(γ => inv(PermutationGroups.Perm([A[op(l, γ)] for l in A])) for γ in Γ),
)
end
function AlphabetPermutation(A::Alphabet, W::Constructions.WreathProduct, op)
return AlphabetPermutation(
Dict(
w => inv(PermutationGroups.Perm([A[op(op(l, w.p), w.n)] for l in A])) for
w in W
),
)
end
function SymbolicWedderburn.action(
act::AlphabetPermutation,
γ::Groups.GroupElement,
g::Groups.AbstractFPGroupElement,
)
G = parent(g)
w = SymbolicWedderburn.action(act, γ, word(g))
return G(w)
end
function SymbolicWedderburn.action(
act::AlphabetPermutation,
γ::Groups.GroupElement,
w::Groups.AbstractWord,
)
return w^(act.perms[γ])
end

26
src/actions/autfn_conjugation.jl Normal file
View File

@ -0,0 +1,26 @@
## Particular definitions for actions on Aut(F_n)
function _conj(
t::Groups.Transvection,
σ::PermutationGroups.AbstractPerm,
)
return Groups.Transvection(t.id, t.i^inv(σ), t.j^inv(σ), t.inv)
end
function _flip(t::Groups.Transvection, g::Groups.GroupElement)
isone(g) && return t
return Groups.Transvection(t.id === :ϱ ? :λ : :ϱ, t.i, t.j, t.inv)
end
function _conj(
t::Groups.Transvection,
x::Groups.Constructions.DirectPowerElement,
)
@assert Groups.order(Int, parent(x).group) == 2
i, j = t.i, t.j
t = ifelse(isone(x.elts[i] * x.elts[j]), t, inv(t))
return _flip(t, x.elts[i])
end
action_by_conjugation(sautfn::Groups.AutomorphismGroup{<:Groups.FreeGroup}, Σ::Groups.Group) =
AlphabetPermutation(alphabet(sautfn), Σ, _conj)

20
src/actions/sln_conjugation.jl Normal file
View File

@ -0,0 +1,20 @@
## Particular definitions for actions on SL(n,ℤ)
function _conj(
t::MatrixGroups.ElementaryMatrix{N},
σ::PermutationGroups.AbstractPerm,
) where {N}
return MatrixGroups.ElementaryMatrix{N}(t.i^inv(σ), t.j^inv(σ), t.val)
end
function _conj(
t::MatrixGroups.ElementaryMatrix{N},
x::Groups.Constructions.DirectPowerElement,
) where {N}
@assert Groups.order(Int, parent(x).group) == 2
just_one_flips = xor(isone(x.elts[t.i]), isone(x.elts[t.j]))
return ifelse(just_one_flips, inv(t), t)
end
action_by_conjugation(sln::Groups.MatrixGroups.SpecialLinearGroup, Σ::Groups.Group) =
AlphabetPermutation(alphabet(sln), Σ, _conj)

27
src/actions/spn_conjugation.jl Normal file
View File

@ -0,0 +1,27 @@
## Particular definitions for actions on Sp(n,ℤ)
function _conj(
t::MatrixGroups.ElementarySymplectic{N,T},
σ::PermutationGroups.AbstractPerm,
) where {N,T}
@assert iseven(N)
@assert degree(σ) == N ÷ 2 "Got degree = $(degree(σ)); N = $N"
i = mod1(t.i, N ÷ 2)
ib = i == t.i ? 0 : N ÷ 2
j = mod1(t.j, N ÷ 2)
jb = j == t.j ? 0 : N ÷ 2
return MatrixGroups.ElementarySymplectic{N}(t.symbol, i^inv(σ) + ib, j^inv(σ) + jb, t.val)
end
function _conj(
t::MatrixGroups.ElementarySymplectic{N,T},
x::Groups.Constructions.DirectPowerElement,
) where {N,T}
@assert Groups.order(Int, parent(x).group) == 2
@assert iseven(N)
just_one_flips = xor(isone(x.elts[mod1(t.i, N ÷ 2)]), isone(x.elts[mod1(t.j, N ÷ 2)]))
return ifelse(just_one_flips, inv(t), t)
end
action_by_conjugation(sln::Groups.MatrixGroups.SymplecticGroup, Σ::Groups.Group) =
AlphabetPermutation(alphabet(sln), Σ, _conj)

287
src/blockdecomposition.jl
View File

@ -1,287 +0,0 @@
###############################################################################
#
# BlockDecomposition
#
###############################################################################
struct BlockDecomposition{T<:AbstractArray{Float64, 2}, GEl<:GroupElem, P<:Generic.Perm}
orbits::Vector{Vector{Int}}
preps::Dict{GEl, P}
Uπs::Vector{T}
dims::Vector{Int}
end
function BlockDecomposition(RG::GroupRing, autS::Group; verbose=true)
verbose && @info "Decomposing basis of RG into orbits of" autS
@time orbs = orbit_decomposition(autS, RG.basis, RG.basis_dict)
@assert sum(length(o) for o in orbs) == length(RG.basis)
verbose && @info "The action has $(length(orbs)) orbits"
verbose && @info "Finding projections in the Group Ring of" autS
@time autS_mps = Projections.rankOne_projections(GroupRing(autS, collect(autS)))
verbose && @info "Finding AutS-action matrix representation"
@time preps = perm_reps(autS, RG.basis[1:size(RG.pm,1)], RG.basis_dict)
@time mreps = matrix_reps(preps)
verbose && @info "Computing the projection matrices Uπs"
@time Uπs = [orthSVD(matrix_repr(p, mreps)) for p in autS_mps]
multiplicities = size.(Uπs,2)
dimensions = [Int(p[one(autS)]*Int(order(autS))) for p in autS_mps]
if verbose
info_strs = ["",
lpad("multiplicities", 14) * " =" * join(lpad.(multiplicities, 4), ""),
lpad("dimensions", 14) * " =" * join(lpad.(dimensions, 4), "")
]
@info join(info_strs, "\n")
end
@assert dot(multiplicities, dimensions) == size(RG.pm,1)
return BlockDecomposition(orbs, preps, Uπs, dimensions)
end
function decimate(od::BlockDecomposition; verbose=true)
nzros = [i for i in 1:length(od.Uπs) if !isempty(od.Uπs[i])]
Us = sparsify!.(od.Uπs, eps(Float64) * 1e4, verbose = verbose)[nzros]
#dimensions of the corresponding Uπs:
dims = od.dims[nzros]
return BlockDecomposition(od.orbits, od.preps, Array{Float64}.(Us), dims)
end
function orthSVD(M::AbstractMatrix{T}) where {T<:AbstractFloat}
fact = svd(convert(Matrix{T}, M))
M_rank = sum(fact.S .> maximum(size(M)) * eps(T))
return fact.U[:, 1:M_rank]
end
orbit_decomposition(
G::Group,
E::AbstractVector,
rdict = GroupRings.reverse_dict(E);
op = ^,
) = orbit_decomposition(collect(G), E, rdict; op=op)
function orbit_decomposition(elts::AbstractVector{<:GroupElem}, E::AbstractVector, rdict=GroupRings.reverse_dict(E); op=^)
tovisit = trues(size(E));
orbits = Vector{Vector{Int}}()
orbit = zeros(Int, length(elts))
for i in eachindex(E)
if tovisit[i]
g = E[i]
Threads.@threads for j in eachindex(elts)
orbit[j] = rdict[op(g, elts[j])]
end
tovisit[orbit] .= false
push!(orbits, unique(orbit))
end
end
return orbits
end
###############################################################################
#
# Sparsification
#
###############################################################################
dens(M::SparseMatrixCSC) = nnz(M)/length(M)
dens(M::AbstractArray) = count(!iszero, M)/length(M)
function sparsify!(M::SparseMatrixCSC{Tv,Ti}, tol=eps(Tv); verbose=false) where {Tv,Ti}
densM = dens(M)
droptol!(M, tol)
verbose && @info(
"Sparsified density:",
rpad(densM, 20),
" → ",
rpad(dens(M), 20),
" ($(nnz(M)) non-zeros)"
)
return M
end
function sparsify!(M::AbstractArray{T}, tol=eps(T); verbose=false) where T
densM = dens(M)
clamp_small!(M, tol)
if verbose
@info("Sparsifying $(size(M))-matrix... \n $(rpad(densM, 20)) → $(rpad(dens(M),20))), ($(count(!iszero, M)) non-zeros)")
end
return sparse(M)
end
function clamp_small!(M::AbstractArray{T}, tol=eps(T)) where T
for n in eachindex(M)
if abs(M[n]) < tol
M[n] = zero(T)
end
end
return M
end
function sparsify(U::AbstractArray{T}, tol=eps(T); verbose=false) where T
return sparsify!(deepcopy(U), tol, verbose=verbose)
end
###############################################################################
#
# perm-, matrix-, representations
#
###############################################################################
function perm_repr(g::GroupElem, E::Vector, E_dict)
p = Vector{Int}(undef, length(E))
for (i,elt) in enumerate(E)
p[i] = E_dict[elt^g]
end
return p
end
function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E))
elts = collect(G)
l = length(elts)
preps = Vector{Generic.Perm}(undef, l)
permG = SymmetricGroup(length(E))
Threads.@threads for i in 1:l
preps[i] = permG(PropertyT.perm_repr(elts[i], E, E_rdict), false)
end
return Dict(elts[i]=>preps[i] for i in 1:l)
end
function matrix_repr(x::GroupRingElem, mreps::Dict)
nzeros = findall(!iszero, x.coeffs)
return sum(x[i].*mreps[parent(x).basis[i]] for i in nzeros)
end
function matrix_reps(preps::Dict{T,Generic.Perm{I}}) where {T<:GroupElem, I<:Integer}
kk = collect(keys(preps))
mreps = Vector{SparseMatrixCSC{Float64, Int}}(undef, length(kk))
Threads.@threads for i in 1:length(kk)
mreps[i] = AbstractAlgebra.matrix_repr(preps[kk[i]])
end
return Dict(kk[i] => mreps[i] for i in 1:length(kk))
end
###############################################################################
#
# actions
#
###############################################################################
function Base.:^(y::GroupRingElem, g::GroupRingElem, op = ^)
res = parent(y)()
for elt in GroupRings.supp(g)
res += g[elt] * ^(y, elt, op)
end
return res
end
function Base.:^(y::GroupRingElem, g::GroupElem, op = ^)
RG = parent(y)
result = zero(RG, eltype(y.coeffs))
for (idx, c) in enumerate(y.coeffs)
if !iszero(c)
result[op(RG.basis[idx], g)] = c
end
end
return result
end
function Base.:^(
y::GroupRingElem{T,<:SparseVector},
g::GroupElem,
op = ^,
) where {T}
RG = parent(y)
index = [RG.basis_dict[op(RG.basis[idx], g)] for idx in y.coeffs.nzind]
result = GroupRingElem(sparsevec(index, y.coeffs.nzval, y.coeffs.n), RG)
return result
end
###############################################################################
#
# perm && WreathProductElems actions: MatAlgElem
#
###############################################################################
function Base.:^(A::MatAlgElem, p::Generic.Perm)
length(p.d) == size(A, 1) == size(A, 2) ||
throw("Can't act via $p on matrix of size $(size(A))")
result = similar(A)
@inbounds for i = 1:size(A, 1)
for j = 1:size(A, 2)
result[p[i], p[j]] = A[i, j] # action by permuting rows and colums/conjugation
end
end
return result
end
function Base.:^(A::MatAlgElem, g::WreathProductElem{N}) where {N}
# @assert N == size(A,1) == size(A,2)
flips = ntuple(i -> (g.n[i].d[1] == 1 && g.n[i].d[2] == 2 ? 1 : -1), N)
result = similar(A)
R = base_ring(parent(A))
tmp = R(1)
@inbounds for i = 1:size(A, 1)
for j = 1:size(A, 2)
x = A[i, j]
if flips[i] * flips[j] == 1
result[g.p[i], g.p[j]] = x
else
result[g.p[i], g.p[j]] = -x
end
end
end
return result
end
###############################################################################
#
# perm && WreathProductElems actions: Automorphism
#
###############################################################################
function Base.:^(a::Automorphism, g::GroupElem)
Ag = parent(a)(g)
Ag_inv = inv(Ag)
res = append!(Ag, a, Ag_inv)
return Groups.freereduce!(res)
end
(A::AutGroup)(p::Generic.Perm) = A(Groups.AutSymbol(p))
function (A::AutGroup)(g::WreathProductElem)
isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
parent(g).P.n == length(A.objectGroup.gens) ||
throw("No natural embedding of $(parent(g)) into $A")
elt = one(A)
Id = one(parent(g.n.elts[1]))
for i = 1:length(g.p.d)
if g.n.elts[i] != Id
push!(elt, Groups.flip(i))
end
end
push!(elt, Groups.AutSymbol(g.p))
return elt
end
# fallback:
Base.one(p::Generic.Perm) = Perm(length(p.d))

168
src/certify.jl Normal file
View File

@ -0,0 +1,168 @@
function augment_columns!(Q::AbstractMatrix)
for c in eachcol(Q)
c .-= sum(c) ./ length(c)
end
return Q
end
function _fma_SOS_thr!(
result::AbstractVector{T},
mstructure::AbstractMatrix{<:Integer},
Q::AbstractMatrix{T},
acc_matrix=zeros(T, size(mstructure)...),
) where {T}
s1, s2 = size(mstructure)
@inbounds for k = 1:s2
let k = k, s1 = s1, s2 = s2, Q = Q, acc_matrix = acc_matrix
Threads.@threads for j = 1:s2
for i = 1:s1
@inbounds acc_matrix[i, j] =
muladd(Q[i, k], Q[j, k], acc_matrix[i, j])
end
end
end
end
@inbounds for j = 1:s2
for i = 1:s1
result[mstructure[i, j]] += acc_matrix[i, j]
end
end
return result
end
function _cnstr_sos!(res::StarAlgebras.AlgebraElement, Q::AbstractMatrix, cnstrs)
StarAlgebras.zero!(res)
Q² = Q' * Q
for (g, A_g) in cnstrs
res[g] = dot(A_g, Q²)
end
return res
end
function _augmented_sos!(res::StarAlgebras.AlgebraElement, Q::AbstractMatrix)
A = parent(res)
StarAlgebras.zero!(res)
Q² = Q' * Q
N = LinearAlgebra.checksquare(A.mstructure)
augmented_basis = [A(1) - A(b) for b in @view basis(A)[1:N]]
tmp = zero(res)
for (j, y) in enumerate(augmented_basis)
for (i, x) in enumerate(augmented_basis)
# res += Q²[i, j] * x * y
StarAlgebras.mul!(tmp, x, y)
StarAlgebras.mul!(tmp, tmp, Q²[i, j])
StarAlgebras.add!(res, res, tmp)
end
end
return res
end
function compute_sos(A::StarAlgebras.StarAlgebra, Q::AbstractMatrix; augmented::Bool)
if augmented
z = zeros(eltype(Q), length(basis(A)))
res = StarAlgebras.AlgebraElement(z, A)
return _augmented_sos!(res, Q)
cnstrs = constraints(basis(A), A.mstructure; augmented=true)
return _cnstr_sos!(res, Q, cnstrs)
else
@assert size(A.mstructure) == size(Q)
z = zeros(eltype(Q), length(basis(A)))
_fma_SOS_thr!(z, A.mstructure, Q)
return StarAlgebras.AlgebraElement(z, A)
end
end
function sufficient_λ(residual::StarAlgebras.AlgebraElement, λ; halfradius)
L1_norm = norm(residual, 1)
suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm
eq_sign = let T = eltype(residual)
if T <: Interval
"∈"
elseif T <: Union{Rational,Integer}
"="
else # if T <: AbstractFloat
"≈"
end
end
info_strs = [
"Numerical metrics of the obtained SOS:",
"ɛ(elt - λu - ∑ξᵢ*ξᵢ) $eq_sign $(StarAlgebras.aug(residual))",
"‖elt - λu - ∑ξᵢ*ξᵢ‖₁ $eq_sign $(L1_norm)",
" λ $eq_sign $suff_λ",
]
@info join(info_strs, "\n")
return suff_λ
end
function sufficient_λ(
elt::StarAlgebras.AlgebraElement,
order_unit::StarAlgebras.AlgebraElement,
λ,
sos::StarAlgebras.AlgebraElement;
halfradius
)
@assert parent(elt) === parent(order_unit) == parent(sos)
residual = (elt - λ * order_unit) - sos
return sufficient_λ(residual, λ; halfradius=halfradius)
end
function certify_solution(
elt::StarAlgebras.AlgebraElement,
orderunit::StarAlgebras.AlgebraElement,
λ,
Q::AbstractMatrix{<:AbstractFloat};
halfradius,
augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit))
)
should_we_augment = !augmented && StarAlgebras.aug(elt) == StarAlgebras.aug(orderunit) == 0
Q = should_we_augment ? augment_columns!(Q) : Q
@time sos = compute_sos(parent(elt), Q, augmented=augmented)
@info "Checking in $(eltype(sos)) arithmetic with" λ
λ_flpoint = sufficient_λ(elt, orderunit, λ, sos, halfradius=halfradius)
if λ_flpoint ≤ 0
return false, λ_flpoint
end
λ_int = @interval(λ)
Q_int = [@interval(q) for q in Q]
check, sos_int = @time if should_we_augment
@info("Projecting columns of Q to the augmentation ideal...")
Q_int = augment_columns!(Q_int)
@info "Checking that sum of every column contains 0.0..."
check_augmented = all(0 ∈ sum(c) for c in eachcol(Q_int))
check_augmented || @error(
"Augmentation failed! The following numbers are not certified!"
)
sos_int = compute_sos(parent(elt), Q_int; augmented=augmented)
check_augmented, sos_int
else
true, compute_sos(parent(elt), Q_int, augmented=augmented)
end
@info "Checking in $(eltype(sos_int)) arithmetic with" λ
λ_certified =
sufficient_λ(elt, orderunit, λ_int, sos_int, halfradius=halfradius)
return check && inf(λ_certified) > 0.0, inf(λ_certified)
end

77
src/checksolution.jl
View File

@ -1,77 +0,0 @@
using IntervalArithmetic
IntervalArithmetic.setrounding(Interval, :tight)
IntervalArithmetic.setformat(sigfigs=12)
function fma_SOS_thr!(result::AbstractVector{T}, pm::AbstractMatrix{<:Integer},
Q::AbstractMatrix{T}, acc_matrix=zeros(T, size(pm)...)) where T
s1, s2 = size(pm)
@inbounds for k in 1:s2
let k=k, s1=s1, s2=s2, Q=Q, acc_matrix=acc_matrix
Threads.@threads for j in 1:s2
for i in 1:s1
@inbounds acc_matrix[i,j] = muladd(Q[i, k], Q[j, k], acc_matrix[i,j])
end
end
end
end
@inbounds for j in 1:s2
for i in 1:s1
result[pm[i,j]] += acc_matrix[i,j]
end
end
return result
end
function compute_SOS(pm::AbstractMatrix{<:Integer}, Q::AbstractMatrix)
result = zeros(eltype(Q), maximum(pm));
return fma_SOS_thr!(result, pm, Q)
end
function compute_SOS(RG::GroupRing, Q::AbstractMatrix{<:Real})
result = compute_SOS(RG.pm, Q)
return GroupRingElem(result, RG)
end
function compute_SOS_square(pm::AbstractMatrix{<:Integer}, Q::AbstractMatrix)
result = zeros(eltype(Q), maximum(pm));
for i in 1:size(Q,2)
GroupRings.fmac!(result, view(Q,:,i), view(Q,:,i), pm)
end
return result
end
function compute_SOS_square(RG::GroupRing, Q::AbstractMatrix)
return GroupRingElem(compute_SOS_square(RG.pm, Q), RG)
end
function augIdproj(Q::AbstractMatrix{T}) where T
result = zeros(T, size(Q))
l = size(Q, 2)
Threads.@threads for j in 1:l
col = sum(view(Q, :,j))/l
for i in 1:size(Q, 1)
result[i,j] = Q[i,j] - col
end
end
return result
end
function augIdproj(::Type{Interval}, Q::AbstractMatrix{T}) where {T<:Real}
result = zeros(Interval{T}, size(Q))
l = size(Q, 2)
Threads.@threads for j in 1:l
col = sum(view(Q, :,j))/l
for i in 1:size(Q, 1)
result[i,j] = @interval(Q[i,j] - col)
end
end
check = all([zero(T) in sum(view(result, :, i)) for i in 1:size(result, 2)])
return result, check
end

129
src/constraint_matrix.jl Normal file
View File

@ -0,0 +1,129 @@
"""
ConstraintMatrix{T,I} <: AbstractMatrix{T}
Special type of sparse matrix used to store constraints in SOS problems.
This matrix has in general very few non-zero values which also are multiples of each other.
The constructor accepts
* `nzeros` - a vector of non-zero indices; negative values are used to signify
negative values; repetitions are allowed
* `n`, `m` - the size of matrix
* `val` - the greatest common factor of the values
To iterate efficiently over `A::ConstraintMatrix` use [`nzpairs(A)`](@ref).
# Examples
```julia-repl
julia> ConstraintMatrix{Float64}([-1,2,-1,1,4,2,6], 3,2, π)
3×2 ConstraintMatrix{Float64, Int64}:
-3.14159 3.14159
6.28319 0.0
0.0 3.14159
```
"""
struct ConstraintMatrix{T,I} <: AbstractMatrix{T}
pos::Vector{I} # list of positive indices
neg::Vector{I} # list of negative indices
size::Tuple{Int,Int}
val::T
function ConstraintMatrix{T}(nzeros::AbstractArray{<:Integer}, n, m, val) where {T}
@assert n ≥ 1
@assert m ≥ 1
if !isempty(nzeros)
sort!(nzeros)
a, b = first(nzeros), last(nzeros)
@assert 1 ≤ abs(a) ≤ n * m
@assert 1 ≤ abs(b) ≤ n * m
end
k = searchsortedlast(nzeros, 0)
neg = @view nzeros[begin:k]
pos = @view nzeros[k+1:end]
return new{T,eltype(nzeros)}(pos, -neg, (n, m), val)
end
end
ConstraintMatrix(nzeros::AbstractArray{<:Integer}, n, m, val::T) where {T} =
ConstraintMatrix{T}(nzeros, n, m, val)
Base.size(cm::ConstraintMatrix) = cm.size
__get_positive(cm::ConstraintMatrix, idx::Integer) =
convert(eltype(cm), cm.val * length(searchsorted(cm.pos, idx)))
__get_negative(cm::ConstraintMatrix, idx::Integer) =
convert(eltype(cm), cm.val * length(searchsorted(cm.neg, idx)))
Base.@propagate_inbounds function Base.getindex(
cm::ConstraintMatrix,
i::Integer,
j::Integer,
)
li = LinearIndices(cm)
idx = li[i, j]
pos = __get_positive(cm, idx)
neg = __get_negative(cm, idx)
return pos - neg
end
struct NZPairsIter{T}
m::ConstraintMatrix{T}
end
Base.eltype(::Type{NZPairsIter{T}}) where {T} = Pair{Int,T}
Base.IteratorSize(::Type{<:NZPairsIter}) = Base.SizeUnknown()
# TODO: iterate over (idx=>val) pairs combining vals
function Base.iterate(itr::NZPairsIter, state::Tuple{Int,Int}=(1, 1))
k = iterate(itr.m.pos, state[1])
isnothing(k) && return iterate(itr, state[2])
idx, st = k
return idx => itr.m.val, (st, 1)
end
function Base.iterate(itr::NZPairsIter, state::Int)
k = iterate(itr.m.neg, state[1])
isnothing(k) && return nothing
idx, st = k
return idx => -itr.m.val, st
end
"""
nzpairs(cm::ConstraintMatrix)
Efficiently iterate over non-zero `(idx=>value)` pairs.
If the `cm` was created with repetitions (or contains negative values) there will
be repetitions in the returned sequence of pairs.
# Examples
```julia
julia> ConstraintMatrix{Float64}([-1,2,-1,1,4,2,6], 3,2, π)
3×2 ConstraintMatrix{Float64, Int64}:
-3.14159 3.14159
6.28319 0.0
0.0 3.14159
julia> collect(nzpairs(M))
7-element Vector{Pair{Int64, Float64}}:
1 => 3.141592653589793
2 => 3.141592653589793
2 => 3.141592653589793
4 => 3.141592653589793
6 => 3.141592653589793
1 => -3.141592653589793
1 => -3.141592653589793
```
"""
nzpairs(cm::ConstraintMatrix) = NZPairsIter(cm)
function LinearAlgebra.dot(cm::ConstraintMatrix, m::AbstractMatrix{T}) where {T}
if isempty(cm.pos) && isempty(cm.neg)
isempty(m) && return zero(T)
return zero(first(m) + first(m))
end
pos = isempty(cm.pos) ? zero(first(m)) : sum(@view m[cm.pos])
neg = isempty(cm.neg) ? zero(first(m)) : sum(@view m[cm.neg])
return convert(eltype(cm), cm.val) * (pos - neg)
end

69
src/gradings.jl Normal file
View File

@ -0,0 +1,69 @@
## something about roots
Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} =
Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j)
function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N}
if s.symbol === :A
return Roots.𝕖(N ÷ 2, s.i) - Roots.𝕖(N ÷ 2, s.j)
else#if s.symbol === :B
n = N ÷ 2
i, j = ifelse(s.i <= n, s.i, s.i - n), ifelse(s.j <= n, s.j, s.j - n)
return (-1)^(s.i > s.j) * (Roots.𝕖(n, i) + Roots.𝕖(n, j))
end
end
grading(s::MatrixGroups.ElementarySymplectic) = Roots.Root(s)
grading(e::MatrixGroups.ElementaryMatrix) = Roots.Root(e)
function grading(g::FPGroupElement)
if length(word(g)) == 1
A = alphabet(parent(g))
return grading(A[first(word(g))])
else
throw("Grading is implemented only for generators")
end
end
_groupby(f, iter::AbstractVector) = _groupby(f.(iter), iter)
function _groupby(keys::AbstractVector{K}, vals::AbstractVector{V}) where {K,V}
@assert length(keys) == length(vals)
d = Dict(k => V[] for k in keys)
for (k, v) in zip(keys, vals)
push!(d[k], v)
end
return d
end
function laplacians(RG::StarAlgebras.StarAlgebra, S, grading)
d = _groupby(grading, S)
Δs = Dict(α => RG(length(Sα)) - sum(RG(s) for s in Sα) for (α, Sα) in d)
return Δs
end
function Adj(rootsystem::AbstractDict, subtype::Symbol)
roots = let W = mapreduce(collect, union, keys(rootsystem))
W = union!(W, -1 .* W)
end
return reduce(
+,
(
Δα * Δβ for (α, Δα) in rootsystem for (β, Δβ) in rootsystem if
Roots.classify_sub_root_system(
roots,
first(α),
first(β),
) == subtype
),
init=zero(first(values(rootsystem))),
)
end
function level(rootsystem, level::Integer)
1 ≤ level ≤ 4 || throw("level is implemented only for i ∈{1,2,3,4}")
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 == 3 && return Adj(rootsystem, :A₂) + Adj(rootsystem, Symbol("A₁×C₁"))
level == 4 && return Adj(rootsystem, Symbol("A₁×A₁")) # positive
end

51
src/laplacians.jl
View File

@ -1,51 +0,0 @@
###############################################################################
#
# Laplacians
#
###############################################################################
function spLaplacian(RG::GroupRing, S::AbstractVector, T::Type=Float64)
result = RG(T)
result[one(RG.group)] = T(length(S))
for s in S
result[s] -= one(T)
end
return result
end
function Laplacian(S::AbstractVector{REl}, halfradius) where REl<:Union{NCRingElem, GroupElem}
G = parent(first(S))
@info "Generating metric ball of radius" radius=2halfradius
@time E_R, sizes = Groups.wlmetric_ball(S, radius=2halfradius)
@info "Generated balls:" sizes
@info "Creating product matrix..."
rdict = GroupRings.reverse_dict(E_R)
@time pm = GroupRings.create_pm(E_R, rdict, sizes[halfradius]; twisted=true)
RG = GroupRing(G, E_R, rdict, pm)
Δ = spLaplacian(RG, S)
return Δ
end
function saveGRElem(fname::String, g::GroupRingElem)
RG = parent(g)
JLD.save(fname, "coeffs", g.coeffs, "pm", RG.pm, "G", RG.group)
end
function loadGRElem(fname::String, RG::GroupRing)
coeffs = load(fname, "coeffs")
return GroupRingElem(coeffs, RG)
end
function loadGRElem(fname::String, G::Union{Group, NCRing})
pm = load(fname, "pm")
RG = GroupRing(G, pm)
return loadGRElem(fname, RG)
end
function loadGRElem(fname::String)
pm, G = load(fname, "pm", "G")
RG = GroupRing(G, pm)
return loadGRElem(fname, RG)
end

179
src/roots.jl Normal file
View File

@ -0,0 +1,179 @@
module Roots
using StaticArrays
using LinearAlgebra
export Root, isproportional, isorthogonal, ~, ⟂
abstract type AbstractRoot{N,T} end
struct Root{N,T} <: AbstractRoot{N,T}
coord::SVector{N,T}
end
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))
function isproportional(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
N == M || return false
val = abs(cos_angle(α, β))
return isapprox(val, one(val), atol=eps(one(val)))
end
function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
N == M || return false
val = cos_angle(α, β)
return isapprox(val, zero(val), atol=eps(one(val)))
end
function _positive_direction(α::Root{N}) where {N}
last = -1 / √2^(N - 1)
return Root{N,Float64}(
SVector(ntuple(i -> ifelse(i == N, last, (√2)^-i), N)),
)
end
function positive(roots::AbstractVector{<:Root{N}}) where {N}
# return those roots for which dot(α, Root([½, ¼, …])) > 0.0
pd = _positive_direction(first(roots))
return filter(α -> dot(α, pd) > 0.0, roots)
end
function Base.show(io::IO, r::Root)
print(io, "Root$(r.coord)")
end
function Base.show(io::IO, ::MIME"text/plain", r::Root{N}) where {N}
lngth² = sum(x -> x^2, r.coord)
l = isinteger(sqrt(lngth²)) ? "$(sqrt(lngth²))" : "√$(lngth²)"
print(io, "Root in ℝ^$N of length $l\n", r.coord)
end
𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
"""
classify_root_system(α, β)
Return the symbol of smallest system generated by roots `α` and `β`.
The classification is based only on roots length and
proportionality/orthogonality.
"""
function classify_root_system(α::AbstractRoot, β::AbstractRoot)
lα, lβ = length(α), length(β)
if isproportional(α, β)
if lα ≈ lβ ≈ √2
return :A₁
elseif lα ≈ lβ ≈ 2.0
return :C₁
else
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
end
elseif isorthogonal(α, β)
if lα ≈ lβ ≈ √2
return Symbol("A₁×A₁")
elseif lα ≈ lβ ≈ 2.0
return Symbol("C₁×C₁")
elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2)
return Symbol("A₁×C₁")
else
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
end
else # ⟨α, β⟩ is 2-dimensional, but they're not orthogonal
if lα ≈ lβ ≈ √2
return :A₂
elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2)
return :C₂
else
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
end
end
end
function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root)
k = findfirst(v -> isproportional(α, v), Ω)
if isnothing(k)
error("Line L_α not contained in root system Ω:\n α = $α\n Ω = $Ω")
end
return Ω[k]
end
struct Plane{R<:Root}
v1::R
v2::R
vectors::Vector{R}
end
Plane(α::R, β::R) where {R<:Root} =
Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
function Base.in(r::R, plane::Plane{R}) where {R}
return any(isproportional(r, v) for v in plane.vectors)
end
function classify_sub_root_system(
Ω::AbstractVector{<:Root{N}},
α::Root{N},
β::Root{N},
) where {N}
v = proportional_root_from_system(Ω, α)
w = proportional_root_from_system(Ω, β)
subsystem = filter(ω -> ω in Plane(v, w), Ω)
@assert length(subsystem) > 0
subsystem = positive(union(subsystem, -1 .* subsystem))
l = length(subsystem)
if l == 1
x = first(subsystem)
return classify_root_system(x, x)
elseif l == 2
return classify_root_system(subsystem...)
elseif l == 3
a = classify_root_system(subsystem[1], subsystem[2])
b = classify_root_system(subsystem[2], subsystem[3])
c = classify_root_system(subsystem[1], subsystem[3])
if a == b == c # it's only A₂
return a
end
C = (:C₂, Symbol("C₁×C₁"))
if (a ∈ C && b ∈ C && c ∈ C) && (:C₂ ∈ (a, b, c))
return :C₂
end
elseif l == 4
for i = 1:l
for j = (i+1):l
T = classify_root_system(subsystem[i], subsystem[j])
T == :C₂ && return :C₂
end
end
end
@error "Unknown root subsystem generated by" α β
throw("Unknown root system: $subsystem")
end
end # of module Roots

515
src/sos_sdps.jl
View File

@ -1,260 +1,349 @@
############################################################################### """
# sos_problem_dual(X, [u = zero(X); upper_bound=Inf])
# Constraints Formulate the dual to the sum of squares decomposition problem for `X - λ·u`.
#
###############################################################################
function constraints(pm::Matrix{I}, total_length=maximum(pm)) where {I<:Integer} See also [sos_problem_primal](@ref).
cnstrs = [Vector{I}() for _ in 1:total_length] """
for i in eachindex(pm) function sos_problem_dual(
push!(cnstrs[pm[i]], i) elt::StarAlgebras.AlgebraElement,
order_unit::StarAlgebras.AlgebraElement=zero(elt);
lower_bound=-Inf
)
@assert parent(elt) == parent(order_unit)
algebra = parent(elt)
mstructure = if StarAlgebras._istwisted(algebra.mstructure)
algebra.mstructure
else
StarAlgebras.MTable{true}(basis(algebra), table_size=size(algebra.mstructure))
end end
return cnstrs
# 1 variable for every primal constraint
# 1 dual variable for every element of basis
# Symmetrized:
# 1 dual variable for every orbit of G acting on basis
model = Model()
@variable model y[1:length(basis(algebra))]
@constraint model λ_dual dot(order_unit, y) == 1
@constraint(model, psd, y[mstructure] in PSDCone())
if !isinf(lower_bound)
throw("Not Implemented yet")
@variable model λ_ub_dual
@objective model Min dot(elt, y) + lower_bound * λ_ub_dual
else
@objective model Min dot(elt, y)
end
return model
end end
function orbit_constraint!(result::SparseMatrixCSC, cnstrs, orbit; val=1.0/length(orbit)) 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])
Formulate sum of squares decomposition problem for `X - λ·u`.
Given element `X` of a star-algebra `A` and an order unit `u` of `Σ²A` the sum
of squares cone in `A`, fomulate sum of squares decomposition problem:
```
max: λ
subject to: X - λ·u ∈ Σ²A
```
If `upper_bound` is given a finite value, additionally `λ ≤ upper_bound` will
be added to the model. This may improve the accuracy of the solution if
`upper_bound` is less than the optimal `λ`.
The default `u = zero(X)` formulates a simple feasibility problem.
"""
function sos_problem_primal(
elt::StarAlgebras.AlgebraElement,
order_unit::StarAlgebras.AlgebraElement=zero(elt);
upper_bound=Inf,
augmented::Bool=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(order_unit))
)
@assert parent(elt) === parent(order_unit)
N = LinearAlgebra.checksquare(parent(elt).mstructure)
model = JuMP.Model()
P = JuMP.@variable(model, P[1:N, 1:N], Symmetric)
JuMP.@constraint(model, psd, P in PSDCone())
if iszero(order_unit) && isfinite(upper_bound)
@warn "Setting `upper_bound` together with zero `order_unit` has no effect"
end
A = constraints(parent(elt), augmented=augmented, twisted=true)
if !iszero(order_unit)
λ = JuMP.@variable(model, λ)
if isfinite(upper_bound)
JuMP.@constraint model λ <= upper_bound
end
JuMP.@objective(model, Max, λ)
for b in basis(parent(elt))
JuMP.@constraint(model, elt(b) - λ * order_unit(b) == dot(A[b], P))
end
else
for b in basis(parent(elt))
JuMP.@constraint(model, elt(b) == dot(A[b], P))
end
end
return model
end
function invariant_constraint!(
result::AbstractMatrix,
basis::StarAlgebras.AbstractBasis,
cnstrs::AbstractDict{K,CM},
invariant_vec::SparseVector,
) where {K,CM<:ConstraintMatrix}
result .= zero(eltype(result)) result .= zero(eltype(result))
dropzeros!(result) for i in SparseArrays.nonzeroinds(invariant_vec)
for constraint in cnstrs[orbit] g = basis[i]
for idx in constraint A = cnstrs[g]
result[idx] = val for (idx, v) in nzpairs(A)
result[idx] += invariant_vec[i] * v
end end
end end
return result return result
end end
function orbit_spvector(vect::AbstractVector, orbits) function isorth_projection(ds::SymbolicWedderburn.DirectSummand)
orb_vector = spzeros(length(orbits)) U = SymbolicWedderburn.image_basis(ds)
return isapprox(U * U', I)
for (i,o) in enumerate(orbits)
k = vect[collect(o)]
val = k[1]
@assert all(k .== val)
orb_vector[i] = val
end
return orb_vector
end end
############################################################################### sos_problem_primal(
# elt::StarAlgebras.AlgebraElement,
# Naive SDP wedderburn::WedderburnDecomposition;
# kwargs...
############################################################################### ) = sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)
function SOS_problem_dual(elt::GroupRingElem, order_unit::GroupRingElem; function sos_problem_primal(
lower_bound::Float64=Inf) elt::StarAlgebras.AlgebraElement,
@assert parent(elt) == parent(order_unit) orderunit::StarAlgebras.AlgebraElement,
wedderburn::WedderburnDecomposition;
upper_bound=Inf,
augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)),
check_orthogonality=true
)
RG = parent(elt) @assert parent(elt) === parent(orderunit)
m = Model() if check_orthogonality
if any(!isorth_projection, direct_summands(wedderburn))
y = @variable(m, y[1:length(elt.coeffs)]) error("Wedderburn decomposition contains a non-orthogonal projection")
@constraint(m, λ_dual, dot(order_unit.coeffs, y) == 1)
@constraint(m, psd, [y[i] for i in RG.pm] in PSDCone())
if !isinf(lower_bound)
@variable(m, λ_ub_dual)
expr = dot(elt.coeffs, y) + lower_bound*λ_ub_dual
# @constraint m expr >= lower_bound
@objective m Min expr
else
@objective m Min dot(elt.coeffs, y)
end
return m
end
function SOS_problem_primal(X::GroupRingElem, orderunit::GroupRingElem;
upper_bound::Float64=Inf)
N = size(parent(X).pm, 1)
m = JuMP.Model();
JuMP.@variable(m, P[1:N, 1:N])
# SP = Symmetric(P)
JuMP.@SDconstraint(m, sdp, P >= 0)
if iszero(aug(X)) && iszero(aug(orderunit))
JuMP.@constraint(m, augmentation, sum(P) == 0)
end
if upper_bound < Inf
λ = JuMP.@variable(m, λ <= upper_bound)
else
λ = JuMP.@variable(m, λ)
end
cnstrs = constraints(parent(X).pm)
@assert length(cnstrs) == length(X.coeffs) == length(orderunit.coeffs)
x, u = X.coeffs, orderunit.coeffs
JuMP.@constraint(m, lincnstr[i=1:length(cnstrs)],
x[i] - λ*u[i] == sum(P[cnstrs[i]]))
JuMP.@objective(m, Max, λ)
return m
end
###############################################################################
#
# Symmetrized SDP
#
###############################################################################
function SOS_problem_primal(X::GroupRingElem, orderunit::GroupRingElem, data::BlockDecomposition; upper_bound::Float64=Inf)
Ns = size.(data.Uπs, 2)
m = JuMP.Model();
Ps = Vector{Matrix{JuMP.VariableRef}}(undef, length(Ns))
for (k,s) in enumerate(Ns)
Ps[k] = JuMP.@variable(m, [1:s, 1:s])
JuMP.@SDconstraint(m, Ps[k] >= 0)
end
if upper_bound < Inf
λ = JuMP.@variable(m, λ <= upper_bound)
else
λ = JuMP.@variable(m, λ)
end
@info "Adding $(length(data.orbits)) constraints..."
@time addconstraints!(m, Ps, X, orderunit, data)
JuMP.@objective(m, Max, λ)
return m, Ps
end
function constraintLHS!(M, cnstr, Us, Ust, dims, eps=1000*eps(eltype(first(M))))
for π in 1:lastindex(Us)
M[π] = dims[π].*PropertyT.clamp_small!(Ust[π]*(cnstr*Us[π]), eps)
end
end
function addconstraints!(m::JuMP.Model,
P::Vector{Matrix{JuMP.VariableRef}},
X::GroupRingElem, orderunit::GroupRingElem, data::BlockDecomposition)
orderunit_orb = orbit_spvector(orderunit.coeffs, data.orbits)
X_orb = orbit_spvector(X.coeffs, data.orbits)
UπsT = [U' for U in data.Uπs]
cnstrs = constraints(parent(X).pm)
orb_cnstr = spzeros(Float64, size(parent(X).pm)...)
M = [Array{Float64}(undef, n,n) for n in size.(UπsT,1)]
λ = m[:λ]
for (t, orbit) in enumerate(data.orbits)
orbit_constraint!(orb_cnstr, cnstrs, orbit)
constraintLHS!(M, orb_cnstr, data.Uπs, UπsT, data.dims)
x, u = X_orb[t], orderunit_orb[t]
JuMP.@constraints m begin
x - λ*u == sum(dot(M[π], P[π]) for π in eachindex(data.Uπs))
end end
end end
return m
end
function reconstruct(Ps::Vector{Matrix{F}}, data::BlockDecomposition) where F feasibility_problem = iszero(orderunit)
return reconstruct(Ps, data.preps, data.Uπs, data.dims)
end
function reconstruct(Ps::Vector{M}, model = JuMP.Model()
preps::Dict{GEl, P}, Uπs::Vector{U}, dims::Vector{Int}) where if !feasibility_problem # add λ or not?
{M<:AbstractMatrix, GEl<:GroupElem, P<:Generic.Perm, U<:AbstractMatrix} λ = JuMP.@variable(model, λ)
JuMP.@objective(model, Max, λ)
lU = length(Uπs) if isfinite(upper_bound)
transfP = [dims[π].*Uπs[π]*Ps[π]*Uπs[π]' for π in 1:lU] JuMP.@constraint(model, λ <= upper_bound)
tmp = [zeros(Float64, size(first(transfP))) for _ in 1:lU] if feasibility_problem
@warn "setting `upper_bound` with zero `orderunit` has no effect"
Threads.@threads for π in 1:lU
tmp[π] = perm_avg!(tmp[π], transfP[π], values(preps))
end
recP = sum(tmp)./length(preps)
return recP
end
function perm_avg!(result, P, perms)
lp = length(first(perms).d)
for p in perms
# result .+= view(P, p.d, p.d)
@inbounds for j in 1:lp
k = p[j]
for i in 1:lp
result[i,j] += P[p[i], k]
end end
end end
end end
return result
end
############################################################################### P = map(direct_summands(wedderburn)) do ds
# dim = size(ds, 1)
# Low-level solve P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric)
# @constraint(model, P in PSDCone())
############################################################################### P
function setwarmstart!(m::JuMP.Model, warmstart)
if solver_name(m) == "SCS"
primal, dual, slack = warmstart
m.moi_backend.optimizer.model.optimizer.data.primal = primal
m.moi_backend.optimizer.model.optimizer.data.dual = dual
m.moi_backend.optimizer.model.optimizer.data.slack = slack
else
@warn "Setting warmstart for $(solver_name(m)) is not implemented! Ignoring..."
end end
return m
end
function getwarmstart(m::JuMP.Model) begin # preallocating
if solver_name(m) == "SCS" T = eltype(wedderburn)
warmstart = ( M = zeros.(T, size.(P))
primal = m.moi_backend.optimizer.model.optimizer.data.primal, M_orb = zeros(T, size(parent(elt).mstructure)...)
dual = m.moi_backend.optimizer.model.optimizer.data.dual, tmps = SymbolicWedderburn._tmps(wedderburn)
slack = m.moi_backend.optimizer.model.optimizer.data.slack
)
else
@warn "Saving warmstart for $(solver_name(m)) is not implemented!"
return (primal=Float64[], dual=Float64[], slack=Float64[])
end end
return warmstart
X = convert(Vector{T}, StarAlgebras.coeffs(elt))
U = convert(Vector{T}, StarAlgebras.coeffs(orderunit))
# defining constraints based on the multiplicative structure
cnstrs = constraints(parent(elt), augmented=augmented, twisted=true)
@info "Adding $(length(invariant_vectors(wedderburn))) constraints"
for iv in invariant_vectors(wedderburn)
x = dot(X, iv)
u = dot(U, iv)
M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv)
M = SymbolicWedderburn.diagonalize!(M, M_orb, wedderburn, tmps)
SymbolicWedderburn.zerotol!.(M, atol=1e-12)
M_dot_P = sum(dot(M[π], P[π]) for π in eachindex(M) if !iszero(M[π]))
if feasibility_problem
JuMP.@constraint(model, x == M_dot_P)
else
JuMP.@constraint(model, x - λ * u == M_dot_P)
end
end
return model, P
end end
function solve(m::JuMP.Model, with_optimizer::JuMP.OptimizerFactory, warmstart=nothing) function reconstruct(Ps, wd::WedderburnDecomposition)
N = size(first(direct_summands(wd)), 2)
P = zeros(eltype(wd), N, N)
return reconstruct!(P, Ps, wd)
end
set_optimizer(m, with_optimizer) function group_of(wd::WedderburnDecomposition)
# this is veeeery hacky... ;)
return parent(first(keys(wd.hom.cache)))
end
# TODO: move to SymbolicWedderburn
SymbolicWedderburn.action(wd::WedderburnDecomposition) =
SymbolicWedderburn.action(wd.hom)
function reconstruct!(
res::AbstractMatrix,
Ps,
wedderburn::WedderburnDecomposition,
)
G = group_of(wedderburn)
act = SymbolicWedderburn.action(wedderburn)
@assert act isa SymbolicWedderburn.ByPermutations
for (π, ds) in pairs(direct_summands(wedderburn))
Uπ = SymbolicWedderburn.image_basis(ds)
# LinearAlgebra.mul!(tmp, Uπ', P[π])
# LinearAlgebra.mul!(tmp2, tmp, Uπ)
tmp2 = Uπ' * Ps[π] * Uπ
if eltype(res) <: AbstractFloat
SymbolicWedderburn.zerotol!(tmp2, atol=1e-12)
end
tmp2 .*= SymbolicWedderburn.degree(ds)
@assert size(tmp2) == size(res)
for g in G
p = SymbolicWedderburn.induce(wedderburn.hom, g)
for c in axes(res, 2)
for r in axes(res, 1)
res[r, c] += tmp2[r^p, c^p]
end
end
end
end
res ./= Groups.order(Int, G)
return res
end
##
# Low-level solve
setwarmstart!(model::JuMP.Model, ::Nothing) = model
function setwarmstart!(model::JuMP.Model, warmstart)
constraint_map = Dict(
ct => JuMP.all_constraints(model, ct...) for
ct in JuMP.list_of_constraint_types(model)
)
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])
JuMP.set_dual_start_value.(idx, warmstart.dual[ct])
end
return model
end
function getwarmstart(model::JuMP.Model)
constraint_map = Dict(
ct => JuMP.all_constraints(model, ct...) for
ct in JuMP.list_of_constraint_types(model)
)
primal = value.(JuMP.all_variables(model))
slack = Dict(k => value.(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)
end
function solve(m::JuMP.Model, optimizer, warmstart=nothing)
JuMP.set_optimizer(m, optimizer)
MOIU.attach_optimizer(m) MOIU.attach_optimizer(m)
if warmstart != nothing m = setwarmstart!(m, warmstart)
setwarmstart!(m, warmstart)
end
optimize!(m) JuMP.optimize!(m)
Base.Libc.flush_cstdio() Base.Libc.flush_cstdio()
status = termination_status(m) status = JuMP.termination_status(m)
return status, getwarmstart(m) return status, getwarmstart(m)
end end
function solve(solverlog::String, m::JuMP.Model, with_optimizer::JuMP.OptimizerFactory, 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)
Base.Libc.flush_cstdio()
status, warmstart = open(solverlog, "a+") do logfile status, warmstart = open(solverlog, "a+") do logfile
redirect_stdout(logfile) do redirect_stdout(logfile) do
status, warmstart = PropertyT.solve(m, with_optimizer, warmstart) status, warmstart = solve(m, optimizer, warmstart)
status, warmstart status, warmstart
end end
end end

134
src/sqadjop.jl
View File

@ -1,98 +1,80 @@
isopposite(σ::Generic.Perm, τ::Generic.Perm, i=1, j=2) = import SymbolicWedderburn.PermutationGroups.AbstractPerm
σ[i] ≠ τ[i] && σ[i] ≠ τ[j] &&
σ[j] ≠ τ[i] && σ[j] ≠ τ[j]
isadjacent(σ::Generic.Perm, τ::Generic.Perm, i=1, j=2) = # move to Groups
(σ[i] == τ[i] && σ[j] ≠ τ[j]) || # first equal, second differ Base.keys(a::Alphabet) = keys(1:length(a))
(σ[j] == τ[j] && σ[i] ≠ τ[i]) || # second equal, first differ
(σ[i] == τ[j] && σ[j] ≠ τ[i]) || # first σ equal to second τ
(σ[j] == τ[i] && σ[i] ≠ τ[j]) # second σ equal to first τ
Base.div(X::GroupRingElem, x::Number) = parent(X)(X.coeffs.÷x) ## the old 1812.03456 definitions
function Sq(RG::GroupRing, N::Integer) isopposite(σ::AbstractPerm, τ::AbstractPerm, i=1, j=2) =
S₂ = generating_set(RG.group, 2) i^σ ≠ i^τ && i^σ ≠ j^τ && j^σ ≠ i^τ && j^σ ≠ j^τ
ℤ = Int64
Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ);
Alt_N = [g for g in SymmetricGroup(N) if parity(g) == 0] isadjacent(σ::AbstractPerm, τ::AbstractPerm, i=1, j=2) =
(i^σ == i^τ && j^σ ≠ j^τ) || # first equal, second differ
(j^σ == j^τ && i^σ ≠ i^τ) || # second equal, first differ
(i^σ == j^τ && j^σ ≠ i^τ) || # first σ equal to second τ
(j^σ == i^τ && i^σ ≠ j^τ) # second σ equal to first τ
sq = RG() function _ncycle(start, length, n=start + length - 1)
for σ in Alt_N p = PermutationGroups.Perm(Int8(n))
GroupRings.addeq!(sq, *(Δ₂^σ, Δ₂^σ, false)) @assert n ≥ start + length - 1
for k in start:start+length-2
p[k] = k + 1
end end
return sq÷factorial(N-2) p[start+length-1] = start
return p
end end
function Adj(RG::GroupRing, N::Integer) alternating_group(n::Integer) = PermutationGroups.PermGroup([_ncycle(i, 3) for i in 1:n-2])
S₂ = generating_set(RG.group, 2)
ℤ = Int64
Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ);
Alt_N = [g for g in SymmetricGroup(N) if parity(g) == 0] function small_gens(G::MatrixGroups.SpecialLinearGroup)
Δ₂s = Dict(σ=>Δ₂^σ for σ in Alt_N) A = alphabet(G)
adj = RG() S = map([(1, 2), (2, 1)]) do (i, j)
k = findfirst(l -> (l.i == i && l.j == j), A)
for σ in Alt_N @assert !isnothing(k)
for τ in Alt_N G(Groups.word_type(G)([k]))
if isadjacent(σ, τ)
GroupRings.addeq!(adj, *(Δ₂s[σ], Δ₂s[τ], false))
end
end
end end
return adj÷factorial(N-2)^2 return union!(S, inv.(S))
end end
function Op(RG::GroupRing, N::Integer) function small_gens(G::Groups.AutomorphismGroup{<:FreeGroup})
if N < 4 A = alphabet(G)
return RG() S = map([(:ϱ, 1, 2), (:ϱ, 2, 1), (:λ, 1, 2), (:λ, 2, 1)]) do (id, i, j)
k = findfirst(l -> (l.id == id && l.i == i && l.j == j), A)
@assert !isnothing(k)
G(Groups.word_type(G)([k]))
end end
S₂ = generating_set(RG.group, 2) return union!(S, inv.(S))
ℤ = Int64 end
Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ);
Alt_N = [g for g in SymmetricGroup(N) if parity(g) == 0] function small_laplacian(RG::StarAlgebras.StarAlgebra)
Δ₂s = Dict(σ=>Δ₂^σ for σ in Alt_N) G = StarAlgebras.object(RG)
op = RG() S₂ = small_gens(G)
return length(S₂) * one(RG) - sum(RG(s) for s in S₂)
end
for σ in Alt_N function SqAdjOp(A::StarAlgebras.StarAlgebra, n::Integer, Δ₂=small_laplacian(A))
for τ in Alt_N @assert parent(Δ₂) === A
alt_n = alternating_group(n)
G = StarAlgebras.object(A)
act = action_by_conjugation(G, alt_n)
Δ₂s = Dict(σ => SymbolicWedderburn.action(act, σ, Δ₂) for σ in alt_n)
sq, adj, op = zero(A), zero(A), zero(A)
tmp = zero(A)
for σ in alt_n
StarAlgebras.add!(sq, sq, StarAlgebras.mul!(tmp, Δ₂s[σ], Δ₂s[σ]))
for τ in alt_n
if isopposite(σ, τ) if isopposite(σ, τ)
GroupRings.addeq!(op, *(Δ₂s[σ], Δ₂s[τ], false)) StarAlgebras.add!(op, op, StarAlgebras.mul!(tmp, Δ₂s[σ], Δ₂s[τ]))
end
end
end
return op÷factorial(N-2)^2
end
for Elt in [:Sq, :Adj, :Op]
@eval begin
$Elt(RG::GroupRing{AutGroup{N}}) where N = $Elt(RG, N)
$Elt(RG::GroupRing{<:MatAlgebra}) = $Elt(RG, nrows(RG.group))
end
end
function SqAdjOp(RG::GroupRing, N::Integer)
S₂ = generating_set(RG.group, 2)
ℤ = Int64
Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ);
Alt_N = [σ for σ in SymmetricGroup(N) if parity(σ) == 0]
sq, adj, op = RG(), RG(), RG()
Δ₂s = Dict(σ=>Δ₂^σ for σ in Alt_N)
for σ in Alt_N
GroupRings.addeq!(sq, *(Δ₂s[σ], Δ₂s[σ], false))
for τ in Alt_N
if isopposite(σ, τ)
GroupRings.addeq!(op, *(Δ₂s[σ], Δ₂s[τ], false))
elseif isadjacent(σ, τ) elseif isadjacent(σ, τ)
GroupRings.addeq!(adj, *(Δ₂s[σ], Δ₂s[τ], false)) StarAlgebras.add!(adj, adj, StarAlgebras.mul!(tmp, Δ₂s[σ], Δ₂s[τ]))
end end
end end
end end
k = factorial(N-2) k = factorial(n - 2)
return sq÷k, adj÷k^2, op÷k^2 return sq ÷ k, adj ÷ k^2, op ÷ k^2
end end

243
test/1703.09680.jl
View File

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

322
test/1712.07167.jl
View File

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

364
test/1812.03456.jl
View File

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

261
test/actions.jl
View File

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

21
test/constratint_matrices.jl Normal file
View File

@ -0,0 +1,21 @@
@testset "ConstraintMatrix" begin
@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, π)
@test cm == Float64[
-π π
2π 0.0
0.0 π
]
@test collect(PropertyT.nzpairs(cm)) == [
1 => 3.141592653589793
2 => 3.141592653589793
2 => 3.141592653589793
4 => 3.141592653589793
6 => 3.141592653589793
1 => -3.141592653589793
1 => -3.141592653589793
]
end

86
test/graded_adj.jl Normal file
View File

@ -0,0 +1,86 @@
@testset "Adj via grading" begin
@testset "SL(n,Z) & Aut(F₄)" begin
n = 4
halfradius = 1
SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=halfradius, twisted=true)
Δ = RSL(length(S)) - sum(RSL(s) for s in S)
Δs = let ψ = identity
PropertyT.laplacians(
RSL,
S,
x -> (gx = PropertyT.grading(ψ(x)); Set([gx, -gx])),
)
end
sq, adj, op = PropertyT.SqAdjOp(RSL, n)
@test PropertyT.Adj(Δs, :A₁) == sq
@test PropertyT.Adj(Δs, :A₂) == adj
@test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
halfradius = 1
G = SpecialAutomorphismGroup(FreeGroup(n))
RG, S, sizes = PropertyT.group_algebra(G, halfradius=halfradius, twisted=true)
Δ = RG(length(S)) - sum(RG(s) for s in S)
Δs = let ψ = Groups.Homomorphism(Groups._abelianize, G, SL)
PropertyT.laplacians(
RG,
S,
x -> (gx = PropertyT.grading(ψ(x)); Set([gx, -gx])),
)
end
sq, adj, op = PropertyT.SqAdjOp(RG, n)
@test PropertyT.Adj(Δs, :A₁) == sq
@test PropertyT.Adj(Δs, :A₂) == adj
@test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
end
@testset "Symplectic group" begin
genus = 3
halfradius = 1
SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true)
Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
Δ = RG(length(S)) - sum(RG(s) for s in S)
Δs = PropertyT.laplacians(
RG,
S,
x -> (gx = PropertyT.grading(ψ(x)); Set([gx, -gx])),
)
Δ, Δs
end
@testset "Adj correctness: genus=$genus" begin
all_subtypes = (
:A₁, :C₁, Symbol("A₁×A₁"), Symbol("C₁×C₁"), Symbol("A₁×C₁"), :A₂, :C₂
)
@test PropertyT.Adj(Δs, :A₂)[one(SpN)] == 384
@test iszero(PropertyT.Adj(Δs, Symbol("A₁×A₁")))
@test iszero(PropertyT.Adj(Δs, Symbol("C₁×C₁")))
@testset "divisibility by 16" begin
for subtype in all_subtypes
subtype in (:A₁, :C₁) && continue
@test isinteger(PropertyT.Adj(Δs, subtype)[one(SpN)] / 16)
end
end
@test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) == Δ^2
end
end
end

54
test/optimizers.jl Normal file
View File

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

36
test/runtests.jl
View File

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