Sage for Specht Modules

Appropos of nothing, these were some cute monkeys playing on the water tank this morning. I actually got a very close look at them, but by the time the camera was out, they were scurrying onto the roof.

As the strike drags on, I’ve had some time to actually do some maths.  In particular, I’m preparing to run a weekly seminar on using representation theory for certain statistical problems.  The plan is to work from some old lecture notes by Persi Diaconis entitled ‘Group Representations in Probability and Statistics.’ These deal with, for example, how long one should apply a random shuffling process before the thing which is being shuffled is well-mixed.  Particular examples include the question of how many times one should shuffle a deck of cards, and how long one should let a random walk on Z_n run before we can be reasonably sure that every point has been reached.  There are numerous real-world applications of the results, and it uses a lot of first-rate representation theory along the way!

As I’ve been reading, I thought it would be really nice to have some sage code cooked up to demonstrate the ideas in the text.  I then realized that there doesn’t seem to be a nice implementation of the irreducible representations of the symmetric group (aka, Specht modules) in Sage, though many of the pieces are there.  You can get the character table for the symmetric group, but this is actually rather less than a full combinatorial implementation of the Specht modules!  That existing character table code was nice to have, since it makes it easy to check that my algorithms are actually producing sensible results.

So I spent the last day or so writing some simple code to do the trick.  It’s basically an implementation of chapter two of Bruce Sagan’s (excellent) book on the Symmetric Group, using a bunch of the pre-existing machinery in Sage to handle things like standard tableaux.

One nice thing that I implemented which could be worth putting into main sage was the dominance order on tableaux, which is useful for iterating over the polytabloids that make up a basis of the Specht modules.  This order has some nice triangularity properties which make it easy to get the coefficients that appear in a matrix for a permutation acting on the Specht module.

The reason for implementing all this was to get those matrices; in the notes, there’s a notion of a Fourier transform of a probability function with respect to a given representation.  The formula includes a sum over all of the matrices for the permutations in the given representation.  There maybe a smarter way to do it (well, almost certainly given how slow my code is just for S_5), and I’ll update as it becomes clear…

Ok, that’s all for now!  Here’s what I’ve got so far.

def StandardTableauxDominancePoset(lam):
    #This poset is used for triangularity in the Specht Module.
    def composition_tuple(s):
        #s a standard tableau of shape lam.
        #last composition is always lam, so no need to include it here...
        Mu=[]
        for i in range(lam.size()):
            Mui = [ len([j for j in r if j<=i+1]) for r in list(s)]
            Mu.append( Composition( Mui ) )
        return tuple(Mu)

    def dominates(s,t):
        #Checks if s dominates t.
        cs=composition_tuple(s)
        ct=composition_tuple(t)
        for i in range(len(cs)):
            cs_partial_sums=cs[i].partial_sums()
            ct_partial_sums=ct[i].partial_sums()
            #print i, 't', cs_partial_sums, 't', ct_partial_sums
            for j in range(len(cs_partial_sums)):
                if not cs_partial_sums[j]>=ct_partial_sums[j]: return False
        return True

    st=StandardTableaux(lam)
    G=DiGraph()
    if st.cardinality()==1:
        G.add_vertex(st[0])
        return Poset(G)
    for s in st:
        for t in st:
            if dominates(s,t): 
                G.add_edge((s,t))
                #print s,t
    return Poset(G)

class SpechtModule(CombinatorialFreeModule):
    def __init__(self, shape, R=QQ):
        self.shape=Partition(shape)
        self.size=shape.size()
        self.S=SymmetricGroup(self.size)
        self.A=self.S.algebra()
        self.st=StandardTableaux(lam)
        L=StandardTableauxDominancePoset(lam).linear_extension()
        #L.reverse()
        self.ordered_st=[t.element for t in L]

        CombinatorialFreeModule.__init__(self, R, Tableaux(shape=shape))

    def __repr__(self):
        return "Tableau module of shape " + str(self.shape)

    def _permutation_action_on_tableau(self, p, t):
        #t a tableau
        #p a permutation
        s=t.to_list()
        n=t.size()
        for i in [1..n]:
            for c in t.cells_containing(i):
                s[c[0]][c[1]]=p(i)
        return self(Tableau(s))

    def _algebra_action_on_tableau(self, x, t):
        #x an element of the group algebra
        #t a tableau
        action=self._permutation_action_on_tableau
        return sum([x.coefficient(p)*action(p,t) for p in x.support()])

    def action(self, x, t):
        #x an element of the Sn group algebra
        #t an element of self.
        x=x.parent().algebra()(x)
        algAct=self._algebra_action_on_tableau
        return sum([t.coefficient(s)*algAct(x, s) for s in t.support()])

    @cached_method
    def polytabloid(self, t):
        cs=t.column_stabilizer()
        ct=sum([a.sign()*self.A(a) for a in cs])
        rs=t.row_stabilizer()
        rt=sum([self.A(a) for a in rs])
        b=rt*ct
        return self._algebra_action_on_tableau(b, t)

    @cached_method
    def spechtBasis(self):
        return [self.polytabloid(Tableau(t)) for t in self.ordered_st]

    def gens(self):
        return self.spechtBasis()

    def matrix_of_permutation_on_basis(self,p):
        m=[]
        for b in self.spechtBasis():
            a=self.action(p,b)
            coeffs=[]
            for t in self.ordered_st:
                #print t, a.coefficient(t)
                c=a.coefficient(t)
                coeffs.append(c)
                if c!=0: a=a-c*self.polytabloid(t)
            assert a==0, "Something's gone terribly wrong."
            m.append(coeffs[:])
        #m.reverse()
        return matrix(m)

    def character(self):
        ch=[]
        for pi in self.S.conjugacy_classes_representatives():
            ch.append( self.matrix_of_permutation_on_basis(pi).trace() )
        return ch

    def fourierTransform(self, f):
        #f a function from the permutation group to the base field.
        mat=self.matrix_of_permutation_on_basis
        return sum([f(p)*mat(p) for p in self.S])

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

n=lam.size()
S=SymmetricGroup(n)
x=var('x')

def P_reflect(s):
    assert s in S, "Need permutations. "+str(s) + str(S)
    if s==S.one(): return x
    if s^2==S.one(): return (1-x)/binomial(n,2)
    return 0

convolve = lambda P,Q: lambda s: sum([P(s.inverse()*t)*Q(t) for t in S])
Advertisements