Iterative Residual Fitting for Spherical Harmonic
Transform of Band-Limited Signals on the Sphere:
Generalization and Analysis
Usama Elahi∗ , Zubair Khalid† , Rodney A. Kennedy∗ and Jason D. McEwen‡
† Research

School of Engineering, The Australian National University, Canberra, ACT 2601, Australia
of Science and Engineering, Lahore University of Management Sciences, Lahore 54792, Pakistan
‡ Mullard Space Science Laboratory, University College London, Surrey RH5 6NT, UK
usama.elahi@anu.edu.au, zubair.khalid@lums.edu.pk, rodney.kennedy@anu.edu.au, jason.mcewen@ucl.ac.uk

arXiv:1709.02503v1 [] 8 Sep 2017

∗ School

Abstract—We present the generalized iterative residual fitting (IRF) for the computation of the spherical harmonic transform (SHT) of band-limited signals on the sphere. The proposed
method is based on the partitioning of the subspace of bandlimited signals into orthogonal subspaces. There exist sampling
schemes on the sphere which support accurate computation of
SHT. However, there are applications where samples (or measurements) are not taken over the predefined grid due to nature of the
signal and/or acquisition set-up. To support such applications, the
proposed IRF method enables accurate computation of SHTs of
signals with randomly distributed sufficient number of samples.
In order to improve the accuracy of the computation of the SHT,
we also present the so-called multi-pass IRF which adds multiple
iterative passes to the IRF. We analyse the multi-pass IRF for
different sampling schemes and for different size partitions.
Furthermore, we conduct numerical experiments to illustrate
that the multi-pass IRF allows sufficiently accurate computation
of SHTs.
Index Terms—Spherical harmonics, basis functions, spherical
harmonic transform, residual fitting, band-limited signals, 2sphere (unit sphere).

I. I NTRODUCTION
Signals are defined on the sphere in a variety of fields including geodesy [1], computer graphics [2], cosmology [3], astrophysics [4], medical imaging [5], acoustics [6] and wireless
communication [7] to name a few. Spherical harmonic (SH)
functions [8] are a natural choice of basis functions for
representing the signal on the sphere in all these applications.
Analysis on the sphere is done in both spatial (spherical)
and spectral (spherical harmonic) domains. The transformation
between the two domains is enabled by the well known
spherical harmonic transform (SHT) [8], [9]. For harmonic
analysis and signal representation (reconstruction), the ability
to accurately compute the SHT of a signal from its samples
taken over the sphere is of great importance.
Sampling schemes have been devised in the literature for the
accurate and efficient computation of SHTs [10], [11]. However, the samples may not be available, in practice (e.g., [5],
Usama Elahi, Zubair Khalid and Rodney A. Kennedy are supported by
Australian Research Council’s Discovery Projects funding scheme (project no.
DP150101011). Jason D. McEwen is partially supported by the Engineering
and Physical Sciences Research Council (grant number EP/M011852/1).

[12]), over the grid defined by these sampling schemes. To support the computation of SHTs in applications where samples or
data-sets are not available on the pre-defined grid, least squares
fitting (LSF) methods have been investigated for efficient
computation of the SHTs [12]–[18]. LSF methods formulate
a large linear system of basis functions and then attempt to
solve it efficiently. However, due to memory overflow, it is not
suitable for systems with large band-limits, L > 1024 [19]–
[21]. To solve this problem, an iterative residual fitting (IRF)
method has been proposed in [19] as an extension of LSF
and incorporates a divide and conquer technique for the
computation of SHTs. The basic idea of IRF is to divide
the subspace spanned by all spherical harmonics into smaller
partitions and then perform least squares on each partition
iteratively. Although IRF is fast, it creates a less accurate
reconstruction [19] as the size of the harmonic basis increases
for large band-limits. To improve the reconstruction accuracy,
a multi-pass IRF approach is used which includes multiple
passes for fitting. This is same as IRF but it involves multiple
IRF operations rather than one. A variant of this scheme is
presented in [19], where reconstruction for 3D surfaces is
carried out by taking large number of samples.
In this paper, we present an IRF method for the computation
of the SHT of a band-limited signal in a general setting that
partitions the subspace of band-limited signals into orthogonal
subspaces, where each orthogonal subspace can be spanned
by different numbers of basis functions. We also formulate
multi-pass IRF to improve the accuracy of computation of
the SHT. We analyze multipass IRF for different choices of
partitioning of the subspace and sampling schemes [10], [11],
[19], [22] and show that the computation of the SHT converges
in all cases. We also show that the convergence is fast for the
partition choice considered in this work.
The remainder of this paper is organized as follows. Section
2 provides the necessary mathematical background and notation required to understand the work. IRF and multi-pass IRF
methods are formulated in section 3. In section 4, we carry out
accuracy analysis of the proposed IRF method for different
partition choices and sampling schemes. Finally, concluding
remarks are presented in section 5.

II. M ATHEMATICAL BACKGROUND
A. Signals on the Sphere
A point v̂
=
v̂(θ, φ) on the unit sphere
S2 , {v̂ ∈ R3 : |v̂| = 1}, is parameterized by
[sin θ cos φ, cos θ cos φ, cos θ]T ∈ S2 ⊂ R3 , where (.)T
represents the transpose, θ ∈ [0, π] represents the co-latitude
and φ ∈ [0, 2π) denotes the longitude. The space of square
integrable complex functions of the form g(θ, φ), defined
on the unit sphere, form a complex separable Hilbert space,
denoted by L2 (S2 ), with inner product defined as by [8]
Z
hg, hi ,
g(θ, φ)h(θ, φ) sin θ dθ dφ, g, h ∈ L2 (S2 ), (1)
S2

where (·) represents the complex conjugate operation. The
functions with finite induced norm kgk , hg, gi1/2 are referred
to as signals on the sphere.
B. Spherical Harmonics
Spherical harmonic (SH) functions, denoted by Y`m (θ, φ)
for integer degree ` ≥ 0 and integer order |m| ≤ `, serve as
complete basis for L2 (S2 ) [8]. Due to the completeness of the
SH functions, any function g on the sphere can be expanded
as
`
∞ X
X
m
(g)m
(2)
g(θ, φ) =
` Y` (θ, φ),
`=0 m=−`

(g)m
`

where
are the SH coefficients of degree ` and order m
and form the spectral domain representation of the signal g,
given by the spherical harmonic transform (SHT) defined as
Z
m
(g)m
,
f,
Y
=
f (θ, φ)Y`m (θ, φ) sin θ dθ dφ. (3)
`
`
S2

The signal g is band-limited at degree L if (g)m
` = 0 for all
` ≥ L, |m| ≤ `. A set of band-limited signals forms an L2
dimensional subspace of L2 (S2 ), denoted by HL .
III. G ENERALIZED I TERATIVE R ESIDUAL F ITTING
Here we present the generalization of the IRF method [5],
[19] for the computation of the SHT of the band-limited signal
g ∈ HL from its samples.
A. Iterative Residual Fitting (IRF) – Formulation
The IRF method is based on the idea to partition the
subspace HL into smaller subspaces and carry out leastsquares estimation on these partitions iteratively. In this way, a
large linear problem is divided into manageable small subsets
of linear problems. The subspace HL has graphical representation of the form shown in Fig. 1, which also represents the
SH (spectral) domain formed by the SH coefficients of the
band-limited signal in HL . We partition HL into K orthogonal
k
subspaces HL
, k = 1, 2, · · · , K, each of dimension Nk . We
analyse different choices for partitioning later in the paper.
k
We index the SH functions that span the subspace HL
as
Ykj , j = 1, 2, · · · , Nk . We also define (g)kj = g, Ykj .

Fig. 1: Spherical harmonic domain representation of a bandlimited signal in HL .

Given M samples (measurements) of the band-limited signal g ∈ HL , we wish to compute SH coefficients. By defining
a vector

G , g(θ1 , φ1 ), . . . , g(θM , φM )]T ,
(4)
of M measurements (samples) of the signal g ∈ HL on
the sphere and the matrix Yk , with entries {Yk }p,q =
Ykq (θp , φq ), of size M × Nk containing SH functions that
k
evaluated at M sampling points, the
span the subspace HL
vector gk = [gk1 , gk2 , · · · , gkNk ]T of SH coefficients can be
iteratively computed (estimated) in the least-squares sense as
g̃k = (YkH Yk )−1 YkH rk ,

(5)

where (.)H represents the Hermetian of a matrix and
rk = G −

k−1
X

Yk0 g̃k0 ,

r0 = G

(6)

k0 =1

is the residual between the samples of the signal and the
signal obtained by using the coefficients g̃k0 for k 0 =
1, 2, . . . , k − 1 and the estimation of coefficients is carried out
iteratively for k = 1, 2, . . . , K. We note that the computational
complexity for (5) for each k would be of the order of
max(O(M Nk2 ), O(Nk3 )) = O(M Nk2 ). The computational
complexity to compute (6) is O(M L2 ). We later analyse the
estimation accuracy of the IRF method for different sampling
schemes on the sphere and different partitions of the subspace
HL of band-limited signals. For a special case of partitioning
k
the subspace HL into L subspaces HL
based on the degree
of spherical harmonics ` = 0, 1, . . . , L − 1, it has been shown
that the iterative residual fitting allows sufficiently accurate
estimation of SH coefficients [19].
The proposed IRF method enables accurate computation
of the SHT of signals with a sufficient number of randomly
distributed samples. The IRF algorithm finds significance use
in applications where samples on the sphere are not taken over
a predefined grid. For example, the samples are taken over
the cortical surface in medical imaging [5], where IRF allows
sufficient accurate parametric modeling of cortical surfaces.

B. Multi-Pass IRF and Residual Formulation
To improve the estimation accuracy, we employ the socalled multi-pass IRF [19] which is based on the use of
IRF method in an iterative manner. In multi-pass IRF, the
IRF algorithm is run for a number of iterations, denoted
by i = 1, 2, . . .. To clarify the concept, we incorporate the
iteration index i in the formulation in (5) and (6) as
g̃k (i) = (YkH Yk )−1 YkH rk (i),
rk (i) = G −

i−1 k−1
X
X

(7)

Yk0 g̃k0 (i0 ),

i0 =1 k0 =1

r0 (i) = rK (i − 1), r0 (1) = G.

(8)

After i-th iteration, g̃k can be computed for each k =
1, 2, . . . , K as
g̃k (i) =

i
X

g̃k (i0 ).

(9)

i0 =1

By defining
−1

Ak , (YkH Yk )

YkH ,

Ck , Yk Ak ,

the residual after the i-th iteration is given by
!i
K
Y
rK (i) =
(1 − Ck ) G.

(10)

(11)

k=1

In general, the residual in (11) depends on the distribution
of sampling points and nature of partitioning of HL . In the
next section, we show that the residual converges to zero for
a variety of sampling schemes and different partitions.
IV. A NALYSIS OF M ULTI -PASS IRF
A. Partition Choices
In order to understand the partitioning of HL , we refer to
the graphical representation of HL shown in Fig. 1, which
describes the position of spherical harmonic coefficients with
respect to degree ` ∈ (0, 1, . . . , L − 1) and order m ≤ |`|.
We give numbers to the spectral harmonic coefficients (basis
functions) shown in Fig. 1 from 1 to L2 in a way that we
start the domain from ` = 0, m = 0 and then traverse the
whole domain by m = −` to m = ` for increasing values of
`. In a similar way, we can also traverse the whole domain by
fixing m for all values of `. We analyse four different type of
partitions, whose sizes vary with the increasing or decreasing
values of degrees ` and orders m. The size of each partition
is denoted by Nk . In all the partitions, the generalized IRF is
run for all values of k and for a fixed value of i.
Partition Choice 1: We first consider the partitioning of
HL based on the spherical harmonic degree [19]. We take
k
K = L partitions HL
each for degree ` = k − 1 such that
k
the subspace HL is spanned by spherical harmonics of degree
k − 1. Consequently, the dimension of each subspace is Nk =
2k −1. As mentioned earlier, the IRF has been applied already
for this choice of partition [19]. We show through numerical
experiments that alternative choices for partitioning result in
faster convergence and more accurate computation of the SHT.

Partition Choice 2: For partition choice 2, we combine the
k-th partition choice 1 and K − k + 1-th partition choice 1,
partitions for even or odd band-limit L
to obtain L2 or L+1
2
k
respectively. For even L, each partition 2 HL
is of size Nk =
L
partitions
2L for k = 1, 2, . . . , 2 . For odd L, we have L+1
2
with Nk = 2L for k = 1, 2, . . . , L−1
and
one
partition
of size
2
N L+1 L.
2
Partition Choice 3: Here, we consider partitioning with
respect to each order |m| < L (see Fig. 1). Consequently, we
have 2L−1 partitions, one for each order |m| < L and spanned
by SH functions of order m.
Partition Choice 4: Partition choice 4 is obtained by
combining the partitions in partition choice 3. We obtain L
partitions by combining partition choice 3 for m and −(L−m)
for m = 1, 2, . . . L − 1. With such combining, we have L
partitions of HL each of size L.
B. Analysis
Here we analyse the accuracy of the computation of the
SHT, that is, the computation of SH coefficients, of the bandlimited signal sampled over different sampling schemes. For
the distribution of samples on sphere, we consider equiangular
sampling [10] and optimal-dimensionality sampling [11] in our
analysis as these schemes support the accurate computation
of the SHT for band-limited signals. Among the sampling
schemes on the sphere, which do not support the highly
accurate computation of the SHT, we consider the HEALPix
sampling scheme [22] and random samples with uniform
distribution with respect to the differential measure sin θdθdφ.
In order to analyse accuracy, we take a test signal g ∈ HL
by first generating the spherical harmonic coefficients (g)m
`
with real and imaginary part uniformly distributed in [−1, 1]
and using (2) to obtain the signal over the samples for each
sampling scheme. For a meaningful comparison, we take
approximately the same number of points for each sampling
scheme. We apply the proposed multi-pass IRF for each
choice of partition and each sampling scheme to compute the
estimate of SH coefficients (g̃)m
` and record the maximum
error between reconstructed and original SH coefficients given
by
max ,

max

`<L, |m|≤`

m
|(g)m
` − (g̃)` |,

(12)

which is plotted in logarithmic scale in Fig. 2 for band-limit
L = 15. Different partition choices and different sampling
schemes (see caption for number of samples for each sampling
scheme) against the number of iterations of the proposed
multi-pass IRF are plotted, where it can be observed that
1) the error converges to zero (10−16 , double precision) for
all partition choices and sampling schemes, and 2) the error
converges quickly for partition choice 4. We also validate
the formulation of the residual in (11) by computing after
each iteration of the multi-pass IRF. To illustrate the effect
of the number of samples on the accuracy of the proposed
multi-pass IRF, we have taken 2L2 , 4L2 and 6L2 samples of
optimal dimensionality sampling [11] and plot the error max

5

5
Partition1
Partition2
Partition3
Partition4

5
Partition1
Partition2
Partition3
Partition4

0
-5

-5

-10

-10

-10

-15

-15

-15

-20

-20
0

20

40

60

80

100

-20
0

20

No of Iterations, i

40

60

80

100

0

0

ǫ max

-5

-5

-5

-10

-10

-15

-15

-15

-20

-20
40

60

80

100

No of Iterations, i

100

Partition1
Partition2
Partition3
Partition4

0

-10

20

80

5
Partition1
Partition2
Partition3
Partition4

ǫ max

0

60

(c)

5
Partition1
Partition2
Partition3
Partition4

40

No of Iterations, i

(b)

5

0

20

No of Iterations, i

(a)

ǫ max

Partition1
Partition2
Partition3
Partition4

0

ǫ max

-5

ǫ max

ǫ max

0

-20
0

20

40

60

No of Iterations, i

(d)

(e)

80

100

0

20

40

60

80

100

No of Iterations, i

(f)

Fig. 2: Maximum reconstruction error max , given in (12), between the original and reconstructed SH coefficients of a bandlimited signal with L = 15. Reconstructed SH coefficients are obtained using the proposed multi-pass IRF, where the samples
of the signal are taken as (a) 991 samples of the Equiangular sampling scheme, (b) 972 samples of the HEALpix sampling
scheme, (c) 900 random samples (d) 450 (e) 900 and (f) 1350 samples of the optimal dimensionality sampling scheme.

in Fig. 2(d)-(f), where it is evident that the error converges
quickly for a greater number of samples. The convergence of
the error is in agreement with the formulation of the residual in
(11), however, convergence changes with the sampling scheme
and nature of the partition of the subspace of band-limited
signals. This requires further study and is the subject of future
work.
V. C ONCLUSIONS
We have presented the generalized iterative residual fitting (IRF) method for the computation of the spherical harmonic transform (SHT) of band-limited signals on the sphere.
Proposed IRF is based on partitioning the subspace of bandlimited signals into orthogonal spaces. In order to improve the
accuracy of the transform, we have also presented a multi-pass
IRF scheme and analysed it for different sampling schemes
and for four different size partitions. We have performed
numerical experiments to show that accurate computation of
the SHT is achieved by multi-pass IRF. For different partitions
and different sampling distributions, we have analysed the
residual (error) and demonstrated the convergence of the
residual to zero. Furthermore, it has been demonstrated that the
rate of convergence of error depends on the sampling scheme
and choice of partition. Rigorous analysis relating the nature
of partitions and convergence of the proposed method and
application of proposed method in medical imaging, computer
graphics and beyond are subjects of future work.
R EFERENCES
[1] A. Amirbekyan, V. Michel, and F. J. Simons, “Parametrizing surface
wave tomographic models with harmonic spherical splines,” Geophys.
J. Int.l, vol. 174, no. 2, pp. 617–628, Jan. 2008.

[2] R. Ng, R. Ramamoorthi, and P. Hanrahan, “Triple product wavelet
integrals for all-frequency relighting,” ACM Trans. Graph., vol. 23, no. 3,
pp. 477–487, Aug. 2004.
[3] Y. Fantaye, C. Baccigalupi, S. Leach, and A. Yadav, “Cmb lensing
reconstruction in the presence of diffuse polarized foregrounds,” J.
Cosmol. Astropart. Phys., vol. 2012, no. 12, p. 017, July 2012.
[4] N. Jarosik, C. L. Bennett, J. Dunkley, B. Gold, M. R. Greason,
M. Halpern, R. S. Hill, G. Hinshaw, A. Kogut, E. Komatsu, D. Larson,
M. Limon, S. S. Meyer, M. R. Nolta, N. Odegard, L. Page, K. M. Smith,
D. N. Spergel, G. S. Tucker, J. L. Weiland, E. Wollack, and E. L.
Wright, “Seven-year wilkinson microwave anisotropy probe (wmap)
observations: Sky maps, systematic errors, and basic results,” ApJS, vol.
192, no. 2, p. 14, Jan. 2011.
[5] M. Chung, K. Dalton, L. Shen, A. Evans, and R. Davidson, “Weighted
fourier series representation and its application to quantifying the amount
of gray matter,” IEEE Trans. Med. Imag., vol. 26, no. 4, pp. 566–581,
Apr. 2007.
[6] W. Zhang, M. Zhang, R. A. Kennedy, and T. D. Abhayapala, “On
high resolution head-related transfer function measurements: An efficient sampling scheme,” IEEE/ACM Trans. Audio, Speech, Language
Process., vol. 20, no. 2, pp. 575–584, Dec. 2012.
[7] Y. F. Alem, Z. Khalid, and R. A. Kennedy, “3D spatial fading correlation
for uniform angle of arrival distribution,” IEEE Commun. Lett., vol. 19,
pp. 1073–1076, June 2015.
[8] R. A. Kennedy and P. Sadeghi, Hilbert Space Methods in Signal
Processing. Cambridge, UK: Cambridge University Press, March 2013.
[9] J. J. Sakurai, Modern Quantum Mechanics, 2nd ed. Reading, MA:
Addison Wesley Publishing Company, Inc., 1994.
[10] J. D. McEwen and Y. Wiaux, “A novel sampling theorem on the sphere,”
IEEE Trans. Signal Process., vol. 59, no. 12, pp. 5876–5887, Dec 2011.
[11] Z. Khalid, R. A. Kennedy, and J. D. McEwen, “An optimaldimensionality sampling scheme on the sphere with fast spherical
harmonic transforms,” IEEE Trans. Signal Process., vol. 62, no. 17,
pp. 4597–4610, Sep. 2014.
[12] N. Sneeuw, “Global spherical harmonic analysis by least-squares and
numerical quadrature methods in historical perspective,” Geophys. J.
Int., vol. 118, no. 3, pp. 707–716, Sep. 1994.
[13] D. M. Healy, Jr., D. Rockmore, P. J. Kostelec, and S. S. B. Moore,
“FFTs for the 2-sphere - improvements and variations,” J. Fourier Anal.
and Appl., vol. 9, pp. 341–385, 2003.

[14] J. Blais and M. Soofi, “Spherical harmonic transforms using quadratures
and least squares,” in Computational Science ICCS 2006, ser. Lecture
Notes in Computer Science. Springer Berlin Heidelberg, 2006, vol.
3993, pp. 48–55.
[15] K. Ivanov and P. Petrushev, “Irregular sampling of band-limited functions on the sphere,” Appl Comput Harmon Anal., vol. 37, no. 3, pp.
545 – 562, Nov. 2014.
[16] J. Keiner, S. Kunis, and D. Potts, “Efficient reconstruction of functions
on the sphere from scattered data,” J. Fourier Anal. Appl., vol. 13, no. 4,
pp. 435–458, May 2007.
[17] S. Kunis and D. Potts, “Fast spherical fourier algorithms,” J. Comput.
Appl. Math., vol. 161, no. 1, pp. 75 – 98, Dec. 2003.
[18] ——, “Stability results for scattered data interpolation by trigonometric
polynomials,” SIAM J. Sci. Comput., vol. 29, no. 4, pp. 1403–1419, Feb.
2007.
[19] L. Shen and M. K. Chung, “Large-scale modeling of parametric surfaces
using spherical harmonics,” in 3D Data Processing, Visualization, and
Transmission, Third International Symposium on, June 2006, pp. 294–
301.
[20] C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations.
Society for Industrial and Applied Mathematics, 1995.
[21] R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra,
V. Eijkhout, R. Pozo, C. Romine, and H. van der Vorst, Templates for
the Solution of Linear Systems: Building Blocks for Iterative Methods.
Society for Industrial and Applied Mathematics, 1994.
[22] K. M. Grski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen,
M. Reinecke, and M. Bartelmann, “Healpix: A framework for highresolution discretization and fast analysis of data distributed on the
sphere,” Astrophys. J., vol. 622, no. 2, p. 759, Apr. 2005.

