central projections

[kalmarek]

@tweisser just in case you're bored over the next week:

function central_projection(chi::SymbolicWedderburn.AbstractClassFunction{T}) where T
    ccls = conjugacy_classes(chi)
    d = degree(one(first(first(ccls))))

    result = zeros(T, d, d)

    for cc in ccls
        val = chi(first(cc))
        for g in cc
            for i in 1:d
                result[i, i^g] += val
            end
        end
    end
    deg = degree(chi)
    ordG = sum(length, ccls)

    return deg//ordG, result
end

G = PermGroup([perm"(2,3)(4,5)"])
ccG = conjugacy_classes(G)

chars = SymbolicWedderburn.characters_dixon(ccG)

then

julia> χ = chars[1]
Character over Cyclotomics.Cyclotomic{Int64,SparseArrays.SparseVector{Int64,Int64}}
()^G         →   +1*E(1)^0
(2,3)(4,5)^G →   +1*E(1)^0


julia> c, U = SymbolicWedderburn.central_projection(χ);

julia> b, _ = SymbolicWedderburn.row_echelon_form(1.0*U);

julia> b
5×5 Array{Cyclotomics.Cyclotomic{Float64,SparseArrays.SparseVector{Float64,Int64}},2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  1.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> χ = chars[2]
Character over Cyclotomics.Cyclotomic{Int64,SparseArrays.SparseVector{Int64,Int64}}
()^G         →   +1*E(1)^0
(2,3)(4,5)^G →   -1*E(1)^0


julia> c, U = SymbolicWedderburn.central_projection(χ);

julia> b, _ = SymbolicWedderburn.row_echelon_form(1.0*U);

julia> b
5×5 Array{Cyclotomics.Cyclotomic{Float64,SparseArrays.SparseVector{Float64,Int64}},2}:
 0.0  1.0  -1.0  0.0   0.0
 0.0  0.0   0.0  1.0  -1.0
 0.0  0.0   0.0  0.0   0.0
 0.0  0.0   0.0  0.0   0.0
 0.0  0.0   0.0  0.0   0.0

which is the invariant basis from the example.

[tweisser]

I'm slightly confused. Where did the 1 for the constant polynomial go? Also shouldn't we get x[i]+x[i+1] and x[i]-x[i+1] for i=1,3?

When I want to minimize x+y over 1-x^2-y^2>=0, the symmetric SOS formulation should look like

max t
x+y-t = [1 x+y x-y] A [1 x+y x-y]  + a * (1-x^2-y^2)
A psd, a >=0 

[kalmarek]

it's in the first row of the first b?

[tweisser]

Yes, but where did it go in the second b?

[kalmarek]

it's an orthogonal decompostion into invariant spaces, it can't be in both of them at the same time...

Your A will consist now of 2×2 block (corresponding to the trivial repr) and 1×1 (sign repr)

[tweisser]

We will need to transform the matrices into complex numbers as the matrices we are looking for are

R1, R2+R3, and (R2-R3)/(sqrt(-1))

in julia we get

julia>R1
4×4 Array{Cyclotomic{Float64,SparseArrays.SparseVector{Float64,Int64}},2}:
 1.0  1.0  1.0  0.0
 0.0  0.0  0.0  1.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> R2+R3
4×4 Array{Cyclotomic{Float64,SparseArrays.SparseVector{Float64,Int64}},2}:
 2.0  -1.0  -1.0  0.0
 0.0   0.0   0.0  0.0
 0.0   0.0   0.0  0.0
 0.0   0.0   0.0  0.0

julia>  1/im.*convert.(ComplexF64, R2-R3)
4×4 Array{Complex{Float64},2}:
 0.0+0.0im  1.73205-6.66134e-16im  -1.73205+6.66134e-16im  0.0+0.0im
 0.0+0.0im      0.0+0.0im               0.0+0.0im          0.0+0.0im
 0.0+0.0im      0.0+0.0im               0.0+0.0im          0.0+0.0im
 0.0+0.0im      0.0+0.0im               0.0+0.0im          0.0+0.0im

So using Julia's complex numbers might be a bad decision. However, we know that we will only substract complex conjugate pairs (and divide by the imaginary unit).

[kalmarek]

Inducing Actions

I came up with this design for inducing actions:

abstract type AbstractAction{T} end

struct OnPoints{T} <: AbstractAction{T}
    features::Vector{T}
    reversef::Dict{T, Int}
end

Base.getindex(act::OnPoints, i::Integer) = act.features[i]
Base.getindex(act::OnPoints{T}, f::T) where T = act.reversef[f]
(act::OnPoints)(p::Perm{I}) where I = Perm(I[act[f^p] for f in act.features])

function OnPoints(basis::Union{MonomialVector, AbstractVector{<:Monomial}})
    # action on monomials = action on exponent vectors in this case
    basis_exps = Vector{Vector{Int}}(undef, length(basis))
    basis_dict = Dict{Vector{Int}, Int}()
    sizehint!(basis_dict, length(basis))

    for (i, b) in enumerate(basis)
        e = MP.exponents(b) # so that we allocate exponents only once
        basis_exps[i] = e
        basis_dict[e] = i
    end

    return OnPoints(basis_exps, basis_dict)
end

and this is how it works in practice:

C3 = PermGroup(perm"(1,2,3)")
chars = SymbolicWedderburn.characters_dixon(C3)
mvec = reverse(collect(monomials(x, 0:2)), 1)

let G = C3, chars = chars, basis = mvec
    induced_action = OnPoints(basis)

    S = gens(G)
    generators_large = Vector{Perm{Int}}(undef, length(S))
    for (i, s) in enumerate(S)
        generators_large[i] = induced_action(s)
    end

    G_large = PermGroup(generators_large)

    ccG_large = conjugacy_classes(G_large)

    # TODO: permute ccG_large if necessary to match order of
    # conjugacy_classes(charse[i])

    for i in 1:length(chars)
        conjugacy_classes(chars[i]) .= ccG_large
    end
end

---

Then after defining

c1, U1 = SymbolicWedderburn.central_projection(chars[1])
c2, U2 = SymbolicWedderburn.central_projection(chars[2])
c3, U3 = SymbolicWedderburn.central_projection(chars[3])

we get:

julia> R1, ids1 = let U = U1
           SymbolicWedderburn.row_echelon_form(float.(U))
       end; R1
4×4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> R2, ids2 = let U = U2 + conj(U2)
           SymbolicWedderburn.row_echelon_form(float.(U))
       end; R2
4×4 Array{Float64,2}:
 0.0  1.0  0.0  -1.0
 0.0  0.0  1.0  -1.0
 0.0  0.0  0.0   0.0
 0.0  0.0  0.0   0.0

julia> R3, ids3 = let U =-E(4)*(U2 - conj(U2)) # -E(4) = -im = (im)^-1
           SymbolicWedderburn.row_echelon_form(float.(U))
       end; R3
4×4 Array{Float64,2}:
 0.0  1.0  0.0  -1.0
 0.0  0.0  1.0  -1.0
 0.0  0.0  0.0   0.0
 0.0  0.0  0.0   0.0

U2 + U3 corresponds to the class function char[2] + char[3] which in turn correspond to the direct sum of representations defined by char[2] and char[3]; Actually I should allow such arithmetic on the level of characters so that we could write c2, U2 = SymbolicWedderburn.central_projection(chars[2] + chars[3]) ;)
Since char[2] and char[3] are conjugate we have U3 = conj(U2), and U2+conj(U2) is guaranteed to be real.
By a similar trick to (x,y) → (x+y, x-y) we can guarantee that (chars[2] - conj(chars[2])) is purely imaginary, so U3 = -E(4)*(chars[2] - conj(chars[2])) is real.

It seems that we'd need to look for conjugate pairs of representations when it comes to defining real projections.

Anyway, since U2 + conj(U2) is actually 2-dimensional representation the third one should be either 0-dimensional, or the same (and this is what we see in R3: it's projection onto the same subspace).

I don't know if this makes much sense to You, but it's late here and I'm going to bed ;)

btw. I didn't want to step on your toes and push to your branch, so I have my own with central_projection: https://github.com/kalmarek/SymbolicWedderburn.jl/tree/enh/central_projections (feel free to merge it/to it whatever you wish).

And I get a strange error:

julia> msym = SOSModel(SCS.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: SCS

julia> @variable   msym t
t

julia> @objective  msym Max t
t

julia> @variable   msym sos1 SOSPoly(FixedPolynomialBasis(R1[1:length(ids1),:]*mvec))
GramMatrix{VariableRef,FixedPolynomialBasis{Polynomial{true,Float64},Array{Polynomial{true,Float64},1}},GenericAffExpr{Float64,VariableRef}}(VariableRef[noname noname; noname noname], FixedPolynomialBasis{Polynomial{true,Float64},Array{Polynomial{true,Float64},1}}(Polynomial{true,Float64}[1.0, x₁ + x₂ + x₃]))

julia> @variable   msym sos2 SOSPoly(FixedPolynomialBasis(R2[1:length(ids2),:]*mvec))
GramMatrix{VariableRef,FixedPolynomialBasis{Polynomial{true,Float64},Array{Polynomial{true,Float64},1}},GenericAffExpr{Float64,VariableRef}}(VariableRef[noname noname; noname noname], FixedPolynomialBasis{Polynomial{true,Float64},Array{Polynomial{true,Float64},1}}(Polynomial{true,Float64}[x₁ - x₃, x₂ - x₃]))

julia> @constraint msym sum(x.^2) - t == sos1 + 0.5*sos2
ERROR: ArgumentError: reducing over an empty collection is not allowed
Stacktrace:
 [1] _empty_reduce_error() at ./reduce.jl:299
 [2] mapreduce_empty(::Function, ::Function, ::Type{T} where T) at ./reduce.jl:342
 [3] reduce_empty(::Base.MappingRF{MultivariateBases.var"#9#11"{SymMatrix{VariableRef},FixedPolynomialBasis{Polynomial{true,Float64},Array{Polynomial{true,Float64},1}},Int64},typeof(MutableArithmetics.add!)}, ::Type{Int64}) at ./reduce.jl:329
 [4] reduce_empty_iter at ./reduce.jl:355 [inlined]
 [5] mapreduce_empty_iter(::Function, ::Function, ::UnitRange{Int64}, ::Base.HasEltype) at ./reduce.jl:351
 [6] _mapreduce(::MultivariateBases.var"#9#11"{SymMatrix{VariableRef},FixedPolynomialBasis{Polynomial{true,Float64},Array{Polynomial{true,Float64},1}},Int64}, ::typeof(MutableArithmetics.add!), ::IndexLinear, ::UnitRange{Int64}) at ./reduce.jl:400
 [7] _mapreduce_dim(::Function, ::Function, ::NamedTuple{(),Tuple{}}, ::UnitRange{Int64}, ::Colon) at ./reducedim.jl:318
 [8] mapreduce(::Function, ::Function, ::UnitRange{Int64}; dims::Function, kw::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at ./reducedim.jl:310
 [9] mapreduce at ./reducedim.jl:310 [inlined]
 [10] polynomial(::SymMatrix{VariableRef}, ::FixedPolynomialBasis{Polynomial{true,Float64},Array{Polynomial{true,Float64},1}}, ::Type{T} where T) at /home/kalmar/.julia/packages/MultivariateBases/KKez4/src/fixed.jl:26
 [11] polynomial at /home/kalmar/.julia/packages/SumOfSquares/61JEU/src/gram_matrix.jl:77 [inlined]
 [12] polynomial at /home/kalmar/.julia/packages/SumOfSquares/61JEU/src/gram_matrix.jl:74 [inlined]
 [13] _multconstant(::Float64, ::Function, ::GramMatrix{VariableRef,FixedPolynomialBasis{Polynomial{true,Float64},Array{Polynomial{true,Float64},1}},GenericAffExpr{Float64,VariableRef}}) at /home/kalmar/.julia/packages/MultivariatePolynomials/CZWzD/src/operators.jl:248
 [14] multconstant at /home/kalmar/.julia/packages/MultivariatePolynomials/CZWzD/src/operators.jl:250 [inlined]
 [15] *(::Float64, ::GramMatrix{VariableRef,FixedPolynomialBasis{Polynomial{true,Float64},Array{Polynomial{true,Float64},1}},GenericAffExpr{Float64,VariableRef}}) at /home/kalmar/.julia/packages/MultivariatePolynomials/CZWzD/src/operators.jl:46
 [16] promote_operation(::typeof(*), ::Type{T} where T, ::Type{T} where T) at /home/kalmar/.julia/packages/MutableArithmetics/NuiNA/src/interface.jl:21
 [17] promote_operation(::typeof(MutableArithmetics.sub_mul), ::Type{T} where T, ::Type{T} where T, ::Type{T} where T) at /home/kalmar/.julia/packages/MutableArithmetics/NuiNA/src/shortcuts.jl:53
 [18] mutability(::Type{T} where T, ::Function, ::Type{T} where T, ::Type{T} where T, ::Type{T} where T) at /home/kalmar/.julia/packages/MutableArithmetics/NuiNA/src/interface.jl:132
 [19] mutability at /home/kalmar/.julia/packages/MutableArithmetics/NuiNA/src/interface.jl:138 [inlined]
 [20] operate!(::typeof(MutableArithmetics.sub_mul), ::Polynomial{true,GenericAffExpr{Float64,VariableRef}}, ::Float64, ::GramMatrix{VariableRef,FixedPolynomialBasis{Polynomial{true,Float64},Array{Polynomial{true,Float64},1}},GenericAffExpr{Float64,VariableRef}}) at /home/kalmar/.julia/packages/MutableArithmetics/NuiNA/src/rewrite.jl:70
 [21] top-level scope at /home/kalmar/.julia/packages/MutableArithmetics/NuiNA/src/rewrite.jl:227
 [22] top-level scope at /home/kalmar/.julia/packages/JuMP/YXK4e/src/macros.jl:416

??? (0.5 is the culprit)

[kalmarek]

I'm not really sure that what I wrote is correct with respect to U2 + conj(U2), etc, since it doesn't work for e.g. C5. (btw. I don't even know what to expect for this group??) but anyway ;)

---

I must have accidentally hit the right design with OnPoints, since it can act as a homomorphism, so that we can:

abstract type AbstractAction{T} end

struct OnPoints{T} <: AbstractAction{T}
    features::Vector{T}
    reversef::Dict{T, Int}
end

Base.getindex(act::OnPoints, i::Integer) = act.features[i]
Base.getindex(act::OnPoints{T}, f::T) where T = act.reversef[f]
(act::OnPoints)(p::Perm{I}) where I = Perm(I[act[f^p] for f in act.features])

function OnPoints(basis::Union{MonomialVector, AbstractVector{<:Monomial}})
    basis_exps = Vector{Vector{Int}}(undef, length(basis))
    basis_dict = Dict{Vector{Int}, Int}()
    sizehint!(basis_dict, length(basis))

    for (i, b) in enumerate(basis)
        e = MP.exponents(b) # so that we allocate exponents only once
        basis_exps[i] = e
        basis_dict[e] = i
    end

    return OnPoints(basis_exps, basis_dict)
end

function (act::OnPoints)(orb::O) where {T, O<:AbstractOrbit{T, Nothing}}
    elts = act.(orb)
    dict = Dict(e=>nothing for e in elts)
    return O(elts, dict)
end

---

And now a very natural way of using OnPoints as a homomorphism:

include("test/smallgroups.jl")
G = SmallPermGroups[6][1]
@assert length(conjugacy_classes(G)) == 3
chars = SymbolicWedderburn.characters_dixon(G)

@polyvar x[1:6]
mvec = reverse(monomials(x, 0:2))

chars_large = let chars = chars, basis = mvec

    @assert all(χ.inv_of == first(chars).inv_of for χ in chars)

    induced_action = OnPoints(basis)
    ccG_large = incuded_action.(conjugacy_classes(first(chars)))

    # just to double check:
    let ccls = ccG_large, large_gens = induced_action.(gens(G))
        G_large = PermGroup(large_gens)
        ccG_large = conjugacy_classes(G_large)
        @assert all(Set.(collect.(ccG_large)) .== Set.(collect.(ccls)))
    end

    [SymbolicWedderburn.Character(χ.vals, χ.inv_of, ccG_large) for χ in chars]
end

U = [SymbolicWedderburn.central_projection(χ) for χ in chars_large]

R = map(U) do c_u
    u = last(c_u)
    if all(isreal, u)
        SymbolicWedderburn.row_echelon_form(float.(u))
    else
        throw("Not Implemented")
    end
end

which results in

julia> pr, ids = R[1]; pr[1:length(ids), :]*mvec
18-element Array{Polynomial{true,Float64},1}:
 -x₂ + x₆
 -x₁ + x₅
 -x₂ + x₄
 -x₁ + x₃
 -x₂² + x₆²
 -x₁x₂ + x₅x₆
 -x₁² + x₅²
 -x₂x₄ + x₄x₆
 -x₁x₆ + x₄x₅
 -x₂² + x₄²
 -x₁x₄ + x₃x₆
 -x₁x₃ + x₃x₅
 -x₁x₂ + x₃x₄
 -x₁² + x₃²
 -x₂x₄ + x₂x₆
 -x₁x₄ + x₂x₅
 -x₁x₆ + x₂x₃
 -x₁x₃ + x₁x₅

julia> pr, ids = R[2]; pr[1:length(ids), :]*mvec # truly invariant polynomials
7-element Array{Polynomial{true,Float64},1}:
 1.0
 x₁ + x₂ + x₃ + x₄ + x₅ + x₆
 x₁² + x₂² + x₃² + x₄² + x₅² + x₆²
 x₁x₂ + x₃x₄ + x₅x₆
 x₁x₃ + x₁x₅ + x₂x₄ + x₂x₆ + x₃x₅ + x₄x₆
 x₁x₆ + x₂x₃ + x₄x₅
 x₁x₄ + x₂x₅ + x₃x₆

julia> pr, ids = R[3]; pr[1:length(ids), :]*mvec # alternating polynomials
3-element Array{Polynomial{true,Float64},1}:
 -x₁ + x₂ - x₃ + x₄ - x₅ + x₆
 -x₁² + x₂² - x₃² + x₄² - x₅² + x₆²
 -x₁x₃ - x₁x₅ + x₂x₄ + x₂x₆ - x₃x₅ + x₄x₆

[tweisser]

I'm working on it