000
001002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
038
039
040
041
042
043
044
045
046
047
048
049
050
051
052
053
054

Supplementary material for: Harmonic Exponential Families

This document provides additional background and implementation details for the paper “Harmonic ExponentialFamilies on Manifolds”.
Code will be released on the project webpage.
A simple test script for the harmonic density on the circle(the generalized von-Mises) is included. This density is notstudied experimentaly in the paper, but the script shows thesimplicity and effectiveness of the algorithm.

1. Numerical stability: the log - Fourier - exp
trick
When the distribution is highly peaked, the unnormalized
probabilities ϕ(g) = exp (η · T (g)) can become extremely
large, causing numerical instability or overflow. A similar
problem is regularly encountered in probilistic
computaP
tions where one wishes to compute ln exp xi , for somelog-probabilities xi . The standard solution is to subtract themaximum x∗ = max xi before exponentiating:XXx∗ + ln
exp (xi − x∗ ) = ln
exp (xi )(1)
iiThis is stable because all numbers exp (xi − x∗ ) are between zero and one.
Our situation is slightly different, because we want to compute the Fourier transform of exponentially large quantities.
The Fourier transform is a linear transformation but not aplain sum, so the outputs (the moments) can be negativeeven when the inputs are all positive. The logarithm of anegative number is a complex number, so one way to makethe computation stable is to use the log-Fourier-exp trickwith a complex logarithm function.
Alternatively, we may notice that all we need is the ratio oflarge quantities, so we need not compute the logarithm atall. The algotihm becomes:
1. Compute ln ϕ = F −1 η.
2. Compute m = maxi ln ϕ(gi ).
3. Compute ϕ̄ = exp (ln ϕ − m).
4. Compute M = F ϕ̄.
0
5. Compute Ep(g|η) [T (g)] = M/M00.

The factor e−m cancels in the division.

2. Fast and stable evaluation of sphericalharmonicsIn order to compute the empirical moments one must evaluate the sufficient statistics at the data. For a spherical harmonic density, these sufficient statistics are spherical harλmonics Yml (denoted Tm0in the paper) of potentially highdegree l. We found that the spherical harmonics implementation in SciPy has acceptable performance for evaluations of all spherical harmonics up to degree 50 – 80,
but was much too slow to be used for orders up to 200. We
also found that for large orders the SciPy function produces
NaN values.
Faster and more stable algorithms for the evaluation ofspherical harmonics have been developed recently, but they
are quite complicated, involving code generators or extended floating point numbers (Sloan, 2013; Fukushima,2011). We developed a very simple, stable, and fast algorithm for evaluation of all spherical harmonics up to order L. Our method uses the algorithm of Pinchon & Hoggan (2007) for computing real-valued Wigner D-functions
λλ(g), and the fact that spherical har(α, β, γ) = Tmn
Dmn
λmonics Ym are the n = 0 column of this matrix:
lYml (φ, θ) = Dm0(φ, θ, 0).

(2)(up to an arbitrary scaling constant)Pinchon and Hoggan decompose the D matrix as
Dl (α, β, γ) = X l (α)J l X l (β)J l X l (γ),(3)

where X are matrices with sinusoids on the diagonal andanti-diagonal, and J l is a pre-computed symmetric orthogonal block matrix with 4 non-zero blocks.
Computing the entire matrix
D(φ, θ, 0) = blockdiag(D0 (φ, θ, 0), . . . , DL (φ, θ, 0)),and then selecting the n = 0 column from each block is
not feasible because it involves the multiplication of O(L)blocks of dimension O(L), which makes the complexity
O(L4 ).
A better approach is to construct a binary coefficient vector
c with ones in the positions corresponding to m = 0 (one055
056
057
058
059
060
061
062
063
064
065
066
067
068
069
070
071
072
073
074
075
076
077
078
079
080
081
082
083
084
085
086
087
088
089
090
091
092
093
094
095
096
097
098
099
100
101
102
103
104
105
106
107
108
109

Harmonic Exponential Families

110
111112113114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164

4. Action of SO(3) on the spherical harmonics
expansionfor each value of l), and then computeD(φ, θ, 0)c = X l (α)J l X l (β)J l c(4)in the right-associative order.
The vector c for L = 2 would be:c = (1; 0, 1, 0; 0, 0, 1, 0, 0)(the semicolons are there for readability only). One way
to think of c is as selecting the n = 0 columns from theblocks of D when left multiplied by D. Another way to understand this is to think of c = Y (0, 0) as the evaluation of
the spherical harmonics at the north-pole of the sphere, andthen D(φ, θ, 0)c rotates this vector to become the spherical harmonics evaluated at some other point (φ, θ) on thesphere.
By exploiting right-associativity, each multiplication is amatrix-vector multiplication with cost O(L2 ) instead ofmatrix-matrix multiplication with cost O(L3 ). Furthermore, the products involving X(·) take only linear time because X is very sparse (Pinchon & Hoggan, 2007). The
total cost of computing spherical harmonics is thus O(L2 )per degree or O(L3 ) for all degrees up to L. Since thenumber of spherical harmonics grows quadratically withthe maximum degree L, the complexity of this algorithm
when measured in terms of the number of spherical harmonics is O(N 3/2 ) Stability follows from the orthogonality of X and J, as demonstrated by Pinchon & Hoggan(2007).

3. A basis for L2 (H)We claim in the paper that one can find a subset of matrix elements of irreducible unitary representations to span
the space of square-integrable functions on a homogeneousspace. This will be clarified in this section.
Let G λ be the finite-dimensional subspace of L2 (G)λspanned by the matrix elements Umn. Let Hλ be the subλset of G consisting of right invariant functions (i.e. forf ∈ Hλ , we have f (gk) = f (g), ∀g ∈ G, k ∈ K). First,
notice that Hλ is indeed a subspace of G λ : Let f 0 , f 00 ∈ Hλand set f (g) = αf 0 (g) + βf 00 (g) for scalars α, β. Thenf (gk) = αf 0 (gk) + βf 00 (gk) = αf 0 (g) + βf 00 (g) = f (g).
So f ∈ Hλ and hence Hλ is a vector space.
Since Hλ is a finite-dimensional subspace one can obtaina basis for it by taking linear combinations of matrix eleλments Umn. In fact, one may choose the matrix elements
for the group in such a way that a subset of these form abasis for Hλ .

We claim in the paper that the rotation group SO(3) actsirreducibly on the spherical harmonics coefficients. Thatis to say, if x is a function on the sphere and x̂ = Fx
are its Fourier coefficients, then the Fourier coefficients of
x0 = R(g)x(p) = x(g −1 p) are given by x̂0 = Fx0 =U (g)x̂, where U (g) is a block-diagonal matrix with irreducible representations U λ as blocks.
Recall that the spherical Fourier transform is the expansionof a function on the sphere in terms of spherical harmonics
λYmλ (p) = Um0(gp ), where U λ are the irreducible unitary
representation of SO(3), and gp ∈ SO(3) is an elementin the coset that represents p. In the following derivation,
we will write a finite-dimensional vector of spherical harλ∗monics evaluated at p as Y (p) = (Y00 (p), . . . , Yλ∗(p))T =
U:,0 (gp ). In full detail:[R(h)x](p) = x(h−1 p)= x̂ · Y (h−1 p)= x̂ · U·,0 (h−1 gp )= x̂ · U (h−1 )U·,0 (gp )= x̂ · U (h−1 )Y (p)T−1= U (h(5))x̂ · Y (p)= U (h)x̂ · Y (p)= U (h)Fx · Y (p)= F −1 U (h)Fx
The spherical harmonics are not always defined equal toλUm0, but can have different scale factors depending on λ.
The result above still holds for other scalings, though.

References
Fukushima, T.
Numerical computation of sphericalharmics by extending exponent of floating point numbers. J. Geodesy, pp. 271–285, 2011. doi: 10.1007/
s00190-011-0519-2.
Pinchon, D. and Hoggan, P. E. Rotation matrices for real
spherical harmonics: general rotations of atomic orbitals
in space-fixed axes. Journal of Physics A: Mathematical and Theoretical, 40(7):1597–1610, February 2007.
ISSN 1751-8113.
Sloan, P.-p. Efficient Spherical Harmonic Evaluation.
Journal of Computer Graphics Techniques, 2(2):84–90,
2013. URL http://jcgt.org.

165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219