SetDirectory[NotebookDirectory[]];
<<RISC`HolonomicFunctions`
<<"ZonalPolynomials.m";

HolonomicFunctions Package version 1.7.3 (21-Mar-2017)
written by Christoph Koutschan
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Calculate zonal polynomials

When we only give a partition, we obtain the zonal polynomial in terms of the symmetric monomial functions:

ZonalPolynomial[{3,2}]

Alternatively, one can give a list of variables as second argument to see the polynomial explicitly:

ZonalPolynomial[{2,1},{a,b,c}]

The coefficients (Subscript[c, κ,λ]) of the zonal polynomials are computed recursively, by the following command:

ZonalCoefficient[{2,1},{1,1,1}]

Let’s compute a zonal coefficient for larger partitions:

ZonalCoefficient[{8,6,6,3},{7,7,5,3,1}]//Timing

We can create a table of all zonal polynomial coefficients indexed by partitions of n=4:

With[{p=IntegerPartitions[4]},Outer[ZonalCoefficient,p,p,1]]//TableForm

For a simpler (and slightly faster) way to create this table, there is the following command:

ZonalCoefficientTable[4]//TableForm

For the zonal coefficients (Subscript[c, κ,λ]) in the upper left corner, we can derive general closed forms for arbitrary n. This refers to the situation when κ and λ are of the form (n-i,π) where π is a partition of i and where n is symbolic. Here n is assumed to be sufficiently large, so that n-i is greater than or equal to the largest part of π).

ZonalCoefficientN[{n-3,2,1},{n-4,2,2}]

This is the table that appears in the paper:

part=Flatten[Table[Prepend[#,n-i]&/@IntegerPartitions[i],{i,0,3}],1];
tab=Outer[ZonalCoefficientN,part,part,1];
trad=TraditionalForm/@part;
TableForm[Join[{Prepend[trad,"κ\λ"],{}},Join[Transpose[{trad}],tab,2]]]

Experimentally verify definition with Wishart distribution

Mathematica provides a command for sampling random matrices according to the Wishart distribution.
Here we show the distribution of the largest eigenvalue of a 2*2 Wishart matrix with parameter ν = 3.

Histogram[Table[Max[Eigenvalues[RandomVariate[WishartMatrixDistribution[3,IdentityMatrix[2]]]]],{10000}]]

However, it is also possible (and easy) to define the Wishart distribution in terms of the normal distribution:

Histogram[Table[Max[Eigenvalues[With[{X=Table[RandomVariate[NormalDistribution[0,1]],{3},{2}]},Transpose[X].X]]],{10000}]]

(* Symmetric functions U_lambda. They form a basis. *)
SymU[λ_,y_]:=If[Length[λ]>Length[y],0,With[{lam=Append[λ,0]},Product[SymmetricPolynomial[i,y]^(lam[[i]]-lam[[i+1]]),{i,Length[λ]}]]];
BasisU[n_,y_]:=DeleteCases[SymU[#,y]&/@IntegerPartitions[n],0];
BasisU[3,{y1,y2}]

(* The result of the map tau, when we use only a single randomly generated matrix. *)
τ[ν_,U_]:=
Module[{y=Variables[U],W,ev},
W=RandomVariate[WishartMatrixDistribution[ν,IdentityMatrix[Length[y]]]];
W=Round[W*10^6]/10^6;
ev=Eigenvalues[DiagonalMatrix[y].W];
Return[N[Expand[U/.Thread[y->ev]]]];
];
τ[3,SymU[{3,1},{y1,y2}]]

(* Compare exact result and Monte-Carlo simulation for two variables. *)
Test2[n_,ν_,s_:10000]:=
Module[{parts,matXi,Lam,eqns,sim},
parts=Select[IntegerPartitions[n],Length[#]≤2&];
matXi=Table[xi[i,j],{i,Length[parts]},{j,Length[parts]}];
eqns=matXi.BasisU[n,{y1,y2}]-(With[{k=Length[#]},Product[(2*#[[i]]+k-i)!,{i,k}]*(2n)!/(2^n*n!*Product[2*#[[i]]-2*#[[j]]-i+j,{i,k-1},{j,i+1,k}])]*ZonalPolynomial[#,{y1,y2}]&)/@parts;
matXi=matXi/.First[Solve[Thread[Flatten[CoefficientList[#,{y1,y2}]&/@eqns]==0],Flatten[matXi]]];
Lam=DiagonalMatrix[2^n*Product[Pochhammer[(ν+1-i)/2,#[[i]]],{i,Length[#]}]&/@parts];
Print["Exact: ",Expand[Inverse[matXi].Lam.matXi.BasisU[n,{y1,y2}]]];
Print["Simul: ",sim=Expand[Sum[τ[ν,BasisU[n,{y1,y2}]],{s}]/s]];
Return[sim];
]

Test2[2,2];

Test2[2,5];

Test2[3,4];

Timing[Test2[4,3,10^5];]

Zonal polynomials are the eigenfunctions of the Laplace-Beltrami operator

We do an experimental verification that the zonal polynomials are the eigenfunctions of the Laplace-Beltrami operator.

LaplaceBeltrami[{x,y,z}]

Table[
With[{y=Array[x,n]},
Function[λ,Together[
ApplyOreOperator[LaplaceBeltrami[y],ZonalPolynomial[λ,y]]/
((ρ[λ]+(Length[y]-1)*Total[λ])*ZonalPolynomial[λ,y])
]]/@Select[Join@@Table[IntegerPartitions[n],{n,2,10}],Length[#]≤Length[y]&]
],{n,5}]

Proofs

Proof of Proposition 6.2

fj=(b-2j)/d/(2b-2d+1)*Binomial[b,j]*Pochhammer[1/2,j]/Pochhammer[b-j+1/2,j]

ct=Factor[CreativeTelescoping[fj,S[j]-1,{S[b],S[d]}]]

gj=Simplify[ApplyOreOperator[-ct[[2,1]],fj]]

(* Sanity check: verify the identity g(j+1)-g(j)=f(j) *)
FullSimplify[((gj/.j->j+1)-gj)/fj]

gj/.{{j->0},{j->d}}

Proof of Theorem 6.3

(* closed form given in Theorem 6.3 *)
cf[a_,b_,d_]:=Binomial[b,d]*(b+1/2)*(2a-b)!/(a-b)!/b!*Pochhammer[1/2,d]/Pochhammer[b-d+1/2,a-b+d+1];
cf[a_,b_,d_]:=(b+1/2)*(2a-b)!*Pochhammer[1/2,d]/d!/(b-d)!/(a-b)!/Pochhammer[b-d+1/2,a-b+d+1];
cf[a,b,d]

Clear[ff]
SetDelayed@@{ff[a_,b_,d_],Simplify[cf[a+d,b+2d,d]]/((2a-b)!/a!/(a-b)!)}
ff[a,b,d]

ct=CreativeTelescoping[ff[a,b,d],S[d]-1,{S[b],S[a]}]//Factor

g1=ApplyOreOperator[-ct[[2,1]],ff[a,b,d]]
g2=ApplyOreOperator[-ct[[2,2]],ff[a,b,d]]

(* Some manual simplifications on the g's (formulas from the paper) *)
Clear[g1,g2];
g1[d_]:=-((a! (a-b)! Pochhammer[(1/2),d])/((d-1)! (b+d-1)!(a-b-d+1)! Pochhammer[(1/2)+b+d,1+a-b]));
g2[d_]:=(a! (a-b-1)! Pochhammer[(1/2),d])/((d-1)! (b+d)!(a-b-d)! Pochhammer[(1/2)+b+d,a-b]);
(* Sanity check: verify identities (6.6) and (6.7). *)
FullSimplify[(ff[a+1,b,d]-ff[a,b,d])/(g1[d+1]-g1[d])]
FullSimplify[(ff[a,b+1,d]-ff[a,b,d])/(g2[d+1]-g2[d])]

g1[a-b+1]-g1[0]+ff[a+1,b,a-b+1]

Together[%/.Pochhammer[3/2+a,a-b+2]->(5/2+2a-b)*Pochhammer[3/2+a,a-b+1]]

g2[a-b+1]-g2[0]-ff[a,b+1,a-b]