# Lie-algebra representation powers - plethysms

#### lpetrich

I've done some more work on my http://homepage.mac.com/lpetrich/Science/SemisimpleLieAlgebras.zip [Broken] package, adding decompositions of representation powers. It was only recently that I discovered that they are called "plethysms".

Rep powers can be decomposed by symmetry types, a feature that can be important in their applications. For power p, the types are symmetric (p >=1), antisymmetric (p >= 2), and mixed (p >= 3), with more than one mixed type for p >= 4.

These are related to the irreps in (SU(n) fundamental rep)p, and they share multiplicities. The symmetry types and irreps can be denoted with Young diagrams having p boxes.

With a lot of experimentation and trial and error, I have worked out how to do rep-power decompositions, but the algorithm is hideously complicated, Though I've successfully tested it for low powers, I don't have a rigorous proof for it, and I've been unsuccessful in finding an algorithm in the professional literature. I've also had to create versions for reducible reps and algebra-product reps, because rep-power decompositions do not have any simple decompositions of those in terms of the rpd's of the component irreps of the component algebras.

-

As an example, consider flavor and spin symmetries of the three light quarks: SU(3) and SU(2) respectively. They can be combined into SU(6), and one can undo this combination into SU(3)*SU(2) again with

Mathematica: SubalgMultSU[ldqksp,{3,2}]; res = DoBranching[ldqksp,{1,0,0,0,0}]
Python: ldqksp = SubalgMultSU((3,2)); res = ldqksp.DoBranching((1,0,0,0,0))

Breakdown of cube of fundamental rep of SU(6):

Mathematica: MakeLieAlgebra[laqksp,{1,5}]; res = DecomposeRepPower[laqksp,{1,0,0,0,0},3]
Python: res = DecomposeRepPower((1,5), (1,0,0,0,0), 3)

Count, Irrep:
S: (1, (3, 0, 0, 0, 0))
M: (1, (1, 1, 0, 0, 0))
A: (1, (0, 0, 1, 0, 0))

Breakdown of cube of corresponding rep of SU(3)*SU(2):

Mathematica: MakeLieAlgebra[laqk, {1, 2}]; MakeLieAlgebra[lasp, {1, 1}]; res = DecomposeAlgProdRepPower[{laqk, lasp}, {{1, 0}, {1}}, 3]
Python: res = DecomposeAlgProdRepPower(((1,2),(1,1)), ((1,0),(1,)), 3)

Count, SU(3) Irrep, SU(2) Irrep
S: (1, ((3, 0), (3,))) . (1, ((1, 1), (1,)))
M: (1, ((1, 1), (3,))) . (1, ((3, 0), (1,))) . (1, ((1, 1), (1,))) . (1, ((0, 0), (1,)))
A: (1, ((1, 1), (1,))) . (1, ((0, 0), (3,)))

-

I'll now translate these results into typical physicist language, using rep multiplicities and spins:

Each quark: (3,1/2)
S: (10,3/2) + (8,1/2)
M: (8,3/2) + (10,1/2) + (8,1/2) + (1,1/2)
A: (8,1/2) + (1,3/2)

The symmetric one is the one that matches the light-baryon spectrum, and not the mixed or the antisymmetric ones. This is contrary to Fermi-Dirac statistics, and it took QCD to resolve that discrepancy. The quarks are antisymmetric in color, which gives them overall antisymmetry, thus preserving FD statistics.

Last edited by a moderator:
Related Linear and Abstract Algebra News on Phys.org

#### Simon_Tyler

Hi Ipetrich,

Glad to see you're still working on this - it's looking good (although it's getting beyond the stuff that I've played with). I like how you're maintaining the Mathematica and python code and always giving examples from QFT.

Branching rules from plethysms are implemented in sage.
It's described here (with references):
http://www.sagemath.org/doc/thematic_tutorials/lie/branching_rules.html#plethysms
and the code can be found here:

This should provide some good comparisons/test cases for your code.

Last edited by a moderator:

#### lpetrich

Thanx. I checked on that, and it does not seem to list plethysm decompositions by symmetry type.

That was part of the algorithms for branching rules, and I'd implemented all of them:
• Levi -- root demotion. Sage's code does not seem to calculate the resulting U(1) factors, however.
• Dynkin-diagram extension splitting.
• Division of orthogonal and symplectic groups: SO(m+n)->SO(m)*SO(n), Sp(2(m+n)) -> Sp(2m)*Sp(2n). All but SO(even) -> SO(odd)*SO(odd) can be derived from root demotion and extension splitting, however, using SO(2) ~ U(1).
• Symmetric subgroups: SU(n) -> SO(n), SU(2n) -> Sp(2n), SO(2n) -> SO(2n-1), E(6) -> F(4)
• Tensor products: SU(m*n) -> SU(m)*SU(n), SO/Sp(m*n) -> SO/Sp(m)*SO/Sp(n) (SO -> SO*SO or Sp*Sp, Sp -> SO*Sp or Sp*SO)
• Symmetric powers: I implement what seems like a superset of their special case: anything -> SU(2). Root height -> SU(2) root value.
The plethysm rule seemed like a subcase of one of Dynkin's theorems. I dug up a translation of a collection of his papers, and I found various useful theorems.

I've also implemented all the maximal subalgebras for the exceptional ones; Slansky lists all the maximal subalgebras for algebras with rank <= 8.

These special cases do not encompass all the subalgebras of classical algebras, but Dynkin has some theorems that supply constraints on the remaining ones:

* At least one "simple" irrep must map onto a single term in the subalgebra, but that term need not have a "simple" irrep. Here are the "simple" ones:

A(n) / SU(n-1): vector and conjugate vector (dim 2n+1)
B(n) / SO(2n+1): vector (dim 2n+1) and spinor (dim 2^n)
C(n) / Sp(2n): vector (dim 2n)
D(n) / SO(2n): vector (dim 2n) and the two spinors (dim 2^(n-1))

* The algebra's adjoint rep must map onto subalgebra's adjoint rep, though usually with extra irreps. The adjoint can be constructed with the help of powers and symmetry types (plethysms) of the defining irreps for the original and the subalgebra.

* The original algebra's irrep and the subalgebra's one must have the same "reality":
• Real -- self-conjugate with even height -- adjoint in antisymmetric part of square
• Pseudoreal or quaternionic -- self-conjugate with odd height -- adjoint in symmetric part of square
• Complex -- not self-conjugate -- adjoint in higher power
In general, (irrep)*(its conjugate) always contains the adjoint rep.

I could not find any general algorithm for calculating projection matrices in this case, or any subcases, so I left unimplemented most of these ones that Slansky had listed. Some examples of this case:

SU(6) -> SU(3): (1,0,0,0,0) 6 -> (2,0) 6
SO(7) -> G(2): (1,0,0) 7 -> (0,1) 7

#### lpetrich

I'd found branching-rule projection matrices in ad hoc fashion; part of the problem is finding which subalgebra rep root an original-algebra rep root maps onto.

Turning to Sage's D7 -> G2 plethysm example, I note that D7's vector rep (1,0,0,0,0,0,0) maps onto G2's adjoint rep (1,0), both with dimension 14. For G2, Sage uses reverse order from my code.

In the example's final output, the first five entries were all antisymmetric products, and they agree with what my code produces. I tried my code on D7 also, and the antisymmetric products were all single irreps with one weight value 1 and the others zero. Going to the 6th power yielded (0,0,0,0,0,1,1) -- the product of the two spinors. So mapping an individual spinor onto G2 will be tricky.

In my code, I identify the symmetry types with Young diagrams, and these map directly onto the irreps in powers of the fundamental rep of SU(n) / A(n+1).
Power p symmetric: (p)
Power p antisymmetric: (p 1's)
Power p mixed: (other values that add up to p)

Symmetric: all index permutations equally weighted
Antisymmetric: like symmetric, but all odd permutations have reverse sign
Mixed: hard to find any simple formulas

Getting back to Young diagrams, one can easily find SU(n) rep dimensions with them, and there's a moderately complicated algorithm called the Littlewood-Richardson rule for rep products. Are there similar techniques for SO(n) and Sp(2n)? For SO(n), product irreps of the vector rep have to be traceless, and one has to take that into account. Like:

SU(n) -- vector * vector = (symmetric tensor) + (antisymmetric tensor)
SO(n) -- the symmetric tensor becomes (symmetric traceless tensor) + (scalar)
Sp(2n) -- the antisymmetric tensor becomes (antisymmetric sort-of-traceless tensor) + (scalar)

I say sort-of-traceless, because it's an operation that involves subtracting out the symplectic form for the algebra:
J = {{0, I}, {-I, 0}}
in the fashion of tracelessness, subtracting out the identity matrix.

#### lpetrich

I have updated my Lie-algebra code yet again, this time to support U(1) factors in the algebra-list versions of my rep powers. U(1) factors can be any integers or rational numbers; they follow the max-weight vectors. It's possible to have only U(1) factors if one so desires. Here, I'll give some results in Python form - it should be easy to translate ino Mma form.

DecomposeAlgProdRepPower( ((1,1),), ((1,),42), 2)
yields
((1, ((2,), 84)),)
((1, ((0,), 84)),)
(adding two spin-1/2's to make a spin 1 (symmetric) and a spin 0 (antisymmetric))

Now to raise vectors in various dimensions in various powers. The outputs are the symmetry-type decomposition lists, whose members have format (count, rep).

2D vector: SO(2) = U(1) -- helicity. Demo of using no semisimple algebras, just U(1)'s
DecomposeAlgProdRepListPower,( (), ((1,),(-1,)), p)
p = 1
((1, (1,)), (1, (-1,)))
p = 2
((1, (2,)), (1, (0,)), (1, (-2,)))
((1, (0,)),)
p = 3
((1, (3,)), (1, (1,)), (1, (-1,)), (1, (-3,)))
((1, (1,)), (1, (-1,)))
()
p = 4
((1, (4,)), (1, (2,)), (1, (0,)), (1, (-2,)), (1, (-4,)))
((1, (2,)), (1, (0,)), (1, (-2,)))
((1, (0,)),)
()
()

3D vector: SO(3) = B(1)
DecomposeRepPower( (2,1), (2,), p)
p = 1
((1, (2,)),)
p = 2
((1, (4,)), (1, (0,)))
((1, (2,)),)
p = 3
((1, (6,)), (1, (2,)))
((1, (4,)), (1, (2,)))
((1, (0,)),)
p = 4
((1, (8,)), (1, (4,)), (1, (0,)))
((1, (6,)), (1, (4,)), (1, (2,)))
((1, (4,)), (1, (0,)))
((1, (2,)),)
()

4D vector: SO(4) = D(2)
DecomposeRepPower( (4,2), (1,1), p)
p = 1
((1, (1, 1)),)
p = 2
((1, (2, 2)), (1, (0, 0)))
((1, (0, 2)), (1, (2, 0)))
p = 3
((1, (3, 3)), (1, (1, 1)))
((1, (1, 3)), (1, (3, 1)), (1, (1, 1)))
((1, (1, 1)),)
p = 4
((1, (4, 4)), (1, (2, 2)), (1, (0, 0)))
((1, (2, 4)), (1, (4, 2)), (1, (2, 2)), (1, (2, 0)), (1, (0, 2)))
((1, (2, 2)), (1, (4, 0)), (1, (0, 4)), (1, (0, 0)))
((1, (2, 2)), (1, (2, 0)), (1, (0, 2)))
((1, (0, 0)),)

The second one for p = 2 is the antisymmetric 2-tensor.
For 2D; scalar: 0
For 3D: vector: (2,)
For 4D: self-dual and anti-self-dual parts: (0, 2) and (2, 0)
Duality: (1/2)*e.A = +/- A

For the Faraday tensor, the self-dual and anti-self-dual parts are, in 3+1 notation, E+i*B and E-i*B

The third one for p = 4 has the symmetries of the Riemann tensor. Its Young diagram is (2,2)
For 2D, it's given by the Ricci scalar: 0
For 3D, it's given by the Ricci tensor (symmetric 2-tensor): (4,) and (0,)
For 4D it's given by the Ricci tensor + the Weyl tensor:
(2, 2) and (0, 0)
(4, 0) and (0, 4)
The latter has self-dual and anti-self-dual parts.

#### lpetrich

I've fixed a bug in the Python version, SemisimpleLieAlgebras.py -- DecomposeRepProduct had been returning empty products, but it's now returning correct ones. I've also added counted-list inputs to rep powers, those being lists of (count,maxwts) instead of just maxwts.

I also have some more conjectures in addition to my conjectured rep-power-decomposition (plethysm) formula. I've looked in the literature, but I can't find whether they've been anticipated.

A rep-product one:

Take representation R and some power, and take the part with the symmetry type denoted by Young diagram Y. I will call this part RY (R to power Y).

The rep product RY1 * RY2 = sumY cY1,Y2,Y RY

where cY1,Y2,Y are the Littlewood-Richardson coefficients for how many copies of Y are in the product of Y1 and Y2.

This is correct if R is the fundamental representation of SU(n), and I've verified it in some other cases. I've looked in the literature for a general proof, without success.

A rep-sum one:

For representations R1 and R2 of a Lie algebra, the sum rep R1 + R2 satisfies
(R1 + R2)Y = sumY1,Y2 cY1,Y2,Y R1Y1 * R2Y2

I've verified that one also in some cases. I'm in the design phase of implementing it more generally, and I should be able to implement it before long.

Testing cannot be complete, of course, but it should be feasible for the lower-rank algebras and lower-multiplicity reps.

Others?

I've found counterexamples to (RY1)Y2 = (RY2)Y1, and I've been unable to guess what cx is in

(R1,R2)Y = sumY1,Y2 cxY1,Y2,Y (R1Y1,R2Y2)

other than Y1 and Y2 having the same order as Y for cx to be nonzero.

#### lpetrich

I've updated http://homepage.mac.com/lpetrich/Science/SemisimpleLieAlgebras.zip [Broken] to include code for testing my previous post's conjectures in the Mathematica version. I've run it on several examples, and it has been successful so far. I still haven't found proofs, however.

But both Mathematica and Python versions now have code for doing purely-symmetric or purely-antisymmetric representation-power decompositions (plethysms).

Example in Python:
Code:
# Lie algebras for color, 3-quark flavor, and spin
# Outer product because a quark has all 3 degrees of freedom
lax = ((1,2), (1,2), (1,1))

# Quark rep: color 3 * flavor 3 * spin mult. 2 (spin 1/2)
qkwt = ((1,0), (1,0), (1,))

# The 3 is for there being 3 quarks in a baryon
# The -1 is for being completely antisymmetric (Fermi-Dirac statistics)
baryonall = DecomposeAlgProdRepPwrSym(lax,qkwt,3,-1)

# Extract the color singlets
baryon = [b for b in baryonall if b == (0,0)]

# Display the baryon reps
for b in baryon: print b

# Multiplicity, (color rep, flavor rep, spin rep)
(1, ((0, 0), (3, 0), (3,)))
(1, ((0, 0), (1, 1), (1,)))

# Equivalent to flavor 10 (triangle), spin 3/2 and flavor 8 (hexagon), spin 1/2
The Sage results that Simon_Tyler had linked to I could duplicate with (Mma)

MakeLieAlgebra[lag2, {7,2}]
DecomposeRepPwrSym[lag2, {1,0}, pwr, -1]

with pwr = 1,2,3,4,5
Sage's code uses the reverse order of roots/weights from mine, which is Robert Cahn's order.

Last edited by a moderator:

#### lpetrich

I've added some stuff on Schur functions and related ones to my Young-diagram Mma notebook in that archive. I've done so because Schur functions are related to the irreps of SU(n) / A(n-1) / SL(n,R). They are even specified in the same way, with partitions / Young diagrams. In general:

Schur function: S(part / YD, vector)

The smallest ones:
S({},{x1,x2,x3}) = 1
S({1},{x1,x2,x3}) = x1+x2+x3
S({2},{x1,x2,x3}) = x1^2+x1*x2+x2^2+x1*x3+x2*x3+x3^2
S({1,1},{x1,x2,x3}) = x1*x2+x1*x3+x2*x3

They form a basis for the symmetric polynomials in their vector arg.

Products of them give linear combinations of them specified with the Littlewood-Richardson rule, familiar from products of SU(n) irreps.

Plethysms on Schur functions are defined with this algorithm:

To calculate S(a)[S(b)] the plethysm of S(a) over S(b), do
For a vector x, decompose S(b,x) as a sum of monic monomials in products and powers of the x components: x(1)^2*x(2), etc.
Make them into a vector, and evaluate S(a,(this monic-monomial vector))

You can now decompose this value into a sum of S(c,x) values. In short,
S(a)[S(b)] = sum over c of pl(a,b,c) S(c)

One can replace S(a) with some other sort of symmetric function on a vector, and still get (Schur function) -> (combination of Schur functions).

However, it's been difficult for me to find an algorithm for calculating pl(a,b,c) comparable to the Littlewood-Richardson rule. One can do it with the algorithm here, or with A(n) Lie-algebra irreps, but I wish to avoid these algorithms, because they can require a very large amount of calculation.

#### lpetrich

I've succeeded in finding an algorithm for pl(a,b,c) when a and b are partitions / Young diagrams for symmetric or antisymmetric SU(n) reps. That is, {n} or {n 1's} / {1^n}

This algorithm sidesteps calculating reps or working with Schur-function vector args, thus making it much faster. It's essentially the algorithm in EUROCAL '85. European Conference on ... - Bob F. Caviness - Google Books

Some special cases of Schur functions:
S({n},x) = homogeneous symmetric polynomial in x of order n
Terms have the form (x1^p1)*(x2^p2)*... powers summing to n
S({1^n},x) = elementary symmetric polynomial in x of order n
Terms have the form x1*x2*... powers no greater than 1 and summing to n
All coefficients are 1

Calculate S({n},x^m) and S({1^n},x^m} from power sums of x using (power-sum to homogeneous) and (power-sum to elementary) formulas.

Calculate the plethysm by treating the above S's as power sums in the (power-sum to homogeneous) and (power-sum to elementary) formulas.

Along the way, I found a formula for power sums in terms of Schur functions:
sum(x^n) = sum over k: 0 to n-1 of (-1)^k * S({n-k, k 1's},x)

I was unable to prove it, though I was able to verify it for small n's.

#### lpetrich

Various identities involving plethysms: rep-power decompositions.

http://www.sbfisica.org.br/bjp/files/v32_641.pdf
The Plethysm Technique Applied to the Classification of Nuclear States, by Alcaras, Tambergs, Krasta, Ruza, Katkevicius

plt(a,x) is the plethysm of a on x, with a the plethysmer and x the plethysmed one, an irrep. a is a partition / Young diagram. A power p decomposes into diagrams a with length p.

Inner product: a.b = sum of ipc(a,b,c) * c

ipc(a,b,c) = 1/N sum over diagrams r of N(r)*X(a,r)*X(b,r)*X(c,r)
X(a,r) = characters of symmetric group for irrep a and class b (always integer), N(r) is order of class with part/diag r, N is total order.

Outer product: a*b = sum of opc(a,b,c) * c

opc(a,b,c) one calculates uses the Littlewood-Richardson rule.

Reps x,y have sum (sum of coefficients of irreps in it), difference (likewise with differences), product (rep product)

The identities:

plt(a +- b, x) = plt(a,x) +- plt(b,x)
plt(a, x + y) = sum over b,c of opc(b,c,a) * plt(b,x)*plt(c,y)
plt(a, x - y) = sum over b,c of opc(b,c,a) * (-1)^length(c) * plt(b,x)*plt(c,y)
plt(a*b,x) = plt(a,x)*plt(b,x)
plt(a,x*y) = sum over b,c of ipc(b,c,a) * plt(b,x)*plt(c,y)
conjg(plt(a,x)) = plt(a,conjg(x))

I earlier went through a lot of trouble to find formulas for plt(a,x+y) and plt(a,x*y), and while I successfully guessed the former one by extrapolating from small examples, I had much more difficulty with the latter one.

Last edited:

#### lpetrich

I now have plethysms / rep-power decompositions working on Young diagrams, without calculating Lie-algebra reps or Schur functions of vectors.

SelYDPlethysm with args plethysmer selector, order, plethysmed one selector, order
Selectors:
sym = {order}, sym all = {p} up to {order}
ats = {1^order}, ats all = {1^p} up to {1^order}
gen = (all YD's with size = order), gen all = (all YD's with size <= order)

As far as I can tell, it works.

#### lpetrich

conjg(plt(a,x)) = plt(a,conjg(x)) is correct for conjugate representations, but for conjugate Schur functions / Young diagrams, one finds this:

conjg(plt(a,x)) =
plt(a,conjg(x)) if len(x) is even
plt(conjg(a),conjg(x)) if len(x) is odd

len(x) is the length or order or total power of x.

Another identity: plt(a,plt(b,x)) = plt(plt(a,b),x)

Some special cases: symmetric and antisymmetric, {n} and {1^n}. I call them that because of the corresponding SU(N) / SL(N) irreps and how they are composed from the fundamental one: the symmetrized or the antisymmetrized product.

plt({n},x+y) = sum over k of plt({n-k},x)*plt({k},y)
plt({1^n},x+y) = sum over k of plt({1^(n-k)},x)*plt({1^k},y)

plt({n},x*y) = sum over a of plt(a,x)*plt(a,y) where len(a) = n
plt({1^n},x*y) = sum over a of plt(a,x)*plt(conjg(a),y) where len(a) = n

#### lpetrich

I think that I ought to give the algorithm that I conjecture for plethysms of Lie-algebra representations. I start with a rep given by a list of sets (n,v) where n is the multiplicity/degeneracy and v is a root vector. The root vectors are all distinct in the list.

For power p, the results are lists of multiplicity-vector sets (N(ni,1,ni,2,...ni,p), vi,1 + vi,2 + ... + vi,p)
I then find irreps from it in standard fashion, by repeatedly looking for the highest weight and subtracting out the irrep with that highest weight.

The problem is, of course, finding the N functions. For the pure symmetric and antisymmetric cases, it is easy. First, find the multiplicities m of each n value, a list of
(mk,nk)

and then one finds
N = product over the (m,n)'s of n*(n+s)*(n+2*s)...(n+(m-1)*s)/m!
where s = 1 for symmetric and -1 for antisymmetric. It's rather easy to find this result from symmetry.

But the mixed-symmetry cases are much more difficult. They are plethysms of Young diagrams y on the rep, and I calculate them for all the y's that add up to the power. The y's correspond to irreps of SU(n) that one finds from powers of SU(n)'s fundamental rep, and a plethysm of y on the fundamental rep yields the rep with Young-diagram y.

That's rather easy to show from the correspondence between Schur functions and the characters of SU(n) Lie algebras. In general, a Lie-algebra character for a rep is the sum of n*exp(c.v) for vector c and rep pairs (n,v).

So I calculate all the N(y,{n}) for Young diagrams y and vectors of n: {n}.

-

I first calculate NW(yw,{n}) the number for the SU(n) Weyl orbit with max weight yw, then find
N(y,{n}) = sum over yw of Kostka(y,yw)*NW(yw,{n})

where

(SU(n) irrep for y) = sum over yw of Kostka(y,yw) * (SU(n) Weyl orbit for yw)

or in Schur functions,

S(y,x) = sum over yw of Kostka(y,yw) * SW(yw,x)

where SW(yw,x) is a Schur-Weyl function, the sum over all terms xi1r1*xi2r2*...*xi(n)r(n) for yw = {r1,r2,...,r(n)} and selections xi's from the vector x.

The Kostka matrix I find using the algorithm in "Determinantal Expression and Recursion for Jack Polynomials", by L. Lapointe, A. Lascoux, J. Morse. Are there any recursive algorithms for this matrix? Algorithms that calculate each power's Kostka matrix from the previous ones.

-

NW(yw,{n}) turns out to be relatively simple. From {n}, find {(m,n)}, and sort the list to place the m's in non-increasing order. The m's form the Young diagram ym.

Find all nestings of yw in ym. Row lengths: yw = {rw}, ym = {rm}. These are all possibilities of
rm(i) = sum over k of rw(j(i,k))
where each rw value appears once and only once in the sums and nestings with different arrangements of repeated values are counted as a single nesting.

For each nesting, find the product over the rm(i)'s in ym of:
using power = count of rw(j(i,k))'s
using repeat counts = counts of repeats of rw(j(i,k)) values.
power!/(product of count! for count in repeat counts) * n(i)*(n(i)-1)*...*(n(i)-power+1)/power!

Add up over all nestings. No nestings possible yields zero. That gives NW(yw,{n}). Find for all yw's. Multiply by the Kostka matrix for the overall power to get the N(y,{n})'s.

This algorithm essentially calculates the plethysm of SU(n) Weyl orbits on the rep character, then uses the appropriate Kostka matrix to find the plethysm of SU(n) irreps on the rep character.

This suggests a proof: turn the rep character into a sum of monomials exp(c.v(i)) where each such term is repeated n(i) times. Then make it a vector of exp(c.v(i)) values and use it as the vector arg in the Schur or Schur-Weyl function. I'd have to work out a lot of combinatorics, but it seems like it could be done in the Schur-Weyl case.

#### Simon_Tyler

Hey lpetrich,

I don't have much to say about your recent posts and work, since it's out of my range of experience and I don't have time to study it at the moment.
So I'll just say that it sounds good and keep going!

#### lpetrich

Thanx.

I worked out the combinatorics for the Schur-Weyl case. Let's consider the case where x1, x2, ..., x(n) have the same root value v, the same value of exp(c.v). Let's also consider them divided among the Schur-Weyl row lengths r1, r2, ..., r(p): x1r1*x2r2*...*x(p)r(p) and all permutations. There are thus n*(n-1)*...*(n-p+1) such terms. However, sets of equal r's have the permutations of x's in them counted only once, which means dividing the number of terms by the product (number of r equal to some value)! for all values that the r's equal.

If the full set of terms contains something like x1*x2*x3 + x1*x3*x2 + x2*x1*x3 + x2*x3*x1 + x3*x1*x2 + x3*x2*x1, one must divide by (# of x's here)! = 3! to count only x1*x2*x3.

Which is the final step of the proof.

I note as an aside the Weyl character formula for SU(2) irreps. Using quantum-mechanical notation with angular momentum j, the highest weight becomes {2j}, and the character

X(j,c) = exp((-j/2)*c) + exp((-j/2+1)*c) + ... + exp((j/2)*c)
= sinh((j+1)*c/2)/sinh(c/2)

#### lpetrich

Not surprisingly, large reps take a long time to process. I tried the square plethysm of the fundamental irrep of E8, the "248", and I found about 3 minutes for Mathematica and 1 minute for Python on my 2-GHz Intel Core 2 Duo iMac

starting
(0, 0, 0, 0, 0, 0, 1, 0): 248
symmetric part of square
1 of (0, 0, 0, 0, 0, 0, 2, 0): 27000
1 of (1, 0, 0, 0, 0, 0, 0, 0): 3875
1 of (0, 0, 0, 0, 0, 0, 0, 0): 1
antisymmetric part of square
1 of (0, 0, 0, 0, 0, 1, 0, 0): 30380
1 of (0, 0, 0, 0, 0, 0, 1, 0): 248
Using Robert Cahn's root-numbering convention

I've been considering alternatives to both Mma and Python in the hope of getting a speedup. They have varying advantages and disadvantages.

Here are some benchmarking sites:
http://shootout.alioth.debian.org/u64/which-programming-languages-are-fastest.php [Broken], the most useful that I've found
Programming Languages Benchmarks
The winners are Fortran, C, C++, Ada, ATS, and Java.

Neither Ada nor ATS is very well-known, so I'll ignore them. The next question is whether the languages have convenience features for some of the more common tasks. Associative arrays are good for finding whether a vector is already in a list of vectors, a very common task in my code, and operator overloading is good for implementing fractions and abstract-algebra entities.

Programming Language Comparison has a section on operator overloading, and of the fast languages, only C++ has it. Looking in Wikipedia, I find that of those languages, only C++ and Java have them.

So I choose C++. With its templates and operator overloading, one can implement fractions, and with its Standard Template Library template "map", one can implement associative arrays.

Last edited by a moderator:

#### Simon_Tyler

Have you tried just using compilation to C in both Mathematica and Python?
http://reference.wolfram.com/mathematica/Compile/tutorial/CompilationTarget.html
http://docs.cython.org/src/quickstart/overview.html
The generated C code can sometimes be almost as fast or faster than a hand coded C/C++
http://fredrik-j.blogspot.com/2009/05/report-from-sage-days-15.html <-- this blog is worth reading
http://carlo-hamalainen.net/blog/2008/11/09/cython-vs-c-improved/ [Broken]
It also makes parallelization easier:
http://www.walkingrandomly.com/?p=3635
http://reference.wolfram.com/mathematica/Compile/tutorial/Efficiency.html

The advantage of using these is that you only need to compile the key steps in your code and only need to make minimal modifications to your existing code.

The Mathematica Compile[..., CompilationTarget->"C"] function will fall back to normal virtual machine compilation if no C compiler exists.
And for Cython, you can use python decorators to achieve a similar thing
http://docs.cython.org/src/tutorial/pure.html
This is even easier if you're using sage
http://docs.cython.org/src/quickstart/build.html?highlight=sage

The other thing to try is to use (integer) packed arrays and make sure that you don't accidentally unpack them.
For Mathematica, see
http://www.wolfram.com/technology/guide/PackedArrays/
http://stackoverflow.com/questions/4721171/performance-tuning-in-mathematica
For pure Python see
http://docs.python.org/library/array.html
but maybe it would be easier and faster to use Numpy
http://docs.scipy.org/doc/numpy/user/basics.types.html
http://stackoverflow.com/questions/8385602/why-are-numpy-arrays-so-fast

Of course, if you do completely rewrite your code in C or C++ then you can create interfaces for Mathematica and Python (and most other languages). However, your code will also become less accessible to most people.

Anyway, sorry for another link dump. Hopefully you find something useful in it!

Last edited by a moderator:

#### lpetrich

Thanx for those tips. I had no idea that the Mma Append / AppendTo / Prepend / PrependTo functions would get computationally expensive for long lists -- those functions make a new copy of the list.

[noparse]So I've replaced them with indexing functions. Instead of PosRoots[[ix]], I do PosRoots[ix], and maintain a total count. I don't really need to do so, since I can get all the definitions with DownValues[PosRoots].[/noparse]

Python: 71 s
Mathematica: 74 s
About 2.5 times faster than before, and neck-and-neck with Python.

I've updated my archive with the new version. Also, for Notable Physics Results, I've added a couple results to the the SU(3) results in the Flavor and Gauge Symmetry section. A summary:
Code:
(* Create SU(2) and SU(3) Lie algebras: spin, flavor, color *)
MakeLieAlgebra[spnsu2, {1, 1}]
MakeLieAlgebra[gfssu3, {1, 2}]

(* Spin and flavor together -- all symmetry types *)
DecomposeAlgProdRepPower[{spnsu2, gfssu3}, {{1}, {1, 0}}, 3]

(* Spin, flavor, and color together -- antisymmetry for Fermi-Dirac statistics *)
DecomposeAlgProdRepPwrSym[{spnsu2, gfssu3, gfssu3}, {{1}, {1, 0}, {1, 0}}, 3, -1]

#### Simon_Tyler

Yeah, Leonid's performance tuning tips are good (like all of his answers - and his book!).

I've sent you an invite for the new Mathematica.StackExchange site, which is currently in private beta. If you have any questions about optimizing aspects of your code, this is currently the place to ask them (since everyone is a little overenthusiastic at the moment).

Oh, and +1 for the Notable Physics Results!

#### lpetrich

Thanx.

I updated the Mathematica and Python versions with some reorganization and generalizaiton of the code. The rep products and the branching can now accept reducible reps as inputs, specified as lists of irreducible ones. That makes them like the plethysms.

For a C++ version, I have considered SWIG: Simplified Wrapper and Interface Generator, enabling one to build add-ons written in C/C++ for several programming languages, including Python. However, I've found it to have a rather serious deficiency: it does not transfer arrays back and forth. One would have to do member-by-member transfers.

Another curiosity I've come across is a use of G2 in a context that does not suggest exceptional Lie algebras. Back in 1949, Giulio Racah was trying to analyze f electron states in atoms, and he used the subalgebra chain
U(7) / SU(7) > SO(7) > G2 > SU(2) / SO(3)

With my code, one sets up
SubalgSUSO[ld1,7]
SubalgExtra[ld2,"B3G2"]
SubalgExtra[ld3,"G2A1"]

By calculating the overall projection matrix, I've found that this chain is equivalent to
SubalgHeightA1[ld123,{1,6}]
For each root vector, add the root-vector components and it becomes the root-vector component for SU(2).

The G2->A1 one is equivalent to a related branching-rule generator
SubalgHeightA1[ld3x,{7,2}]

#### lpetrich

I've been working on a C++ version, slowly but surely. I've created lots of template functions for doing the vector and matrix operations, and also template classes for matrices and fractions. I've done that to make them as general as possible. I've also been using the Standard Template Library's classes and functions, because it handles all its memory allocation and deallocation behind the scenes, and because I can create associative arrays without much trouble. I use that to check if some calculated root vector is already present in a list of them.

So far, I've done the algebra-data construction, and the irreducible-representation total degeneracies. I had to use the GNU Multiple Precision library to do bignums because one of my test cases overflowed 64-bit integers during the calculation. 64 bits is the largest size supported by recent C/C++ compilers. Here it is:

Algebra: E8, with Robert Cahn's root order
Max Wts: 0 0 2 0 0 0 0 0
Total Degen: 684252595649750400
It needs only 60 bits, but an intermediate result went past 63 bits, since I was using signed ones.

#### lpetrich

I've got my C++ version working up to where it finds representations as lists of root/weight vectors. Here's C++ vs. Python vs. Mathematica on my 2-GHz Core 2 Duo iMac:

E8
00000010 -- 248
0.056 - 0.406 - 0.466
10000000 -- 3875
0.894 - 5.017 - 5.875
00000020 -- 27000
4.436 - 24.575 - 26.904
00000100 -- 30380
4.232 - 23.876 - 25.963
00000001 -- 147250
14.892 - 82.709 - 86.085

About 6 times faster, while Mma is slightly slower than Python.

#### lpetrich

First, the Python ones were with out-of-the-box Python, without numerical-Python add-ons.

I was compiling without optimization. Let's see what it does with the optimization levels available with gcc: 0, 1, 2, 3.

I've also changed from <sys.time.h> gettimeofday() to <time.h> clock(), for a better measure of CPU time.

Mathematica, Python, C++ 0, C++ 1, C++ 2, C++ 3
Algebra: E8
00000010 -- 248
0.466 - 0.406 - 0.054 - 0.012 - 0.010 - 0.007
10000000 -- 3875
5.875 - 5.017 - 0.811 - 0.183 - 0.163 - 0.113
00000020 -- 27000
26.904 - 24.575 - 3.983 - 0.926 - 0.821 - 0.569
00000100 -- 30380
25.963 - 23.876 - 3.821 - 0.894 - 0.793 - 0.554
00000001 -- 147250
86.086 - 82.709 - 13.308 - 3.219 - 2.857 - 2.022

Speedups relative to the Mma version (last 3 reps):
Python: 1.04 - 1.09
C++ 0: 6.5 - 6.8
C++ 1: 27 - 29
C++ 2: 30 - 33
C++ 3: 43 - 47

#### lpetrich

I'm now getting best-case numbers like 1.6 seconds for E8 00000001 -- ~54 times faster than Mma and Python. It certainly helps to do much of the object management and method dispatch at compile time.

I've updated my archive with my C++ source code and a makefile for a shared library and some test programs. The shared library is what my main code gets compiled into. I have rep powers / plethysms working in it, complicated coding and all, though not branching rules.

I've half-thought of doing a Java version, but since Java does not support operator overloading, a fraction class would be awkward, and I'd have to integerize the matrix inversion. I'm not familiar enough with Perl or Ruby to want to work with those. For a JavaScript demo in a web browser, there's a risk of causing a user's web browser to hang or crash from CPU and memory usage.

#### lpetrich

I'm still in the planning stages for branching rules. The main problem I face is how to code for the special cases in C++, since I can't use the sorts of array values that I can use in Mathematica and Python. I'm thinking of storing the exceptional-algebra special cases in a serialized fashion:

number of algebras, for each subalgebra:
- family, rank, projection matrix expanded as:
- - values in rows and columns, expanded as:
- - - numerator, denominator
This should be straightforward to generate from the Mathematica version.

I could even do all the special cases together in this fashion, giving one big array for the root demotions, one for the extension splittings, and one for the extras. The extras I might specify with enums: SUBALG_EXTRA_B3G2, etc. Again straightforward to generate from the Mathematica version.