replace BlockDecomposition by SymbolicWedderburn

This commit is contained in:
Marek Kaluba 2022-11-07 15:34:30 +01:00
parent 92d9e468d2
commit 9511e34de4
No known key found for this signature in database
GPG Key ID: 8BF1A3855328FC15
4 changed files with 2 additions and 553 deletions

View File

@ -4,23 +4,16 @@ authors = ["Marek Kaluba <kalmar@amu.edu.pl>"]
version = "0.3.2"
[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"
IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
[compat]
AbstractAlgebra = "^0.10.0"
GroupRings = "^0.3.2"
Groups = "^0.5.0"
IntervalArithmetic = "^0.16.0"
JLD = "^0.9.0"
JuMP = "^0.20.0"
SCS = "^0.7.0"
julia = "^1.3.0"

View File

@ -1,22 +1,16 @@
__precompile__()
module PropertyT
using AbstractAlgebra
using LinearAlgebra
using SparseArrays
using Dates
using Groups
using GroupRings
using SymbolicWedderburn
using JLD
using JuMP
import AbstractAlgebra: Group, NCRing
include("laplacians.jl")
include("RGprojections.jl")
include("blockdecomposition.jl")
include("sos_sdps.jl")
include("checksolution.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

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))