From 9511e34de4db6388765884a98edd6c93e1f59873 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 15:34:30 +0100 Subject: [PATCH] replace BlockDecomposition by SymbolicWedderburn --- Project.toml | 9 +- src/PropertyT.jl | 8 +- src/RGprojections.jl | 251 --------------------------------- src/blockdecomposition.jl | 287 -------------------------------------- 4 files changed, 2 insertions(+), 553 deletions(-) delete mode 100644 src/RGprojections.jl delete mode 100644 src/blockdecomposition.jl diff --git a/Project.toml b/Project.toml index 40d907f..155508e 100644 --- a/Project.toml +++ b/Project.toml @@ -4,23 +4,16 @@ authors = ["Marek Kaluba "] 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" diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 2ea9c30..8cba26d 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -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") diff --git a/src/RGprojections.jl b/src/RGprojections.jl deleted file mode 100644 index 50ceb86..0000000 --- a/src/RGprojections.jl +++ /dev/null @@ -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 diff --git a/src/blockdecomposition.jl b/src/blockdecomposition.jl deleted file mode 100644 index d4af161..0000000 --- a/src/blockdecomposition.jl +++ /dev/null @@ -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))