Rotating Homogeneous Polynomials & Representations of SO(3)

Under the guise of rotating solid spherical harmonic functions, I want to investigate the computational aspects of representing rotations of homogeneous polynomials. The connection being that harmonic homogeneous polynomials are solid spherical harmonics. My goal in this note is to rotate solid homogeneous polynomials ( in variables x, y, and z ) by computing representations of SO(3) on Pd, homogeneous polynomials of degree d. One nice aspect of computing representations of SO(3) is that only algebra is required. Differential equations typically used to describe/derive spherical harmonics make no appearance in what follows.

Go straight to the code if your primary interest is the implementation of rotating homogeneous polynomials.

Preliminaries

The type of homogeneous polynomials I want to focus on are functions of the form

Pd={ f:RR | f(rx)=rdf(x) rR}.

SO(3), the group of 3dimensional rotations, acts on Pd by rotating coordinates of the domain R3. Suppose pPd is a homogeneous polynomial then the action of ASO(3) is

Ap=p(Ax).

The result of the action is just another homogeneous polynomial of degree d, because rotations preserve distance.

Ap(rx)=p(Arx)=p(rAx)=rdp(Ax).

Let A, BSO(3) and denote the action instead by ρ. The action is actually a homomorphism.

ρAB(p(x))=p(ABx)=ρA(p(Bx))=(ρAρB)(p(x)).

Individually, each ρA is linear.

ρA(p+q)=(p+q)(Ax)=p(Ax)+q(Ax)=ρA(p)+ρA(q)ρA(λp)=λp(Ax)=λρA(p).

Additionally, ASO(3) is a change of coordinates (and therefore a linear isomorphism of R3) implies each ρA has trivial kernel. To see this let x~=Ax denote the new coordinates of R3. Then 0=p(x~) x~R3 directly implies that p=0. With that ρd is a representation of SO(3) on Pd. Specifically, ρd is a homomorphism

ρd:SO(3)GL(Pd).

What I did recently was compute the eigenvalues of z-axis rotations on P2. All z-axis rotations have the following form

Aθ=(gθ001)SO(3)where gθ=(cosθsinθsinθcosθ)SO(2).

For notational simplicity I prefer to use a and b instead of cosθ and sinθ to describe gθ.

gθ=(abba)where a2+b2=1.

First consider representations of Aθ on P1. The basis functions {x, y, z} are just the coordinates functions of R3, so computing the rotated basis of P1, in terms of {x, y, z}, requires simply applying Aθ.

Aθx=ax+by Aθy=bx+ay Aθz=z

So for polynomial c1x+c2y+c3zP1 which in vector form is [c1c2c3]

Aθc1x+c2y+c3z=[xyz](gθ001)[c1c2c3].

Restricting to the vector representation of P1 one sees that ρAθ=Aθ. Now, that the form of ρ is known for z-axis rotations on P1 the eigenvalues of the representation are easily computed. They are namely {eiθ,1,eiθ}.

To make the calculation of the representation of Aθ on P2 easier I choose the basis to be

{x2y2, xy, x2+y2,zx, zy, z2zP1}.

It’s immediately apparent that the last three basis elements are products of z and the basis of P1. Since the action of Aθ on z is the identity, the representation of Aθ on zP1 is the same as on P1. Rewriting the P2 basis with the zP1 subspace

{x2y2, xy, x2+y2}zP1.

What might be less apparent is that the complement subspace is closed under the action of Aθ.

Aθx2y2=(ax+by)2(bx+ay)2=a2x2+2abxy+b2y2(a2y22abxy+b2x2)=(a2b2)(x2y2)+4ab xy Aθxy=(ax+by)(bx+ay)=ab(x2y2)+(a2b2) xy

Since x2+y2 is unchanged by zaxis rotations, Aθx2+y2=x2+y2.

Therefore, the representation of Aθ on P2 is

ρ(Aθ)=[a2b2ab04aba2b200001ab00ba0001]

After substituting a=cosθ and b=sinθ, the eigenvalues of ρ(Aθ) are easily computed, to be {ei2θ,eiθ,1,eiθ,ei2θ}. All the eigenvalues except for 1 have multiplicity one; eigenvalue 1 has multiplicity two.

Nice…but computing the representation of Aθ on P1 and P2 as I have above is not scalable at all, because there are so many terms to multiply. What I want is to express a SO(3) representation in Pd for any d as a product of matrices involving, necessarily, Aθ.

Computing the Representation of Aθ on Pd

While developing the formula for ρ2, shown at the end of the previous section, the upper summand of ρ2

[a2b2ab04aba2b20001]

looked to me like a sort of square of Aθ. Indeed the portion of P2 transformed by it are the quadratic polynomials in x and y, and the spectrum of the upper summand is a square of all the eigenvalues of Aθ. ρ2 is not the square of Aθ, but the idea of somehow realizing ρ2 as a product of A with itself stuck with me. But what type of product? What properties should it preserve?

Well, the product should mimic the action of A (for brevity I am dropping the subscript θ) on the homogeneous basis polynomials…clearly. Consider

Axy=(Ax)(Ay)=(ax+by)(bx+ay)=ab x2+a2 xyb2 yx+ab y2.

A distributes to x and y in the first line. The actions of A on x and y, linear homogeneous polynomials, are resolved, resulting in a factored expression, which is expanded in the last line. The expansion is determined by the rules of multiplication: multiplication distributes over addition. So, two things are happening here. First the action of A on xy is actually a product of A-actions on lower degree homogeneous polynomials, here x and y. Think recursion. Second, expanding the product of the actions on the lower degree factors expresses the result in the basis of P2 and relies only on the bilinearity of multiplication. Cutting to the chase…the tensor product of linear operators does both these things.

A can be tensored with itself, which is a well defined extension of A to P1P1. Identifying for the moment xy with xy – more on embedding P1 into tensor space P1P1 later – notice that the action of AA on xy is very close in appearance to the action of A on xy.

AAxy=AxAy=(ax+by)(bx+ay)=ab xx+a2 xyb2 yx+ab yy.

AA, in the first line above, factors into actions of A on lower degree homogeneous polynomials in each tensor coordinate: Ax on the left side of the tensor, and Ay on the right of . The action of A on the lower degree factors is no different than before, and neither is the expansion of the factored expression in line 2 because is bilinear. What remains is to collect like terms and project P1P1 into P2, arriving at a finally expression of the action on xy in terms of P2.

Now a approach to computing ρd can be stated, and it is simply

ρ2=S2(AA)E2

and for any degree d the recursive formula

ρd=Sd(Aρd1)Ed.

Ed:PdP1Pd1 and Sd:P1Pd1Pd are the embedding and projection operators, respectively. I will describe their constructions shortly in the next sections. The key though to scalably computing ρ2 is the matrix tensor product. It is a pure matrix operation so it is easily computed/programmed.

Ordering the Basis of Pd

The main function of Sd, the projection operator introduced in the recursive formula at the end of the previous section, is to collect the contributions of like terms. Projecting P1Pd1 into Pd accomplishes this. But what are like terms? Of main interest here are the basis elements of P1Pd1. Referring to the action of AA on xy, it seems sensible that xy and yx should be considered alike.

Let basis tensors be equivalent if the sum of degrees in each variable x, y, z over both tensor coordinates are the same. For clarity represent the degree sums in coordinates (degx,degy,degz). A basis tensor of P1Pd1 is simply a factorization of some degree d homogeneous monomial into degree 1 and degree d1 monomials; and all I am saying is that basis tensors are alike if they are factorizations of the same degree d homogeneous monomial. Consider three basis tensors of P1Pd1:

xxa1ybzc(1+a1,b,c)=(a,b,c)yxayb1zc(a,1+b1,c)=(a,b,c)zxaybzc1(a,b,1+c1)=(a,b,c).

All have the same degree sums coordinate; so xxa1ybzc, yxayb1zc, and zxaybzc1 are all equivalent. With an equivalence in hand a minimal definition of Sd can be stated: Sd maps tensors into their equivalence classes; specifically, Sd maps P1Pd1 into P1Pd1/ like so

Sd(xxa1ybzc)=[zxaybzc1].

Projecting tensors into their equivalence classes is a linear map which has the effect of collecting the contributions of equivalent terms. Take for instance a linear combination of the equivalent terms, then

Sd(p1xxa1ybzc+p2yxayb1zc+p3zxaybzc1)=(p1+p2+p3)[zxaybzc1].

The quotient space P1Pd1/ is actually isomorphic to Pd via the degree sums coordinate mapping. The degree sums coordinate, (a,b,c), represents a unique homogeneous monomial xaybzcPd; so projecting a linear combination of equivalent terms is also

Sd(p1xxa1ybzc+p2yxayb1zc+p3zxaybzc1)=(p1+p2+p3)xaybzc.

The embedding operator is a sort of inverse of projection, and actually with appropriate bases ordering it is a right inverse of the projection. It doesn’t have to be though; here is a minimal defintion: it takes a given homogeneous monomial in Pd to a distinguished basis tensor among those in P1Pd1 which share the same degree sums coordinate as it. This distinguished basis tensor is an equivalence class representative of the basis tensors that project to given degree d homogeneous monomial.

To realize the embedding and projection operators as matrices requires fixing the equivalence class representatives of P1Pd1, and ordering both P1Pd1 and Pd. I could be idiosyncratic and choose class representatives by making some random choices; however that wouldn’t necessarily be scalable or easily scalable. Because recursion would then require storing or preserving the representatives of lower degree objects. Anyway, it’s way easier to just fix an order on Pd, extend it to P1Pd1, and then choose the class representatives in the quotient object based on said order. The later is most important for creating the matrix form of the embedding; since Sd collects contributions of all terms in an equivalence class, all that matters is the ordering of the bases.

The convention I use in plain words is to abstract d-degree homogeneous monomials into a coordinate representation (degx,degy,degz), and then order the monomials first by degz and then degy. For example the basis of P2 in ascending order is

x2(2,0,0)xy(1,1,0)y2(0,2,0)xz(1,0,1)yz(0,1,1)z2(0,0,2).

One of the benefits of the ordering is the basis of subspace zP1, shown above in yellow, is ordered consecutively at the end; in general the pure homogeneous polynomials of degree d in x and y are ordered before zPd1. A direct consequence is that ρd has a nice clear direct sum of matrices structure. P1Pd1 is ordered by coordinate, left to right; below the basis of P1P2 is shown in ascending order.

xx2(3,0,0)xxy(2,1,0)xy2(1,2,0)xxz(2,0,1)xyz(1,1,1)xz2(1,0,2)yx2(2,1,0)yxy(1,2,0)yy2(0,3,0)yxz(1,1,1)yyz(0,2,1)yz2(0,1,2)zx2(2,0,1)zxy(1,1,1)zy2(0,2,1)zxz(1,0,2)zyz(0,1,2)zz2(0,0,3)..

The degree sum coordinate of each basis tensor is written into the column on the right, and these sums aid in identifying equivalent basis tensors. For instance I have highlighted in orange every occurrence of (1,1,1). The corresponding basis tensors are underlined or boxed and belong to the same equivalence class. The boxed tensor, zxy, is maximum element of the class. There are other equivalence classes, and each has a maximum element with respect to the defined ordering. This holds in general, so I define the equivalence class representative to be the maximum tensor among those in the class.

Projecting P1Pd1 into Pd AND Embedding Pd in P1Pd1

The embedding and projection matrices need ordered bases to be defined and constructed. Invoking the just described orderings of P1Pd1 and Pd in the previous section, create a table where:

Since there is a clear bijection between homogeneous monomials in Pd and their degree sum coordinates, the projection Sd is completely determined by the degree sums coordinate mapping. To then determine the set of equivalent basis tensors of P1Pd1, each of which labels a unique row, that map to a given homogeneous monomial of in Pd, each of which its own column, one simply matches the degree sum coordinate of the degree d homogeneous monomial to those in its column. To illustrate I made a table for d=2 and in each column highlighted repeated occurances of the corresponding degree coordinate.

x2xyy2xzyzz2xx(2,0,0)(2,0,0)(2,0,0)(2,0,0)(2,0,0)(2,0,0)xy(1,1,0)(1,1,0)(1,1,0)(1,1,0)(1,1,0)(1,1,0)xz(1,0,1)(1,0,1)(1,0,1)(1,0,1)(1,0,1)(1,0,1)yx(1,1,0)(1,1,0)(1,1,0)(1,1,0)(1,1,0)(1,1,0)yy(0,2,0)(0,2,0)(0,2,0)(0,2,0)(0,2,0)(0,2,0)yz(0,1,1)(0,1,1)(0,1,1)(0,1,1)(0,1,1)(0,1,1)zx(1,0,1)(1,0,1)(1,0,1)(1,0,1)(1,0,1)(1,0,1)zy(0,1,1)(0,1,1)(0,1,1)(0,1,1)(0,1,1)(0,1,1)zz(0,0,2)(0,0,2)(0,0,2)(0,0,2)(0,0,2)(0,0,2)

Boxed and highlighted cells indicate the equivalence class representative–the class maximum. Now the matrix Sd can be constructed from the table by converting table cells into scalars.

Finally, transpose the result to get the matrix form of the projection P1P1 into P2.

#Project tensor space into space of homogeneous polys.
S = Matrix([1, 0, 0, 0, 0, 0, 0, 0, 0],
           [0, 1, 0, 1, 0, 0, 0, 0, 0],
           [0, 0, 0, 0, 1, 0, 0, 0, 0],
           [0, 0, 1, 0, 0, 0, 1, 0, 0],
           [0, 0, 0, 0, 0, 1, 0, 1, 0],
           [0, 0, 0, 0, 0, 0, 0, 0, 1]]).

Likewise, the embedding matrix Ed can be constructed from the table by converting table cells into scalars.

The resulting matrix is the embedding P2P1P1.

#Embed into tensor space.
E = Matrix([
	[1, 0, 0, 0, 0, 0],
	[0, 0, 0, 0, 0, 0],
	[0, 0, 0, 0, 0, 0],
	[0, 1, 0, 0, 0, 0],
	[0, 0, 1, 0, 0, 0],
	[0, 0, 0, 0, 0, 0],
	[0, 0, 0, 1, 0, 0],
	[0, 0, 0, 0, 1, 0],
	[0, 0, 0, 0, 0, 1]]).

Putting it all together, the representation Aθ on P2 is

In [74]: repN = S * TensorProduct(A,A) * E

In [75]: repN
Out[75]:
Matrix([
[ a**2,        -a*b,   b**2, 0,  0, 0],
[2*a*b, a**2 - b**2, -2*a*b, 0,  0, 0],
[ b**2,         a*b,   a**2, 0,  0, 0],
[    0,           0,      0, a, -b, 0],
[    0,           0,      0, b,  a, 0],
[    0,           0,      0, 0,  0, 1]])

The degree 2 representation computed here looks a bit different than the one computed earlier at the end of the preliminaries section because the bases are slightly different; otherwise it is correct. The code I wrote can generate the matrix representation ρd for any degree d. For instance d=3

And because of the ordering ρd has a nice direct sum of matrices structure…for z-axis rotations. The method and code I have described correctly computes representations for any 3D rotation, however no guarantees the representation will decompose nicely, as it does for z-axis rotations.

Eigenvalues of ρd

Earlier drafts of this post were actually motivated by my desire to compute the spectrum of matrix representations of SO(3). Initially, I tried to avoid computing ρd. My hope was to use properties of SO(3) and homogeneous polynomials to show that eigenvalues of a given ASO(3) completely determine the spectrum of the representation of A. Anyway, I couldn’t see a way to do it, and in the meantime I developed an algorithm I spent most of this post describing to compute ρd. I still believe there is a way to compute the spectrum of SO(3) without explicitly computing ρdUpdate: I sort of figured out how. Most of the details are worked out in the final section.

However, since there is a recursive formula for computing ρd, namely

ρd=Sd(Aρd1)Ed,

I can exploit it to make a statement about the spectrum of ρd. And it is this: The eigenvalues of Aρd1, which are easy to compute, are the eigenvalues of ρd. Though not up to multiplicity. Clearly. ρd and Aρd1 have different shapes.

The spectrum of Aρd1 is the product of all eigenvalues A with all the eigenvalues of ρd1. Say, Af=λf and ρd1g=ωg, then fg is an eigenvector of Aρd1 with eigenvalue λω. So then to show λω is a eigenvalue of ρd I can try transforming fg into an eigenvector of ρd. Well, via the projection operator there is one candidate, Sd(fg). First a technical lemma.

Lemma. A Suppose basis tensors β, β~P1Pd1 are equivalent, ββ~. Then

  1. Sd(Aρd1)(β~β)=0
  2. ρdSd=Sd(Aρd1)

Immediately one can see that 2 follows from 1; because EdSd sends β to the equivalence class representative, EdSdββ. To understand why 1 is true, it’s best to consider how A generally acts on degree d homogeneous monomials, given that ρd1 is known. Let fd be a degree d homogeneous monomial; it can be factored into degree 1 and d1 monomials like β and β~. A acts on the degree 1 factor and ρd1 acts on the degree d1 factor, results of theses actions are simply multiplied into a single expression, and finally contributions from like degree d monomial terms are collected. Multiplication is commutative and bilinear so one can take different factorizations of fd and order the operations so that the final results are the same. Sd(Aρd1) is precisely the same procedure with the actions of A and ρd1 combined by tensor product instead of multiplication. Recall from the section above that the tensor product is bilinear, just like multiplication, though not commutative. The projection Sd seen from a different perspective symmetrizes tensors, and has the effect of allowing the x, y, or z, in the P1 coordinate, to commute to its proper place. Therefore Sd(Aρd1) is well defined on equivalent βP1Pd1.

With Lemma A established I can now prove that the spectrums of ρd and Aρd1 are the same.

Lemma. B βP1Pd1 is an eigenvector of Aρd1 with eigenvalue λ if and only if there exists α an eigenvector of ρd with eigenvalue λ.

The implication follows directly from the matrix equality asserted in 2 of Lemma A.

ρdSdβ=Sd(Aρd1)β=λSdβ

In the converse direction suppose αPd is an eigenvector of ρd. Equation 2 of Lemma A and SdEd=Id, the matrix identity of size d, imply:

Combining shows that λ is an eigenvalue of Aρd1.

Sd(Aρd1)Edα=λSdEdα  (Aρd1)Edα=λEdα

And there it is Aρd1 and ρd have the same spectrum,

σ(ρd)={ed,e(d1),,1,,ed1,ed}.

Of course, the multiplicities will be different. But what are those multiplicities?

Dimension of the Eigenspaces: Eigenvalue Multiplicities.

From the proof of Lemma B one sees that Sd projects the λeigenspace of Aρd1 into the λeigenspace of ρd. Peering into the projection however is difficult because ρd1 too has a projection component, and so on. The entire recursion has series of nested projections which hampers the analysis. Instead I can describe ρd as a batch computation. And yes, there is a batch method. Which is, by the way, the definition of how A acts on some monomial pPd: Factor p into d P1 factors and compute the action of A on each factor. Abusing the previously introduced notation

ρd=Sd AAAd copies Ed.

This Sd projection is a true symmetrizing operator on i=1dP1; and just like before it projects the λeigenspace of i=1dA into the λeigenspace of ρd. An eigenvector of i=1dA is built from the eigenvectors of A tensored d times. Let v1, v0, and v+1 be eigenvectors of A for eigenvalues eiθ, 1, and e+iθ, respectively. Then a λeigenvector of i=1dA looks like

v+1v1v+1v1v0,

where each permutation of the elements v in the dtensor is a different λeigenvector. Furthermore, Sd sends each permutation to the same λeigenvector of ρd. What I want is to count the number of λeigenvectors of ρd, and I can build them from the λeigenvectors of i=1dA, but I don’t have to count the permutations; I just need to count the unique combinations that produce a λeigenvector of ρd. To make counting a bit easier I represent a eilθeigenvector of i=1dA as word in +, , and | characters.

Hopefully it is clear that plus is v+1, minus is v1, and pipe is v0. The negative l configuration would be the same drawing with the l plus signs replaced with l minus signs; so e±ilθeigenspaces have the same dimension. Count the number of possible combinations like so

0+ pairs1+ pairsdl2 + pairs

Therefore eigenvalues e±ilθ of ρd have multiplicity

dim(Eig(ρd(A),e±ilθ))=dl2+1 .

Application: Spectrum of SO(3) Representation in Harmonic Polynomials

The representations of SO(3) that I have described so far have been in Pd, homogeneous polynomials of degree d. These are, however reducible representations, because degree d harmonic polynomicals form a SO(3)invariant subspace of Pd. There is certainly more to be said about harmonic representations of SO(3), but my interest right now is to use to the the dimension of the eigenspaces formula to prove degree d harmonic SO(3) representations have eigenvalues

{ed,e(d1),,1,,ed1,ed},

and each has multiplicity one.

Pd+2 is the direct sum of degree d harmonic polynomials Hd+2 and r2Pd; so inside ρd+2 is a Pd representation, ρd. Then to compute the dimension of the e±ilθeigenspaces just subtract the dimension of the e±ilθeigenspace in ρd+2 from the dimension of the e±ilθeigenspace in ρd.

When 0ld

dim(Eig(ρd+2(A),e±ilθ))dim(Eig(ρd(A),e±ilθ))=d+2l2dl2=dl2+1dl2=1

and when d<ld+2

dim(Eig(ρd+2(A),eilθ))=d+2l20+1=1 .

Finally I have an answer/proof of Fact 1.2.1 from these Akshay Venkatesh lecture notes :)