Near Optimal Reconstruction of Spherical Harmonic Expansions
Authors: Amir Zandieh
Insu Han
Haim Avron
Venue: NeurIPS 2023 poster (2023)
Keywords: spherical harmonic expansion
Gegenbauer polynomials
interpolation
leverage score sampling
TL;DR: We present an algorithm for efficiently recovering spherical harmonic expansions of functions on the $d$-dimensional unit sphere, requiring a near-optimal number of function evaluations.

Abstract: We propose an algorithm for robust recovery of the spherical harmonic expansion of functions defined on the $d$-dimensional unit sphere $\mathbb{S}^{d-1}$ using a near-optimal number of function evaluations. We show that for any $f\in L^2(\mathbb{S}^{d-1})$, the number of evaluations of $f$ needed to recover its degree-$q$ spherical harmonic expansion equals the dimension of the space of spherical harmonics of degree at most $q$, up to a logarithmic factor. Moreover, we develop a simple yet efficient kernel regression-based algorithm to recover degree-$q$ expansion of $f$ by only evaluating the function on uniformly sampled points on $\mathbb{S}^{d-1}$. Our algorithm is built upon the connections between spherical harmonics and Gegenbauer polynomials. Unlike the prior results on fast spherical harmonic transform, our proposed algorithm works efficiently using a nearly optimal number of samples in any dimension $d$. Furthermore, we illustrate the empirical performance of our algorithm on numerical examples.

Near Optimal Reconstruction of Spherical Harmonic Expansions
Amir Zandieh Independent Researcher amir@zed512. .ma...mal....... com
Insu Han Yale University insu. haauu.aa@@all...mm edu
Haim Avron
Tel Aviv University
haimavStauex. tau. tau.ac.ii ac i1
Abstract
We propose an algorithm for robust recovery of the spherical harmonic expansion of functions defined on the d-dimensional unit sphere gd using near-optimal number of function evaluations. We show that for any f E (= the number of evaluations of f needed to recover its degree-q spherical harmonic expansion equals the dimension of the space of spherical harmonics of degree at most g. up to logarithmic factor. Moreover, we develop simple yet efficient kernel regression- based algorithm to recover degree-g expansion of f by only evaluating the function uniformly sampled points on gd-1 Our algorithm is built upon the connections between spherical harmonics and Gegenbauer polynomials. Unlike the prior results fast spherical harmonic transform, our proposed algorithm works efficiently using a nearly optimal number of samples in any dimension d. Furthermore, we illustrate the empirical performance of algorithm numerical examples.
1 Introduction
We consider the fundamental problem of recovering a function from a finite number of (noisy) observations. To provide accurate and reliable predictions at unobserved points need to avoid overfitting which is typically achieved through restricting our estimator or interpolant to a family of smooth or structured functions. In this paper, focus interpolating gguare-integrabe functions on the d-dimensional unit sphere, with low-degree spherical harmonies, critical task in scenarios where rotational invariance fundamental property. Spherical harmonies are essential in various theoretical and practical applications, including the representation of electromagnetic fields [Wei95], gravitational potential [Wer97], cosmic microwave background radiation [KKS97] and medical imaging [CLLI5], as well as modelling of 3D shapes in computer graphics [KFR03] and computer vision [2CL+ Further notable real-life applications include molecular/atom systems, where understanding the underlying functions within a spherical context can significantly enhance predictive modeling and simulation accuracy [EEHM17, 211 ZDK 22].
We begin by observing that any function f in L (S(()))), i.e., the family of gquarr-integrable functions on the sphere gd-1, can be uniquely decomposed into orthogonal spherical harmonic components. Specifically, if we denote the space of spherical harmonics of degree in dimension d by Hs
any f ( has unique orthogonal expansion f Σe= fe with ft He (Sd-1) (Lemma 2). With this observation, we aim to find the best spherical harmonic approximation of degree S q to f using minimal number of samples (essentially treating higher order terms in f's expansion as noise).
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Problem 1 (Informal Version of Problem 2). For an unknown function f f ()" and an integer q 2 efficiently (both terms of number of samples from f and computations) learn the first q + 1 spherical harmonic components {fe E H:(Sa-1) 1e- of f which minimize
12
Σ Se(w) - f(w) dw.
8d-1
The angular power spectrum of f commonly obeys power law decay of the form fl S O(E-)), for some $ > 0, depending on the order of differentiability of f. In fact, for any infinitely differentiable f. [.fellgaa. decays asymptotically faster than any rational function of E. Furthermore, for any real analytic f on the sphere, l[fellge decays exponentially. Thus, the first q 1 spherical harmonic components of f typically well approximate f for even modest and answering Problem meaningful for wide range of differentiable functions.
1.1 Our Main Results
We reformulate Problem I as a least-squares regression and then solve it using randomized numerical linear algebra techniques. We first consider orthonormal projection operator that maps functions in L I'((5-1) onto the space of bounded-degres spherical harmonics D'=o He He (Sd-1). Specifically, x:le) is operator that maps any f with expansion Er= fe where ft E He onto f's first q + 1 expansion components, i.e., x()) Ei=o fe, then Problem can be formulated as
|223 g- - -|
min 9EL-(80-11)
However, solving this regression problem with "continuous" cost function is challenging. To avoid such continuous optimizations, we adopt the approach of [AKM* 19] which discretizes the afore- mentioned regression problem according to the leverage scores of the operator X:(") It turns
that if we can draw random samples with probabilities proportional to the leverage scores of K((2) then we can recover the degree-q spherical harmonic expansion of f. i.e. Ei with finite number of observations. Particularly, by exploiting the connections between spherical harmonics and Zonal (Gegenbauer) Harmonics and the fact that zonal harmonics are the reproducing kernels of (e((st))) (Lemma 3), we prove that the leverage scores of K( are constant everywhere on the sphere gd-1, Thus, solving discrete regression problem with uniformly sampled observations yields a near-optimal solution to Problem 1, Informal statements of our results are as follows.
Theorem (Informal Version of Theorem 5). Let Bagd be the dimension of spherical harmonics of degree at most q, i.e. Bqdd = dim (Dr I)) There exists an algorithm that finds (1 + e) -appproxiimaooon to the optimal solution of Problem 1, given (Bg.a a log Pgg. + (gg)) observations of function f at uniformly sampled points on gd-1, with (s² + sh)² runtime.
We also prove that our bound on the number of required samples is optimal up to a logarithmic factor. Theorem 2 (Informal Version of Theorem 6). Any (randomized) algorithm that takes 8 < Bg. samples on any input fails with probability greater than 9/ /1, where Bg.d = dim (+i He 1))
1.2 Related Work
Efficient reconstruction of functions as per Problem 1 has been extensively studied in various fields. Many prior papers considered reconstructing 1-dimensional functions from finite number of samples on a finite interval under smoothness assumption about the underlying function. Notably, the influential line of work of [SP61, LP61, LP22, XRYO1] focused on reconstructing Fourier-bandlimited functions and [CKPS16, KVZ19, EMM20, BKM+ considered interpolating Fourier-sparse signals. Recently, [AKM + 19] unified these reconstruction methods in dimension d = 1 and gave a universal sampling framework for reconstructing nearly all classes of functions with fourier-basee smoothness constraints. One can view 1-dimensional functions on a finite interval as functions on the unit circle
ww 2.3727 is the exponent of the fast matrix multiplication algorithm [Wil12]

---

8', Thus, Problem 1 is indeed generalization of the aforementioned prior work to high dimensions under the assumption that the generalized Fourier series (Lemma 2) of the underlying function only contains bounded-degree spherical harmonics. This degree constraint on spherical harmonic expansions can be viewed as the d-dimensional analog of Fourier- manrier-aandlimitedeesss on circle S!.
Computing the spherical harmonic expansion in dimension d = 3 has received considerable attention in physics and applied mathematics communities. The algorithms for this special case of Problem I are known in the literature as "fast spherical harmonic transform" [SS00, STO2J. Most notably, [RT06] proposed an algorithm for computing spherical harmonic expansion of degree I q to precision E using O($g3) samples and 433 log $q.3 (og(1/=) time. These fast algorithms were developed based on fast Fourier and associated Legendre transforms and make use of colliconditioned) orthogonal basis for Hy (8d-1) ) which happened to be the associated Legendre polynomials when d = 3. However, it is general very difficult to compute an orthogonal basis for spherical harmonics [MNY06], so unlike our Theorem 1, it is inefficient to extend these prior results to higher d.
From the techniques point of view, a related work [GMMM21], which employs harmonic analysis over 3d-1 to analyze the generalization of two-layered neural tangent kernels. They show that an unknown function defined gd-1 can be efficiently recovered using kernel regression w.r.t neural tangent kernel on uniform random samples from the function. However, the number of samples that [GMMM21] requires for recovering bounded degree spherical harmonics, especially when the degrees are high, is sub-optimal and is strictly worse than our result. Additionally, [GMMM21] I] does not guarantee recovery with relative error, while our Theorem provides relative error guarantees.
Furthermore, recent applications of Gegenbauer polynomials, along with other orthonormal poly- nomials like Hermite polynomials, have been found in designing efficient random features for approximating various kernel functions. These applications extend to dot-product kernels [HZA22], Neural Tangent Kernels [ZHA*
+21, 21, H2L+ and Gaussian kernels [WZ22, AKK +20.
2 Mathematical Preliminaries
We denote by gd-1 the unit sphere in d dimension. We use |84-1|=
r(d/2) to denote the surface area of sphere gd-1 and (S(-)1) denote the uniform probability distribution on S6-1, We denote by L (gd-1) the set of all square-integrable real-valued functions sphere gd-1, Furthermore, for any f.9 2 (sd-1) we use the following definition of the inner product on the unit sphere?.
2/2
(fig)ss 1 = Jea-1 f(w)g(w)dv = |s--1|
J
E
[f(w)g(w]]
The function space 2 (Sa-1) is complete with respect to the norm induced by the above inner product, i.e. I|f == (f,f)s so L (Su-1) is Hilbert space.
We often use the term quasi-matrix which is informally defined as a "matrix" in which one dimension is finite while the other is infinite. A quasi-matrix matrix can be tall (or wide) meaning that there is a finite number of columns (or rows) where each one is functional operator. For more formal definition, see [SA22].
Spherical Harmonics are the solutions of Laplace's 's equation in spherical domains and can be thought of as functions defined on Sd-1 employed in solving partial differential equations. Formally,
Definition (Spherical Harmonics). For integers 2 0 and d 2 1, let Pe(d) be the space of degree-E homogeneous polynomials with d variables and real coefficients. Let Held) denote the space of degree-l harmonic polynomials in dimension d, i.e., homogeneous polynomial solutions of Laplace 's equation:
H=(d)==(Pe = (dd): AP
where A 8 du?
82 is the Laplace operator on RRd, Finally, let He (gd-1) be the space of des
(real) Spherical Harmonics of order É in dimension d, i.e. restrictions of harmonic polynomials in
Formally, 2 (sd-1) is a space of equivalence classes of functions that differ at set of points with measure 0. For notational simplicity, here and throughout we use I to denote the specific representative of the equivalence class f L (sd-1). In this way, consider the point-wis wise value f(w) for every gd-1
3

---

(e(d) to the sphere S The dimension of this space, aed = dim (Hl (8 (s+-1)), is
d ( ( 1 - 1000-(12
1 8-2 (
00,d = 1, add aldd=, d, ai1d =
for l 2 2.
2.1 Gegenbauer Polynomials
The Gegenbauer (a.ka. ultraspherical) polynomial of degree 8 2 0 in dimension d 2 is given by
L2/21
(1(1) = oj
Σ
2-22__ (1 - -
(3)
==
where co = and Cy+1 ((-2))((-2)-1) 2(3+1)(-1+22)
Cj for 0, ,..,,..,. 1,/2, 1. These polynomials are
orthogonal on the interval [-1, 1] with respect to the measure (1 - -22 - i.e.
[ PP(().PP
gd-1
Cgg, gd (e=e))-
Zonal Harmonies. The Gegenbauer polynomials naturally provide positive definite dot-product kernels S6-1 known as Zonal Harmonies, which are closely related to the spherical harmonics. The following reproducing property of zonal harmonics plays a crucial role our analysis.
Lemma (Reproducing Property of Zonal Harmonics). Let P.(.) be the Gengenbauer polynomial of degree E in dimension d. For any x, y 1
(1, y)) = orea (B [P] ((x, w)) Pÿ (ty, w))]
Furthermore, for any l + R:
E ( [Pá ((x,w))) Py (1y, (y)
The proof of this and all subsequent results can be found in the appendix. The following useful fact, known as the addition theorem, connects Gegenbauer polynomials and spherical harmonics.
Theorem (Addition Theorem). For every integer 2 0, if { ve..u...... .. } is an orthonormal
basis for He (Sd-1), then for any o, w E Sd-1 we have
agdd Sd-1
Pd w)) - (
3 Reconstructing 1.2 (Sd 1) Functions via Spherical Harmonics
In this section show how to approximate any function f E T,2 (5d-1) by spherical harmonics using the optimal number of samples. We begin with the fact that spherical harmonics form a complete set of orthonormal functions and thus form orthonormal basis for the Hilbert space of square-integrable functions on sphere gd-1 This is analogous to periodic functions, viewed functions defined the circle S1, which functions (sines and cosines) via the Fourier series.
be expressed as linear combination of circular
Lemma 2 (Direct Sum Decomposition of L2(gd-1)). The family of spaces He (Sd-1) yields Hilbert space direct sum decomposition 1 (gd-1)- Di He (sd-1): the summands are closed and pairwise orthogonal, and every f E I (gá-1) is the sum of converging series (in the sense of mean-square convergence with the norm defined in Eq. (2)),
Σ St.
E=
where fe He (1a-1) are uniquely determined functions. Furthermore, given any orthonormal busis { uh ub............... Van of (gd-1) we have St Σg (2 (11

---

The series expansion in Lemma 2 is the analog of the Fourier expansion of periodic func- tions, and is known as "generalized Fourier series" [Pen30] with respect the Hilbert basis {uf j E [ov.a]. 2 0}. We remark that it is in general intractable to compute an orthogonal basis for the space of spherical harmonics [MNY06], which renders the generalized Fourier series expansion in Lemma 2 primarily existential. While finding the generalized Fourier expansion of function f E L (s=-1) is computationally intractable, our goal to answer the next fundamental question, which is about finding the projection of a function f onto the space of spherical harmonics, i.e., the ft's Lemma 2. Concretely, we seek to solve the following problem.
Problem 2. For an integer q 2 0 and a given function f E L2 (sd-1) whose decomposition over the Hilbert sum Dz-o He (Sd-1) is Ee=o St as per Lemma 2, let us define the low-degree expansion of this function as f(a) :- Ei ft. How efficiently can we learn f(a) E +ii o He (sd-1)? More precisely, we want to find set {wi, w2..... Ws) } C gd-1 with minimal cardinality along with an efficient algorithm that given samples (f(w:)) can interpolate f()) with function f(a) E Doo He (s0-1) such that:
2
f(
S se. ga
gd
For case of notation, we denote the Hilbert space of spherical harmonics of degree at most by H(() (s*=1)= == +2=o Me (S To answer Problem 2 we exploit the close connection between the spherical harmonics and Gengenbauer polynomials, and particular the fact that zonal harmonies are the reproducing kernels of the Hilbert spaces He (S=-1).
Lemma (A Reproducing Kernel for He (Su-1)), For every E L2 (8"-1), iffff=Eee fe the unique decomposition of f over Di Hi (Sa-1) as per Lemma 2, then fe is given by
fe(a)= = ag,
E (B
(w ((a,w))) w))] for E Sd-1
Now we define a kernel operator, based on the low-degree Gegenbauer polynomials, which projects functions onto their low-degree spherical harmonic expansion.
Definition (Projection Operator onto H(4) (gd-1)). For any integers q 2 0 and 2 define the kernel operator x:(/) (5u-1) + I (Su-1) as follows: for f T I (Sd-1) and a E S
1 1) (0) = Σ
atd |dd Pf (la, 9a- = =Ears aad E 2=
[f(w) ((a, w))] (5)
aed This is an integral operator with kernel function kq,alo, w) = E'= g--
Pd ((a,w)).
Note that the operator K(()) is self-adjoint and positive semi-definite. Moreover, using the reproducing property of this kernel we can establish that k:(2) is projection operator.
Claim 1. The operator x(g) defined in Definition 2 satisfies the property (x:( k(()
Furthermore, by the addition theorem (Theorem 3), xc(e) is trace-class (i.c., the trace is finite and independent of the choice of basis) because:
trace (Ke(2) = JSa
- knalw, w) dw-
atdd Sd-1
Pd ((w, w)) dw
Iga-1
E=
ard = (d (d+g-1)
d q - 2) - ?) -1.
(6)
By combining Theorem 3 and Lemma 2, and using the definition of the projection operator kc()), follows that for any function I!2 (sd-1) with Hilbert sum decomposition - ft, the low- degree component f(1) El=o ft H(9) (sd-1) can be computed f(a) kC() Equivalently,
5

---

in order to learn f(A), it suffices to solve the following least-squaree regression problem,
k((0) 9
1 -
min g
1)
If g" is an optimal solution to the above regression problem then f(() = x-(1) g" In the next claim we show that solving the least squares problem Eq. (7), even coarse approximation, sufficient to solve our interpolation problem (i.e., Problem 2):
Claim 2. For any f L2 (S--1) any integer 2 O, and any C2 221, 12gifeELL ( satisfies,
|K]9 g- I1p.. SC.
1C 9 1
min (E
and if we let f(z) = K/L) f. where xg)) is defined per Definition 2, then the following holds
|| g - p(w) | |3=_ ga-1 - (C)
1) 11)1119900
gd-1
Claim 2 shows that solving the regression problem in Eq. (7) approximately provides a solution to our spherical harmonics interpolation problem (Problem 2). But how can we solve this squt-squares problem efficiently? Not only does the problem involve a possibly infinite dimensional parameter vector g. but the objective function also involves the continuous domain on the surface of g
3.1 Randomized Discretization via Leverage Function Sampling
We solve the continuous regression in Eq. (7) by randomly discretizing the sphere Sd thereby reducing our problem to regression on a finite set of points W1, .2..... W/ gd-1 In particular, we propose to sample points on gd-1 with probability proportional to the so-called leverage function, specific distribution that has been widely applied in randomized algorithms for linear algebra problems on discrete matrices [LMP13]. We start with the definition of the leverage function for compact operators such as K Definition 3 (Leverage Function). For integers q 2 0 and d > O, we define the leverage function of the operator x((() (see Definition 2) for every Sd-1 as follows,
TT(w)= :
max
|c//
|
||<]0) g)))))()) (w) (w))
(8)
9E 9
ga-1
Intuitively, Ty (w)) is an upper bound how much a function that is spanned by the eigenfunctions of the operator xc(9) can "blow up" at The larger the leverage function (w) implies the higher the probability we will be required to sample w. This ensures that our sample points well reflect any possibly significant components, or "spikes", of the function. Ultimately, the integral fga-1 Tyy Ty(w) dw determines how many samples we require to solve the regression problem Eq. (7) to a given accuracy.
It is known fact that the leverage function integrates to the rank of the operator k:(1) (which is equal the dimensionality of the Hilbert space (Sd-))))))). This ultimately allows us achieve OEE (g,a) sample complexity bound for solving Problem To compute the leverage function, we make use of the following useful alternative characterization of the leverage function.
Lemma 4 (Min Characterization of the Leverage Function). For any w E gd-1 let (w) be the and
leverage function (Definition 3) and define @u E L2d by @w( = Ei=o Te Pd ((a, w)).
We have the following minimization characterization of the leverage function:
Ty (w)
min LP
S.S. xc()) g 2w
Using the min and max characterizations of the leverage function we can find upper and lower bounds on this function. Surprisingly, in this case the upper and lower bounds match, so actually have exact value for the leverage function.
Lemma (Leverage Function Constant), The leverage function given Definition is equal to (T(w) -n B--
for every w E $d-1

---

Algorithm Efficient Spherical Harmonic Expansion
1: Input: accuracy parameter 5 > O, integer y 2 0
2: Set (Bq.d log Sq,d fq.a/e) for sufficiently large fixed constant
3: Sample i.i.d. random points W1, W2, Ws from (((-11)
4: Compute K E Rxx with K E
aea a-18--11 5: Compute f E R with f, = v f(w;) for E
Pa ((w,, wz)) for i,j E [s]
6: Solve the regression by computing 2 = f
7: Return: ((4)( with y(o) = Ei=o Vo-1a-11 Σj=ı 23 Pa ((wg, o)) for s
and
We prove this lemma in Appendix C. The integral of the leverage function, which determines the total samples needed to solve our least-squares regression, is therefore equal to the dimensionality of the Hilbert space H(9) ( (Sd-1).
Corollary 1. The leverage function defined in Definition 3 integrates to the dimensionality of the Hilbert space H(9) (Sd-1), which we denote by Bq di i.e.,
Joc Ty(1)) dw dim (160) = Σ
-
Or
Sq.d-
t=
We now show that the leverage function be used to randomly sample the points the unit sphere to discretize the regression problem in Eq. (7) and solve it approximately.
Theorem (Approximate Regression via Leverage Function Sampling). For any ε > let (Baa log Sy + - f.for (rfforuufff sufficiently large fixed constant c, and let 11; 22. be i.i.d. uniform samples on 8d-1 Define the quasi-matrix P IR - (2d()d-1) as follows, for every Rd:
[P)](o)= Σ
O'g,d
Uy Σ
Pd ((2,,))) for a ES
Vs Sd-
$=1
Also let f R be a vector with J == + [(x;))or for 1, .,,,, and let P be the adjoint of P. If ÿ is optimal solution to the least-squares problem ğ E arg mingget2(3e-1) P f f11 then with probability at least 10- the following holds,
|K g 1|
(1++))
min GEL2(36-1)
ga-1
We prove Theorem Appendix This theorem shows that the function ğ obtained from solving the discretized regression problem provides an approximate solution to Eq. (7).
3.2 Efficient Solution for the Discretized Least-Squares Problem
In this section, we demonstrate how apply Theorem 4 algorithmically approximately solve the regression problem of Eq. (7). To achieve this, we leverage the kernel trick, following a similar approach previous works such as [AKM 19, CP19], which allows us efficiently address the randomly discretized least squares problem as detailed in Algorithm 1. The associated guarantee for this approach provided in Theorem
Theorem 5 (Efficient Spherical Harmonic Interpolation). Algorithm 1 returns function ((d)() such that, with probability at least 1 - 10-4;
- p(o)
2 ga-1 S <e. - ga-1" where ;(a):= :00==Kk0 k((g) f.
Suppose we compute the Gegenbauer polynomial Pf(t) at every point E [-1, in constant time. Then Algorithm queries the function al O (B log 3r.d + Bad points the sphere gd-1 and runs in O(s² d + 5²)) time. This algorithm evaluates y(o) in O(d s) time for any
gd-1
7

---

We provide the proof of this theorem in Appendix D.
Remark (Noise Robustness). In Theorem 5, our method's robustmess is demonstrated under a noise model where the function f is not strictly low-degree spherical harmonic and may include high-degree components. In this scenario, the higher-degree components are considered as noise added to the inpul function.
However, our algorithm robust against alternative noise models, particularly additive i.i.. Normal noise that corrupts the function evaluations f in Algorithm I with iid Normal noise. More precisely, suppose that we observe samples from the function f(g) contaminated by Gaussian noise, i.e., f(;,) = f(a) (wz) + ny in Algorithm for i..dd n1. R2........ B2 (0, 1). The expected value of perturbation's 's norm in the output y of our algorithu caused by this noise is:
|| f())
(1/s) trace (K'KK').
g4-1
O11ss) trace (K'KKt). we By Markov's inequality, with 0.99 probability |ly -f let P be the operator defined in Theorem 4, then we can see that K P'P. Using the properties of the trace, one can see that trace ( K'KK") = 1/>1 + 1/A₂ + ... where X; ('s are the singular values of the operator PP By matrix Chernoff inequalities one can show that all singular values of the operator PP closely approximate the singular values of the projection operator K/()) up to constant factor. So, we have trace KK) O(rank(K)))) 0(3.a). Thus, because 8 2 (1BBe) and using union bound, with 0.98 probability:
2
p(
E.
ga-1
4 Lower Bound on The Number of Required Observations
We conclude by showing that the dimensionality of the Hilbert space (Sa-1) tightly characterizes the sample complexity of Problem 2. Thus, our Theorem is optimal up to a logarithmic factor. Intuitively, there are B.. degrees of freedom for specifying a spherical harmonic. Consequently, any deterministic algorithm attempting to reconstruct such polynomials would need at least Bg samples. We aim to demonstrate that even "randomized" algorithm, which succeeds with only constant probability, must still gather 3g samples. This complements our upper bound, which is established using randomized algorithm. The crucial fact we use for proving the lower bound is that all (non-zero) eigenvalues of the operator K:(4) are equal to one. This fact follows from the addition
theorem presented Theorem 3, i.e., if { (z 1/. Uf..,ee...,,.....
/ }
is an orthonormal basis of He (8d-1),
then for any function f L (Sd-1).
7
(cl)) [](()) (0)
Σ ald 0
E (
PP ((o,w)) [(w)
Σ
Why 1 1 (10)
Theorem (Lower Bound). Consider an error parameter e > O, and any (possibly randomized)
algorithm that solves Problem with probability greater than 1/10 for any input function and makes at most r (possibly adaptive) queries on any input. Then 2 Pq.d.
We describe distribution on input function f on which any deterministic algorithm that takes r < Bad samples fails with probability 2 9/10. The theorem then follows by Yao's principle.
Hard Input Distribution. For integer E < consider an orthonormal basis of He (Sd-1) and thea) Let Y, Raa, - He (sd-1) be the quasi-matrix with sf as its 5th column, i... [Y: ul(o) En u/(o) for any Roe. and gd-1 Let 200)) E Roo.d, geoos,dl (()) (()) E RoLd..... v(a) E R be independent random vectors with i.idd. Gaussian entries: v1() ~ (0, 1). The random input is defined to be En_ Ye. v(()) In other words, f El Yt (0) is random linear combination of the eigenfunctions of xc)(). We prove that accurate
denote it by { yi

---

reconstruction of f drawn from the aforementioned distribution yields accurate reconstruction of random vectors (0))) v(1) ......'e r(9), Since each u(t) is .e..-diimnsional, this reconstruction requires e(Ei (r,d) = (), samples, giving us lower bound for accurate reconstruction of f.
First we show that finding an f(() satisfying the condition of Problem 2 is least as hard as accurately finding all vectors v(0), v()) ..... v/9). The following lemma is proved in Appendix E.
Lemma 6. If a deterministic algorithm solves Problem 2 with probability at least 1/10 over our random input distribution f = -i
f= Yt -v(e), then with probability at least 1/10, the output of the algorithm fla) satisfies fef f( v for all integers 8
Finally, we complete the proof of Theorem by arguing that if f(()) is formed using less than Bg.d queries from f. then oi Y f(s) - -
1 || > 0 with good probability, thus the bound of Lemma cannot hold and flo)) cannot be a solution to Problem 2. Assume for sake of contradiction that there is a deterministic algorithm which solves Problem 2 with probability 2 1/10 over the random input o.()) that makes r = 3q.d d 1 queries on any input (we can modify an algorithm that makes fewer queries to make exactly Be,d 1 queries). For every E Sd-1 and integer S q define 2 = (o))), (/(())((().... v5(0),... Vor.a (o)] Also define == [ug, u ngg E &R and v E R as = (v(0), (1) ..... ol Additionally, define the quasi-matrix Y == [Yo................. IY...... .y).
f =i== Y
Using the above notations and the definition of the hard input instance f. each query to f in fact a query to the random vector v in the form of ((a = (u,, v). Now consider deterministic function Q, that is given input V E gixs,ee (for any positive integer i) and outputs Q(V) E R such that Q(V) has orthonormal rows with the first rows spanning the rows of V. 01, 2i.0.....; .......,... Or E gd-1 denote the points where our algorithm queries the input f, for any integer i E [r], let Q' be orthonormal matrix whose first i rows span the first queries of the algorithm, i... Q' = Q (Iuo, Moz...., 1T, Since the algorithm is deterministic, Q' is deterministic function of input v. The following claim is proved in [AKM* 19]:
Claim (Claim 23 of [AKM* 19]). Conditioned on the queries
f(o,) for <
By the variable IQ" (3n.a) ) is distributed as (0, 1).
Now using Claim we can write,
Pr
Yr f(a)
< Sr @ - =(Qrr=-ju) = Pr Q'u==Q''''' Q"Y
Bala
0=
E Pr QQ 1Q
(Qr - f(c)]
[(01), ((oi).....((.))))
where the expectation in the last line is taken over the randomness of (()).....,),,, Conditioned on f(r1)..... f(a,). [Q ( (Bq,a) is fixed quantity because the algorithm determines fl given the knowledge of the queries f(a))...... ;().)..(((.)) ((r)). Furthermore, by Claim 3, IQ" (B,.a) is This implies that,
random variable distributed as (0, 1). conditioned on
Pr IQ-h
u) (3y.d)
Q²Y f(m) (8y.) f(a.)..... ................ f(a,)) =0. 0.
Thus, Pr Ee ((() Y 7())
2 0 = E((o))... [0] 0. However, we have assumed
that this algorithm solves Problem 2 with probability at least 1/10, and hence, by Lemma 6, 0] /00. This is a contradiction, yiclding Theorem 6.
Pr El ll((0) Y f(a)
5 Numerical Evaluation
Noise-free Setting. For a fixed q, we generate a random function ((o) = ce v)) where
v ~ ( and o 's are i.i.d. samples from (0, 1). Then, f is recovered by running Algorithm with random evaluations of on gd-1, Note that | f-fllll -fllls 0 since f E H(44 (gd)), thus,

---

002
10
602
10
2400
number of samples
sac
523
num Nember of sasples
2000
2000
440
1760
1750
300
2440
0.5
(440
200
24
1120
0,4
1120
200
900
100
degree
degree q
degree
degree
(a) 3
(b) d
Figure I: (Left) Empirical success probabilities of Algorithm I varying the number of samples and the degree spherical harmonic expansion g. (Right) The dimension Bq.d of the Hilbert space ((d-)) as a function of q when (a) and (b) d = 4, respectively.
noise scale 10--6
1.0
avg. Eenorm of errors
Jif
noise scale 1p-6
noise Scale 10-6
1.0
moise scale 10
number of samples
1LS
0.8
3700
avg. Egrnrrm of
360
1440
280
10-1
1120
200
number
112
022
10-11
0.0
200
440
0.0
900 2400 -00
degree 4
number of samples 8
degree 9
number of samples 8
(a) d 3
(b) d d=
Figure 2: (Left) Empirical success probabilities of Algorithm in presence of additive noise, where "success" means the error's energy is below the noise level [] fllsa-1 S [|l]|a--1. (Right) The error's norm (If(s) - -f flls-- as a function of when (a) d = 3 and (b) = 4, respectively.
as shown in Theorem Algorithm can recover f "exactly" using s ( a log By,a) evaluations. where Bq. is the dimension of the Hilbert space ((()) We predict f's value on random test set on gd-1 and consider the algorithm fails if the testing error is greater than 10 12 We count the number of failures among 100 independent random trials with different choices of d E (3, 4}, q E (5.... 22), and E (40, ................................ 2400). The empirical success probabilities for d and are reported in Fig. (a)) and Fig. (bb), respectively. Fig. illustrates that the success probabilities of our algorithm sharply transition to as soon as the number of samples approaches s R 3y.d for wide range of and both d = 3, 4. These experimental results complement our Theorem along with the lower bound analysis in Section 4 and empirically verify the performance of our algorithm.
Noisy Setting. We repeated experiments in the presence of an additive noise which is linear combination random spherical harmonics of degrees + 1 to 2q. More precisely, we let the noise be n(o)= = Eingoe ce PP o)) E (((2))(Sd)) \ H/1) (gd-1) for eg's that are i.i.d. samples from \( 1). We then re-scale the noise to have norm [|m]|ss-- 10-6, Furthermore, the function f is defined as before, and f recovered by Algorithm with random evaluations of f - on Sd-1 The heat-maps in Fig. are generated by considering an instance of our algorithm as "success" if the error's energy is below the noise level,
< ||r||ga 10-6, The success probability transitions less sharply than the noiseless setting but the shift of probabilities starts at B..g samples.
((()) -
ga-1
6 Conclusion
We studied the problem of robustly recovering spherical harmonic expansion of function defined on the sphere. The number of function evaluations needed recover its degree-q expansion the dimension of spherical harmonies of degree at most q, up to logarithmic factor. We develop simple yet efficient kernel regression- evaluating the function on uniformly sampled points on the sphere. Unlike the prior results on fast spherical harmonic transform, our algorithm works efficiently using nearly optimal number of samples in any dimension. We believe our findings would appeal to the readership of the community.
algorithm recover degree-q expansion of the function by only
10

---

References
[AH12]
Kendall Atkinson and Weimin Han. Spherical harmonics and approximations on the unit sphere: introduction. Springer Science & Business Media, 2012.
[AKK*20]
Thomas D Ahle, Michael Kapralov, Jakob BT Knudsen, Rasmus Pagh, Ameya Vel- ingker, David Woodruff, and Amir Zandieh. Oblivious sketching of high-degree polynomial kernels. In Symposium Discrete Algorithms (SODA), 2020.
[AKM* 17] Haim Avron, Michael Kapralov, Cameron Musco, Christopher Musco, Ameya Vel- ingker, and Amir Zandieh. Random Fourier features for kernel ridge regression: Approximation bounds and statistical guarantees. In International Conference on Machine Learning (ICML), 2017.
[AKM+19] Haim Avron, Michael Kapralov, Cameron Musco, Christopher Musco, Ameya Vel- ingker, and Amir Zandieh. A universal sampling method for reconstructing signals with simple fourier transforms. In Symposium the Theory of Compuling (STOC), 2019.
[BKM*23] Karl Bringmann, Michael Kapralov, Mikhail Makarov, Vasileios Nakos, Amir Yagudin, and Amir Zandieh. Traversing the fft computation tree for dimension-indeeende sparse fourier transforms. In Proceedings of the 2023 Anmual ACM-SIAN SIAM Symposium on Discrete Algorithms (SODA), pages 4768-4888 4845. SIAM, 2023.
[CKPS16] Xue Chen, Daniel M Kane, Eric Price, and Zhao Song. Fourier-sparse interpolation without frequency gap. In Foundations of Computer Science (FOCS), 2016.
[CLL15] Pui Tung Choi, Ka Chun Lam, and Lok Ming Lui. FLASH: Fast landmark aligned spherical harmonic parameterization for genus-O closed brain surfaces. SIAM Journal on Imaging Sciences, 2015.
[CP19]
Xue Chen and Eric Price. Active Regression via Linear-Sample Sparsification. In Conference on Learning Theory (COLT), 2019.
[EEHM17]
Michael Eickenberg, Georgios Exarchakis, Matthew Hirn, and Stéphane Mallat. Solid harmonic wavelet scattering: Predicting quantum molecular energy from invariant descriptors of 3d electronic densities. Advances in Neural Information Processing Systems, 30, 2017.
[EMM20] Tamás Erdélyi, Cameron Musco, and Christopher Musco. Fourier sparse leverage scores and approximate kernel learning. Neural Information Processing Systems (NeurIPS), 2020.
[GCL+22] Jan Gerken, Oscar Carlsson, Hampus Linander, Fredrik Ohlsson, Christoffer Petersson, and Daniel Persson. Equivariance versus augmentation for spherical images. In International Conference on Machine Learning, pages 7404-7421. PMIR, 2022.
[GMMM21] Ghorbani, S Mei, T Misiakiewicz, and Montanari. Linearized two-layers neural networks in high dimension. The Annals of Statistics, 2021.
[HZA22] Insu Han, Amir Zandieh, and Haim Avron. Random Gegenbauer Features for Scalable Kernel Methods, In International Conference on Machine Learning, pages -330-8358. PMLR, 2022.
[HZL+22] Insu Han, Amir Zandieh, Jaehoon Lee, Roman Novak, Lechao Xiao, and Amin Karbasi. Fast neural kernel embeddings for general activations. Advances in neural information processing systems, 35:35657-35671, 2022.
[KFR03]
Michael Kazhdan, Thomas Funkhouser, and Szymon Rusinkiewiez. Rotation invariant spherical harmonic representation of 3 d shape descriptors. In Symposium Geometry Processing. 2003.
[KKS97]
Mare Kamionkowski, Arthur Kosowsky, and Albert Stebbins. Statistics of cosmic microwave background polarization. Physical Review D, 1997.
11

---

[KVZ19]
Michael Kapralov, Ameya Velingker, and Amir Zandieh. Dimension-indeeenden sparse fourier transform. In Proceedings of the Thirtieth Annual ACM-SIAA Symposium on Discrete Algorithms, pages 2709-2728. SIAM, 2019.
[Lan12]
Serge Lang. Real and functional analysis. Springer Science & Business Media, 2012.
[LMPI3]
Mu Li, Gary L. Miller, and Richard Peng. Iterative row sampling. In Foundations of Computer Science (FOCS), 2013.
[LP61]
Henry J Landau and Henry Pollak. Prolate spheroidal wave functions, Fourier analysis and uncertainty- II. Bell System Technical Journal, 1961.
[LP62] Henry Landau and Henry O Pollak. Prolate spheroidal wave functions, fourier analysis and uncertainty- III: the dimension of the space of essentially time-and sand-limite signals. Bell System Technical Journal, 1962.
[LWL-21] Yi Liu, Limei Wang. Meng Liu, Yuchao Lin, Xuan Zhang, Bora Oztekin, and Shuiwang Ji. Spherical message passing for 3d molecular graphs. In International Conference on Learning Representations, 2021.
[MNY06]
Ha Quang Minh, Partha Niyogi, and Yuan Yao. Mercer's theorem, feature maps, and smoothing. In Conference on Learning Theory (COLT), 2006.
[Mor98]
Mitsuo Morimoto. Analytic functionals on the sphere. American Mathematical Soc., 1998.
[Pen30]
WO Pennell. A generalized Fourier series representation of function. The American Mathematical Monthly, 1930.
[RT06]
Vladimir Rokhlin and Mark Tygert. Fast algorithms for spherical harmonic expansions. SIAM Journal on Scientific Computing (SISC), 2006.
[SA22]
Paz. Fink Shustin and Haim Avron, Semi-Infinite Linear Regression and Its Applications. SIAM Journal on Matrix Analysis and Applications (SIMAX), 2022.
[SP61]
David Slepian and Henry O Pollak. Prolate spheroidal wave functions, Fourier analysis and uncertainty-I. Bell System Technical Journal, 1961.
[SS00]
Paul N Swarztrauber and William Spotz. Generalized discrete spherical harmonic transforms. Journal of Computational Physics, 2000.
[ST02]
Reiji Suda and Masayasu Takami. A fast spherical harmonics transform algorithm. Mathematics of Computation, 2002.
[Wei95]
DR Weimer. Models of high-latitude electric potentials derived with a least error fit of spherical harmonic coefficients. Journal of Geophysical Research: Space Physics, 1995.
[Wer97]
Robert A Werner. Spherical harmonic coefficients for the potential of a constant-density polyhedron. Computers & Geosciences, 1997.
[Wil12]
Virginia Vassilevska Williams. Multiplying matrices faster than Coppersmith- -inngggdd. In Symposium on the Theory of Computing (STOC), 2012.
[WZ22]
David Woodruff and Amir Zandieh. Leverage score sampling for tensor product matrices in input sparsity time. In International Conference on Machine Learning. pages 23933-23966 PMLR, 2022.
[XRYO]] Hong Xiao, Vladimir Rokhlin, and Norman Yarvin. Prolate spheroidal wavefunctions, quadrature and interpolation. Inverse problems, 2001.
[2DK+22] Larry Zitnick, Abhishek Das, Adeesh Kolluru, Janice Lan, Muhammed Shuaibi, Anuroop Sriram, Zachary Ulissi, and Brandon Wood. Spherical channels for modeling atomie interactions. Advances Neural Information Processing Systems, 35:8054- 8067, 2022.
12

---

[ZHA 21] Amir Zandieh, Insu Han, Haim Avron, Neta Shoham, Chaewon Kim, and Jinwoo Shin. Scaling Neural Tangent Kernels via Sketching and Random Features. In Neural Information Processing Systems (NeurIPS), 2021.
13

---

A Properties of Gegenbauer Polynomials and Spherical Harmonics
In this section we prove the basic properties of the Gegenbauer Polynomials as well the Spherical Harmonics and establish the connection between the two. We start by the direct sum decomposition of the Hilbert space L2 (2 in terms of the spherical harmonics,
Lemma 2 (Direct Sum Decomposition of L²($d-1)). The family of spaces He (St-1) yields Hilbert space direct sum decomposition I,2 (sa-1)- θ He (st-1): the summands are closed and pairwise orthogonal, and every f E L² (sd-1) is the sum of converging series (in the sense of mean- sqcare convergence with the -o-morm defined in Eq. (2)),
ft
t=
where Se He 18 are uniquely determined functions. Furthermore, given any orthonormal
basis { 1 .s...., ..n. Vana) of (Sa-1) we have ft Dj=í (f.
Proof. This in fact a standard result. For example, see [Lan12] for a proof.
Now we show that the Gegenbauer polynomials and spherical harmonics are related through the so called the addition theorem,
Theorem (Addition Theorem). For every integer l 2 0, if { uu ... Waae } is orthonormal -f.
basis for He (SQ-)), then for any a, w E gd-1 we have
aiedd
Pá (10, w) ( vl(oo) v5()).
Proof. The result can be proven analytically, using the properties of the Poisson kernel in the unit ball. This is classic and the proof can be found in [AH12, Theorem 2.9].
Next show that the Gegenbauer kernels can project any function into the space of their corre- sponding spherical harmonics,
if Lemma (A Reproducing Kernel for He (Sd-1)). For every f L2 if Eizo f fe is the unique decomposition of over θi, Ho (gd-i) as per Lemma 2, then ft is given by
fela) = ag,d
E (
[f()) ((a, w))] for sd-1
Proof. This is classic textbook result, see [Mor98].
Now prove that the Gegenbauer kernels satisfy the reproducing property for the Hilbert space Ho (gd
Lemma (Reproducing Property of Zonal Harmonics), Let (-) be the Gengenbauer polyyoomaarollyymmial of degree E in dimension d. For any x, y E sd-1,
Pa((x,v)) and Eun [P] ((x, w)) Pf (Ky, w))]
Furthermore, for any +l:
E B
Pá ((x,))) Ph (1,, w))]
14

---

Proof. This result follows directly from the Tunk-Hecke formula (See [AH12]). However, we provide another proof here. First note that for every z E g the function Pj ((x, -)) E H4 (su-1). Therefore the first claim follow by applying Lemma on function f(o) = Pa ((),o)) which also satisfies ft- = f. On the other hand, P ((v.)) E He (S for every y E gd-1 Thus, for l' # L, using the fact that spherical harmonics are orthogonal spaces of functions, Pd (W. .)) 1 He (su-1). which gives the second claim.
Next we prove that the kernel operator defined in Definition 2 is in fact a projection operator,
Claim 1. The operator K(g) defined in Definition 2 satisfies the property (K!()) x(()
Proof. For every f E L (sd-1) and every o E sd-1, using Definition 2 we have, [(x)9) 2 f (0)
4 apd
(xo9) f, P. ((o,)1)) ge-1
Sd-1
app
R 1)
Pd (10, ((o,w)).
E 1-
[Pd (() f(r)]
8=
E T
ʃ(T)
E w-
[ w)) P4 ((T, [1))
P=
Σ aed
E (
[S(T). P1 ((a,7))
E=0
[x(0
5]
(a),
where the fourth line above follows from Lemma 1. This proves the claim.
B Reducing the Interpolation Problem to a Least-Squares Regression
In this section we show that our spherical harmonic interpolation problem, i.e., Problem 2, can be solved by approximately solving -quurrss problem as claimed in Claim 2. We start by showing that for any function f E 1 xC(") f gives its low-degree component. More precisely, let = Ee fe be the decomposition of f over the Hilbert sum Di He (Sd-1) as per Lemma 2. Now
if we let kC be the kernel operator from Definition 2 and if ( vb.... Vea is an orthonormal basis for He (sa-1), then by Theorem 3 we have,
(()1) 5 (0)
{](0)
Σ
ard
E 1
(lo,w))
t=0
w
asa
od
E [w)
(88-)) 1)
t=
j=1
5o( 1
E
[s(w) (oo))./ff [f(ww]
aid
(S (.,5(1))))r v
f=1
EnO = (o),
15

---

where the the second line above follows from Theorem 3, the fourth line follows from Eq. (2), and the last line follows from Lemma 2. This proves that the low-degree component f( kg f.
Claim 2. For any f E L (sa-1). any integer 2 O, and any C' C22 1, ifg L2 (S6-1) satisfies,
|K' g- f :lle- SC. min 1E12(
| $1
and if we let f(9) = Kg)) S. where Kg)) is defined per Definition 2, then the following holds
||r_000 g- f(
gh-1 13- [2., S -1) -1)-11500) - 8 gM-1
Proof. First, note that g² = f is an optimal solution to the least-squares problem in Eq. (7). Thus we have,
|r -
1|33 su
= =11x00 -- - 11(3- = =1|p(0- -5113-
min SE7P(8 1)
Next, we can write,
k<(1) - -
= | 3 x(( ʃ + (K! ʃ 1) 1
120-1
38-1
| (0 f) + (K(T²) f -- f) 1
gd--1
|| (ā - -f) 2
+ [x(() f -f(
ga
1g
| f(
2
+ -||550) - 8|180-1
ge-1
where the third line follows from the Pythagorean theorem because X/ (g ( H() (sd-1) while x:())
Eeso fe, thus (K(() f f- -) 1 H(w) (5U-1). Combining the two equalities above ga S C' mmiigge 1) |c g ga given in the statement of
with inequality |x/90) 1
f11
the claim, proves Claim 2.
C Approximate Regression via Leverage Score Sampling
In this section we ultimately prove our main result of Theorem 4. We start by proving useful properties of the leverage function given in Definition 3. First, show the fact that the leverage function can be characterized in terms of a least-squares minimization problem, which is crucial for computing the leverage scores distribution. This fact was previously exploited in [AKM t 17] and [AKM* 19] in the context of Fourier operators.
Lemma 4 (Min Characterization of the Leverage Function). For any w E @d-1 let Ty (w) be the leverage function (Definition 3) and define @u E (Sd by @w(a) = Ei=
and Pd (la, w)). gd-1
We have the following minimization characterization of the leverage function:
Ty (w)
11.. 1, S.S. K(9) g Φw (GEL)(
min
We remark that this lemma is in fact an adaptation and generalization of Theorem 5 of [AKM+ 19]. We prove this lemma here for the sake of completeness.
Proof. First we show that the right hand side of Eq. (9) never smaller than the leverage function in Definition 3. Let gh, E L (gd-- be the optimal solution of Eq. (9) for any gd-1, Note that the
16

---

optimal solution satisfies / the = ow. Thus, for any function ʃ E (2(S4-1), using Definition 2, we can write
2
[ r|
(w))
= Σ t=0
OcE
E
[P ((a, w)) f(o)
a-21(8-1)
Wshre f'(ga-1 1
(cll) Per 1) ,)
Koha xC()))
f)
2
(because K2 self-adjoint)
8/80111
11 lzg-1
|
112
(by Cauchy-Schwand ineuuuuulauiii inequality)
184-1
Therefore, for any f 1,2(8(-1) with
|x(0)) r||
> O, we have
84-1
[ f|
(w))
I 11 12a-..
(11)
|
f 11
We conclude the proof by showing that the maximum value attained. First, we show that the optimal solution The of Eq. (9) satisfies the property that k-(9) Fwr Di. Suppose for the sake of contradiction that x(()) Iw # Iw- this case, Claim implies that,
(c(()) Fir - g.m -xc)) I.. KY = x(( g. xc(a) Tu = 0. xc(9) (K2
Thus, the function x((1) Iu satisfies the constraint of the minimization problem in Eq. (9). Now, using the above and the fact that k(()) self-adjoint we can write,
(Keg) Fen k-(9) Ier gin) ga-1 = ( kcla) (NOS) gi sin) ga-1 ) 0.
This shows that k(4) su 1 (Kel)) yt Itu) hence by Pythagorean theorem we have,
| Ig- =
| =1| Fer + || gi. Iir gu
>
| Fer _ gg
= ||al
which is in contrast with the assumption that 9 is the optimal solution of Eq. (9). Therefore, our claim that x-(fi In = the holds.
Now, we show that for f = Iw the maximum value inequality Eq. (11) is attained. For any w E Sd-1 we have the following
(x(()) /) (w) = (KC9)) Dio sa-1 = (Dwr Fulga-n = -11g
On the other hand we have [K((0) /|/3 11 11 Thus, [c()) ga-1 I I which implies that Ty (w) l1 1i and thus proves the lemma.
(AC()) r r](w)] (w))
Next we prove that the leverage function constant.
Lemma (Leverage Function Constant), The leverage function given Definition is equal to (w)
Σ1
azddd
for every E Sd-1
Proof. First prove that Ty(w)= (==0 T.2 (8d-1) be defined as (a)) == Eg
g using the min-characterizzaton. If we let due pf ((a, w)), then by Definition 2, for every a E gd-1
ard
an
[
17

---

we can write,
4 Σ
(( Qu (0)
atd
E 1-
[P] ((o, w)) (w(w)
==
9
ag
Pá ((o,))))) Py ((rw)))
[Sd-1]
t=
=0=0
afg,
Pá ((o, w/) = (w)
(12)
2=
where the third line above follows from Lemma 1. Therefore, the test function == &w satisfies the constraint of the minimization in Eq. (9), i.e., Kid 6n. Thus, Lemma 4 implies that,
Tig (w) S |
1 1:
Σ
S-
where the equality above follows from Lemma along with Eq. (2). This establishes the upper bound on the leverage function that we sought to prove.
Now, using the maximization characterization of the leverage function in Definition 3, we prove that Again, we consider the same test function g = &w and write.
(w) 2 Ei
aed
2 12)00) Owe ga-1
[|xc58 dwl (w))
dw
1ow (w))
1 gd
Pd (ww, 2
and 8
Σl=
(1)
3
ag5
E
S
S4-11'
8=
where the first and second line above follow from Eq. (12) and Lemma I, respectively, Therefore, the max characterization the leverage function Definition implies that,
Ty (w) 2
|
[<
(w))
atd Dd-1
This completes the proof of Lemma 5 and establishes that (w) is uniformly equal to >% 2x
ard
To prove Theorem 4, we need to use prior results about solving linear systems in continuous setting via leverage score sampling. In particular, we use Theorem 6.3 from [CP19], which is restated below,
Theorem 7 (Theorem 6.3 of [CP19]). Consider any dimension linear space F of functions from a domain G to IR. Let be distribution over and f some function from G to R. Also, define the norm with respect to D of any function h G + R [|Ail
Ee-D[[h(z)[]] and let y arg minner || fl Fix any distribution D' over G and let
KD = supred SuPREr (x)
D(z)
[h(2]
D(ze) For i..dd sample query points I1, T2,... Ea D' and weights wi ..D for i E [s), let the weighted ERM estimator f. be defined as F. == arg min Σi= [h(x,) - 1 For any
>> 0 and a sufficiently large fixed constant if the number of queries 2 c (Ko logn Kp
then the weighted ERM estimator J. satisfies,
I1
Pr
S s Ilf w11 21 10-4
18

---

Now we are ready to prove Theorem 4. Our approach is to apply Theorem 7 to the space of degree spherical harmonics H H(9 (su-1) and use the fact that the leverage scores distribution of the operator are uniform on the unit sphere gd-1
Theorem 4 (Approximate Regression via Leverage Function Sampling). For any ε >> 0, let (Bea log a + samples on g Define the quasi matsi-matrix P II" - (2( as follows, for every Rd:
Ao.d = ) ).for sufficiently large fixed constant e, and let æ1, 22,...,, Is be i.d. uniform
PP)[())= 8=
Off,d
Vj Pá (1)),0)) for E gd-1
Sd-1
j=1
Also let f E I2 be vector with fj :- V M(x;) for 1,2,................ and let P be the adjoint of P. 9 is optimal solution the least-square. problem arg min PP' f fllg. then with probability at least 10-4 the following holds,
| 9- -1 10-1 S (1++))
1C/0 g / gd-1
min GEL2(04-1)
Proof. We prove this theorem by invoking Theorem 7. To do so, we first let the space F of function from Sd-1 E be J := H(4) (gd-1). is clear that the space of spherical harmonics linear space of functions because of the existence of the kernel operator xc(0) which is projection operator onto H(9) (gd-1), so F satisfies the first precondition of Theorem 7. Additionally, the dimension of this space of functions is n aq
Also, let D be uniform distribution over the unit sphere 8d-1, For this distribution, the norm defined in Theorem 7 satisfies ||A]| = |sd-1| -1 ||A|| D. for any E (2( Moreover, D' be a uniform distribution over the unit sphere gd-1 as well.
Now we show that for these choices of F, D, D', the condition number KD defined as per Theorem 7 equal to Bg.d- We can write,
D(x)
[xx)]
Kp = sup sup xe3u-1 hEF
D'(x)
[1xx)12
sup sup |gd-1|
[
=184-1 KeT
2
(x)
1gd-11
sup
sup
acge-1 (8
((())
gar-
gd-11
sup
(x)
x
3q.ds
where the second line above follows from the fact that D, D' are both equal to the uniform distribution over gd-1 and | = -- ||A|| The third line above follows from the fact that any function
h E H(0) (s can be expressed as h x(0)) g for some 9 E ( because x<(0) is the projection operator onto 7((0 (Su-1). The fourth line follows from Definition and last line follows from Lemma 5.
19

---

Finally, in order to invoke Theorem 7, we can write that the weighted ERM estimator fo is equal to the following.
Fx = arg min REF
Wi
[h(2z)) - (xx))
||<580 9) (Iri) -((rr))
arg min L
wi
iml
k-(0) 9 (xi) - fi 2
arg min 9EL-(84-1)
v | fl
t=1
arg min (B4-()4-1)
where the third line above uses the definition of fi = t/ (xx)) and the last line follows from the definition of adjoint of the quasi-matrix P. Therefore, the theorem follows by invoking Theorem 7.
D Efficient Algorithm for Spherical Harmonic Interpolation
In this section prove main theorem about spherical harmonic interpolation algorithm. Theorem 5 (Efficient Spherical Harmonic Interpolation). Algorithm 1 returns function (g(() such that, with probability at least 1 - 10-4:
| p()) $5.11500) - 1 where f(a) k(())
Suppose we can compute the Gegenhauer polynomial (t) at every point E [-1, in constant time. Then Algorithm queries the function at O (Bu,d log 9q.d + Byd points the sphere gd-1 and runs in O d 5²) time. This algorithm evaluates y(o) in Old 8) time for any e sd-1
Proof. First note that the random points W1, W2 ., Ws in line 3 of Algorithm are i.i.d. sample with uniform distribution on the surface of go-i Therefore, we can invoke Theorem 4. More specifically, if we let P be the quasi-matrix defined in Theorem 4 corresponding to the random points W1, W2...., we sampled line 3 and if we let f be the vector of function samples defined in line 5 of the algorithm, then with probability at least 10-4, any optimal solution the following least-squares problem
g arg
min (8d-(
P - -f1
(13)
satisfies the following.
f 11001 S (1++)) |K g
1c(0) g get (0u-1)
f 2
min
9
(14)
Now note that the least-squares problem in Eq. (13) has at least one optimal solution ğ which is in the eigenspace of the operator PP* More specifically, there exists vector RR such that ý P is an optimal solution for Eq. (13). Therefore, we can focus on finding this optimal solution by solving the following least-squares problem
arg min zER*
P' f11
and then letting P. z. This is guaranteed to be an optimal solution for Eq. (13), thus it satisfies Eq. (14). We solve the above least-squares problem using the kernel trick. In fact we show that P p is equal to the kernel matrix K computed in line of Algorithm 1. To why, note that for any
20

---

j.j [s] we have,
4
4 1),
pP [P'Pl Plas
01,d
Pd (Wwi)
Σ
atd
Py ((wi.)))
Sd
V
[80-1|
1=
ga-1
Pg ((ws, Ph ((urg) ga
Rd--19 112
(=0
Of
[t ((wi, u)) Pd ((w,, m))
Sd-1
ag
Pá ((wy, wy)) Kigt
Sd-1
where the fourth line above follows from Lemma 1. Therefore, we are interested in the optimal solution of the following least-squares problem
E arg min K ffll2.
The least-squares solution to the above problem is z Kt which exactly what is computed in line 6 of the algorithm. Now note that, the function ğ P x satisfies Eq, (14). Because P. H(1) and because K(() is orthonormal projection operator into H (Sd-1), we have x(())
9 g- = P.2..
2. This together with Eq. (14) imply that,
P fll a1)
min 9E12
|// g s
1)
Now we invoke Claim with C'
=1++ 1 on the above inequality we find that,
p- zp-z- f | ga-1 S $..
ga-1
Finally, one can easily see that the function E H(19) (Sd-1) that Algorithm 1 outputs line is exactly equal P. z. This completes the accuracy bound of the theorem.
Runtime and Sample Complexity. these bounds follow from observing that:
8 d time is needed to generate wi, W2, W2,.... We in line of the algorithm. To do this, we first generate random Gaussian points in RRd and then project then onto $d-1 by normalizing them.
s² d operations are needed to form the kernel matrix K in line 4 of the algorithm.
queries to function f are needed to form the samples vector f in line 5 of the algorithm.
~ time is needed compute the least-squares solution Kt f in line 6 of the algorithm.
s d operations are needed evaluate the output function y in line 7 of the algorithm.
This completes the proof of Theorem 5.
E Lower Bound: Claims and Lemmas
In this section prove the Claims and Lemmas used in lower bound analysis for proving Theorem 6.
Claim 4. Given the random input f - Ei Ye ( generated as described in Section 4, to solve
Problem 2, an algorithm must return function f(() H(9) (gd-1) such that [l fl - 0.
21

---

Proof. Note that Problem 2 requires recovering a function f(o) E H(l0) (s*-1) such that:
f(a) f(s) | 18--1
- 113 18-11
S
(15)
where f(() ~C88 f. Using the definition of the input function f -it=o Ye vil,, write,
f(a) K
Yt v(0)
6=0
q
9
Ye Yē
E=
P==0
Σ Yi 7 f,
where the equality in the second line above follows from Eq. (10) and the addition theorem in Theorem 3, and the third line follows because the operator Y, has orthonormal columns and thus YE Yr loca [(==0). Therefore, plugging this into Eq. (15) gives.
((0) - ga-1 = | 5(0l 13--- Se- (1-11
Sl = --l1s- 0.
Lemma 6. If a deterministic algorithm solves Problem 2 with probability at least 1/10 over our random input distribution f = Σi= Ye vu(),, then with probability at least 1/10, the output of the algorithm fla) satisfies Y, f(1) v( for all integers S
Proof. By Claim 4, the output of the algorithm that solves Problem 2, satisfies [(() - l||
18-1 O.
Therefore, by orthonormality of the columns of the operator Ye, we can write,
4
;ı) Y + Yi( -) Σ Yi ((() v()
P==
22


Reviews:
[{"review_id": "4Qh89GoN6R", "reviewer": "Reviewer_BcP8", "summary": "Consider the $d$-dimensional unit sphere $\\mathbb{S}^{d-1}$, and any function $f:\\mathbb{S}^{d-1}\\rightarrow\\mathbb{R}$ defined on the sphere with bounded L2 norm $\\|\\|f\\|\\|_{\\mathbb{S}^{d-1}}$. This function $f$ can be evaluated at any point $\\vec w \\in \\mathbb{S}^{d-1}$ on the sphere, but it is expensive to evaluate $f$ so we wish to do this as little as possible.\n\nThe goal of this paper is to recover a near-optimal polynomial $\\tilde f$ approximation to $f$, where $\\tilde f$ is constrained to be a multivariate polynomial where the sum-of-degrees of each term is at most $q$. (For example, the sum-of-degrees of each term in $\\tilde f(x,y,z) = x^2y^3 + z + x^5$ is at most 5).\n\nLetting $\\beta_{d,q}$ denote the number of free parameters in any $d$-variate polynomial with terms whose sum-of-degrees is at most $q$, this paper shows that $O(\\beta_{d,q} \\log \\beta_{d,q} + \\frac{\\beta_{d,q}}{\\varepsilon})$ evaluations of $f$ chosen uniformly at random on the sphere suffice to recover such a near-optimal $f$. This is proven in three parts:\n1. The leverage function associated with this polynomial recovery problem is **exactly** the uniform distribution on the sphere\n2. The randomized linear algebra toolkit shows this near-linear sample complexity suffices\n3. A small kernel ridge regression problem can be solved to recover this $\\tilde f$ from uniform random samples on the sphere\n\nThey also prove an $\\Omega(\\beta_{d,q})$ sample complexity lower bound. Lastly, they include some synthetic experiments.", "questions": "I list minor questions, typos, and recommended edits here. These are all soft recommendations, adapt whatever you want to.\n\n1. [Line 53] Consider mentioning a failure probability here.\n1. [Line 87] Replace \"this paper\" with \"[GMMM21]\" since the language is a bit too ambiguous at a glance.\n1. [Line 173] Is this not a bijection? At least for traditional least squares regression, this seems to be a bijection? RandNLA people tend to use the former form, so it's nice to see if the RandNLA form is exactly equivalent to the form of Problem 2.\n1. [Line 213] Consider removing the $\\cdot$ between P and v? This notation feels a bit clunky and odd to me?\n1. [Line 220] Cite [AKM+] or [SA] or [CP from COLT19] here, or some other paper that draws a connection between the semi-infinite regression and kernel ridge regression. This section reads like the idea of the kernel trick here is a new contribution.\n1. [Line 226] Add \"Then,\" before \"Algorithm 1\"\n1. [Line 236] Discuss why this isn't trivial that we need $\\beta_{d,q}$ samples to learn $\\beta_{d,q}$ parameters in a polynomial. Certainly, when we're looking at interpolating a polynomial on the real line, we need at least $q+1$ samples to learn a degree $q$ polynomial exactly. This feels similar to the argument made in this lower bound, but I don't understand what the rest of the formalization is needed for.\n1. [Line 250] Replace \"accurately\" with \"exactly\" (I think that's more technically accurate?)\n1. [Lines 263-274] I got really lost here. I think that If I'd read the lower bound of AKM, then I might be able to stitch together the parts of this lower bound, but the reason any of this construction is made is unclear to me. The notation is also pretty hard to track. I don't really get what the construction is trying to get at. I don't really see why the top rows of Q need to span the queries made so far. This all feels odd to me.\n1. You don't need to do this. But, if you want to, I think you could make a lower bound of $\\Omega(\\frac1\\varepsilon)$ by following the technique in Section 5.4 of [here](https://arxiv.org/pdf/2211.06790.pdf). It'd be interesting to see if this technique holds on the sphere, because at a glance it seems like it should. It'd give an overall lower bound of only $\\Omega(\\beta_{d,q} + \\frac1{\\varepsilon})$ though, which isn't the most compelling rate.", "limitations": "N/A", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper is a nice application of well understood tools from randomized linear algebra (RandNLA) to spherical harmonics, which has not been explored by the RandNLA literature afaik. The real core strength of this paper is its novelty in connecting these two literatures in a simple and elegant way.\n\nSimplicity and elegance really are words that mark this paper. Almost everything is extremely clear and well written. There are basically no typos even! The flow of logic in the paper is very clear, the problem they solve seems important in the spherical harmonics literature (though I'm no expert in that domain), and the proofs are even pretty clean. To that last point, the theorems in this paper are proven either by very clean techniques that are now standard in the RandNLA literature, or by using well established classical facts about polynomials and spherical harmonics. I verified a decent amount of the math, and it all struck me as rather nice and clean.\n\nSince I'm not an expert in spherical harmonics, I can't perfectly speak to the significance of the result. Taking the authors' words at face value, it seems that prior work in spherical harmonics did not use such simple algorithms (Kernel Ridge Regression) and did not achieve optimal sample complexities.\n\nThe paper may not carry a huge amount of new technical ideas to get their sample complexity and simple algorithm, but this result would only exist if someone who knew enough about both spherical harmonics and RandNLA decided to sit down and figure out if this all works together. For that view of novelty, in addition to the simplicity and elegance of the paper, I recommend accepting this paper.\n\nThe experiments are perfectly fine for a theoretical paper. Nothing particularly strong about them, but nothing I'm left looking for.", "weaknesses": "The novelty, simplicity, and elegance is strong. But, the proof techniques are not especially novel. The proof techniques are standard relationships between operators and algorithms as explored by the RandNLA community for a good few years now.\n\nThe only proof that doesn't seem directly tied to something already proven in the RandNLA literature is the proof that the leverage function for the regression problem is uniform on the sphere. Even then, this claim is pretty intuitive, since the problem statement is rotationally invariant. There's no reason for a sampling algorithm to care more about one point of the sphere than another. (This is in contrast to learning on an interval, where an algorithm may prefer to sample near the edges of the interval.) So, even this most novel proof is not terribly surprising, and the proof is short and simple.\n\nThat said, I hesitate to really call this a \"weakness\". This isn't clearly a downside of the paper. While it's certainly nice for a paper to overcome new technical issues and proof difficulties, it's also nice for a paper to show that spherical harmonics smoothly fit into the RandNLA framework.\n\nAll this is to say that the technical weight of the paper truly is in understanding two literatures and connecting them. The effort to connect the literatures seems to be low, but that's given that someone actually understands both literatures well. I'd say this overall comes out as a slight strength for the paper, but it's certainly a tradeoff in the reviewing process.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Many thanks for your valuable feedback on our paper. Below we answer your questions.\n\n\n>`Discuss why this isn't trivial that we need \\beta_{d,q} samples to learn \\beta_{d,q} parameters in a polynomial.`\n\n\nYou're absolutely correct in pointing out that there are $\\beta_{d,q}$ degrees of freedom, thus any deterministic algorithm that reconstructs such polynomials needs at least $\\beta_{d,q}$ samples. Our lower bound proof is showing that even a \"randomized\" algorithm that succeeds with only a constant probability, needs to take $\\beta_{d,q}$ samples. Since our upper bound is established using a randomized algorithm, it was crucial to complement it with a randomized lower bound to create a well-rounded and balanced analysis of the problem.\n\n>`I think that If I'd read the lower bound of [AKM+], then I might be able to stitch together the parts of this lower bound, but the reason any of this construction is made is unclear to me. The notation is also pretty hard to track. I don't really get what the construction is trying to get at. I don't really see why the top rows of Q need to span the queries made so far. This all feels odd to me.`\n\n\nWe apologize for any confusion caused by the clarity of our lower-bound section. The limited space may have contributed to the lack of detailed explanations.\nTo provide further clarification, our hard instance is based on a random vector ${\\bf v}$ following an isotropic Gaussian distribution in dimension $\\beta_{d,q}$. In line 273, our aim is to demonstrate that if an algorithm reconstructs a function $\\tilde{f}^{(q)}$ using only $r <  \\beta_{d,q}$, then even after conditioning on the samples observed by the algorithm, vector ${\\bf v}$ will still possess at least one degree of freedom and will not be entirely deterministic. This intuitively holds true because ${\\bf v}$ consists of $ \\beta_{d,q}$ independent random Gaussian entries. Thus, when conditioning on $r < \\beta_{d,q}$ samples taken by an algorithm, the conditional value of ${\\bf v}$ will remain random, allowing us to extract at least one Gaussian random variable from it using an orthonormal transformation denoted as $Q^r$. \nWe will provide additional details and improve the notation for better readability.\n\n\n> `Regarding minor typos:` \n\nThanks for pointing them out. We will address them all and implement your recommended edits.\n\n\n"}, {"review_id": "uIOLSFMmRP", "reviewer": "Reviewer_7fua", "summary": "To develop kernel regression based algorithm to recover degree-q expansion of f \\in L^2(S^{}d-1| - by only evaluating on uniformly sampled points on S^{d-1}.\n", "questions": "-", "limitations": "-", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The ideas used in this paper is deeply technical and ideas are complicated. It first re-formulates the problem as least squares regression and then uses the sampling based on leverage scores to decide which samples to pick. Then the rest of the paper is proving the estimated bound is satisfied.", "weaknesses": "motivate the problem - how it could be applied in real life\n\n** after rebuttal -- thank you for providing all these papers regarding contributions of this workstream in terms of application side. I am happy to increase my score regarding soundness of the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> `Motivate the problem - how it could be applied in real life?`\n\nMany thanks for your valuable feedback on our paper. \nWe understand the importance of motivating the problem and highlighting its real-life applications. In response to your comment, we would like to elaborate on the practical significance of the problem we address and its potential applications.\n\nOur main focus is recovering (unknown) functions defined in the sphere, a critical task in scenarios where rotational invariance is a fundamental property. In real-life applications, such as machine learning, this property becomes very important as a foundational requirement for modeling 3D point clouds. Notable examples that appeared in top-tier machine learning conferences include molecular/atom systems, where understanding the underlying functions within a spherical context can significantly enhance predictive modeling and simulation accuracy [1,2,3,4]. Other examples are in the field of computer vision, specifically in the recognition, classification, and reconstruction of 3D objects [5,6,7,8].\n\n[1] Eickenberg, Michael, et al. \"Solid harmonic wavelet scattering: Predicting quantum molecular energy from invariant descriptors of 3D electronic densities.\" NeurIPS 2017\n\n[2] Frank, Thorben, et al. \"So3krates: Equivariant attention for interactions on arbitrary length-scales in molecular systems.\" NeurIPS 2022\n\n[3] Zitnick, Larry, et al. \"Spherical channels for modeling atomic interactions.\" NeurIPS 2022\n\n[4] Liu, Yi, et al. \"Spherical message passing for 3d molecular graphs.\" ICLR 2022\n\n[5] Gardner, James, et al. \"Rotation-Equivariant Conditional Spherical Neural Fields for Learning a Natural Illumination Prior.\" NeurIPS 2022\n\n[6] Melnyk, Pavlo, et al. \"Steerable 3D spherical neurons.\" ICML 2022\n\n[7] Gerken, Jan, et al. \"Equivariance versus augmentation for spherical images.\" ICML 2022\n\n[8] Shakerinava, Mehran, and Siamak Ravanbakhsh. \"Equivariant networks for pixelized spheres.\" ICML 2021\n"}, {"review_id": "H8Jujh3xCh", "reviewer": "Reviewer_6ewG", "summary": "The paper studies the approximation of a function $f \\in L_2(\\mathbb{S}^{d-1})$ from its evaluations, via a degree-q spherical harmonic expansion. To this end, an efficient kernel regression based algorithm is proposed which recovers such a degree-q expansion of f, from the evaluations of f on $\\mathbb{S}^{d-1})$. In particular, the number of evaluations needed scales nearly linearly in the dimension of the space of spherical harmonics of degree at most $q$. The main idea is to exploit connections between spherical harmonics and zonal harmonics, and the fact that the zonal harmonics are the reproducing kernels of the space of degree $l$ spherical harmonics. Some numerical simulations are provided on synthetic examples to demonstrate the performance of the algorithm.", "questions": "1. In the optimization problem after line 42, isn’t the solution non-unique? If $g^*$ is a solution, then $g^* + g’$ for any g’ in the null space of the operator is also a solution? Also, isn’t f a solution of this optimization problem?\n\n2. Just to clarify my understanding, my first thought was to simply do linear regression in the basis of the lower degree spherical harmonics. But I suppose this is intractable since computing such a basis is computationally hard for even moderate values of q as noted in the paper. Is this correct? And so, this is why we can only hope for an approximate solution of Problem 1?\n\n3. In Definition 3, couldn’t we define the leverage function for any operator mapping $ L_2(\\mathbb{S}^{d-1})$ to itself? Or is it specifically defined for $\\mathcal{K}_d^{(q)}$?\n", "limitations": "I do not see the limitations discussed anywhere but this probably does not apply to this paper since it is essentially theoretical, and the sample complexity bounds are shown to be nearly optimal. Perhaps the running time which is super-quadratic in the number of samples could be improved upon?", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper is written well overall with a well-defined problem statement, and a clear description of related work.\n\n\n2. The results hold for any dimension $d$ which does not seem to have been handled previously. As noted in the related work, previous results typically applied to small, fixed values of $d$. In that respect, I think the results are quite strong.\n", "weaknesses": "1.  It’s a bit unclear to me whether the results hold only in the noiseless case, or is the method actually robust to noise. For instance, if we the function evaluations are corrupted with iid centred Gaussian noise, what can be said about the recovery error? The abstract mentions that the algorithm provides robust recovery, but I am not sure if this is what is proven.\n\n2.  While I understand the main contributions are theoretical, are there any real examples on which the method can be evaluated? At the moment, experiments are only conducted on synthetical examples.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable feedback on our paper. Below we answer your concerns and questions.\n\n> `It’s a bit unclear to me whether the results hold only in the noiseless case, or is the method actually robust to noise. For instance, if we the function evaluations are corrupted with iid centred Gaussian noise, what can be said about the recovery error? The abstract mentions that the algorithm provides robust recovery, but I am not sure if this is what is proven.`\n\nIn our paper, we consider a specific noise model to address the robustness of our method. We assume that the unknown function $f$, which we can query its values, is not necessarily a low-degree spherical harmonic and may contain high-degree components. These higher-degree components in the spherical harmonic expansion are treated as noise in our model. \n\nTo clarify, the noise we consider in our study is different from the typical iid noise that can corrupt measurements of the function. Instead, we assume that the noise values are drawn from an underlying and unknown $L_2(S^{d-1})$ function. Under this noise model, our algorithm can successfully recover the function up to a $(1+\\epsilon)$ factor of the noise's $2$-norm.\nOur results indeed hold under the above-specified noise model, which treats the higher degree components of $f$ as noise. Within this noise framework, we have proven the robustness of our algorithm.\n\n\n> `While I understand the main contributions are theoretical, are there any real examples on which the method can be evaluated?`\n\nOur method can be applied to any learning problem on the sphere. Spherical functions play a crucial role in problems with rotational-invariant property, where some real-world applications in machine learning is to recover 3D objective detection [1, 4, 5, 6], lighting estimation from images [2], and predicting atomic energies and forces [3].\n\n[1] Melnyk, Pavlo, et al. \"Steerable 3D spherical neurons.\" ICML 2022\n\n[2] Gardner, James, et al. \"Rotation-Equivariant Conditional Spherical Neural Fields for Learning a Natural Illumination Prior.\" NeurIPS 2022\n\n[3] Zitnick, Larry, et al. \"Spherical channels for modeling atomic interactions.\" NeurIPS 2022\n\n[4] Gardner, James, et al. \"Rotation-Equivariant Conditional Spherical Neural Fields for Learning a Natural Illumination Prior.\" NeurIPS 2022\n\n[5] Gerken, Jan, et al. \"Equivariance versus augmentation for spherical images.\" ICML 2022\n\n[6] Shakerinava, Mehran, and Siamak Ravanbakhsh. \"Equivariant networks for pixelized spheres.\" ICML 2021\n\n\n> `Just to clarify my understanding, my first thought was to simply do linear regression in the basis of the lower degree spherical harmonics. But I suppose this is intractable since computing such a basis is computationally hard for even moderate values of q as noted in the paper. Is this correct? And so, this is why we can only hope for an approximate solution of Problem 1?`\n\nYou’re right about solving a linear regression in the basis of the low-degree spherical harmonics. While this approach could solve the problem, it still requires discretizing the linear regression to deal with basis vectors which are continuous functions. \nOne of our main objectives is to minimize the number of samples needed to accurately recover the unknown function. However, solving the exact linear regression demands inner products between the basis vectors and the input function, leading to an infinite number of required samples from the input function.\n\nFurthermore, even though it's possible to construct a basis for spherical harmonics, as shown in Theorem 5.1 (https://arxiv.org/pdf/1304.2585.pdf), these basis functions are exceptionally complex. As a result, numerical calculations, even in moderate dimensions and with a moderate degree q, become quite challenging and intractable.\nOn the other hand, our approach is based on kernel regression, which ensures numerical stability even in high dimensions and for large degree q. This makes our method more practical and easier to implement compared to using spherical harmonics.\n\n> `Can the leverage function be defined for any operator?`\n\nYes, you are right the leverage function can be defined for any compact operator (please see Definition3 in https://arxiv.org/pdf/1812.08723.pdf ). We will modify our Definition 3 to make it clear that the leverage function is defined for any compact operator.\n"}, {"review_id": "JWzek8N7f2", "reviewer": "Reviewer_jKVP", "summary": "A technique is proposed to recover spherical harmonic expansions for functions defined on a d-dimensional sphere from a set of function evaluations.  Spherical harmonic expansions are recovered by solving an optimizaiton problem via a kernal approach, which is accompanied by theoretical guarantees.  Numerical experiments demonstrate phase transitions close to the theoretical bounds.", "questions": "How is the proposed method relevant for the field of deep learning on the sphere?  Can it offer a differentiable transform that can be integrated into approaches such as  [Cohen et al.](https://arxiv.org/abs/1801.10130)?\n\nWhat is meant by the comment: \"it is generally intractable to compute an orthogonal basis for the space of spherical harmonics, which renders the generalized Foureri series expansion in Lemma 2 primarily existential\"?  This comment is repeated twice in the manuscript. \n Could this please be elaborated?  ", "limitations": "No special negative societal impacts.", "rating": 4, "confidence": 5, "soundness": 4, "presentation": 2, "contribution": 2, "strengths": "A solid theoretical analysis is presented to derive a new approach to recovering spherical harmonic expansions from the evaluation of functions on the d-dimensional sphere.  It is shown that functions should be evalued uniformly randomly over the sphere.  The approach presented is accompanied by numerous theoretical results and guarantees.  Experimental results validate the theory presented, with phase transitions in success probabilities close to the theoretical bounds.", "weaknesses": "While the method is interesting, it does not seem that NeurIPS is the appropriate venue for this work.  While the field of deep learning on the sphere is an active area of research (e.g. [Cohen et al.](https://arxiv.org/abs/1801.10130)), this contribution would appear to be somewhat orthogonal to that body of literature.  No connection is made in the submitted manuscript beyond the final, somewhat cryptic, comment: \"We believe our finding would appeal to the readership of the community\".\n\nFurthermore, the connection to previous literature of the topic of fast spherical harmonic transforms and efficient sampling on the sphere is very poor.  No reference is made to the central works of [Driscoll & Healy](https://www.sciencedirect.com/science/article/pii/S0196885884710086) and [McEwen & Wiaux](https://ieeexplore.ieee.org/document/6006544) and the follow-up articles of their groups.\n\nThe use of the term \"near optimal\" in the title and throughout the article is somewhat overstated.  The proposed method is \"near\" optimal up to a logarithmic factor.  I would typically expect near optimal to be up to a constant factor. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable feedback on our paper. We greatly appreciate your insights and suggestions.\n\n> `It does not seem that NeurIPS is the appropriate venue for this work. While the field of deep learning on the sphere is an active area of research (e.g. Cohen et al.), this contribution would appear to be somewhat orthogonal to that body of literature.`\n\nWe appreciate your feedback on our work's relevance to deep learning on the sphere, inspired by research like Cohen et al. Although our work touches on some deep learning aspects, it's not exclusively focused on them. NeurIPS covers a broad range of research areas, with deep learning being just one facet among twelve distinct research domains highlighted in the NeurIPS 2023 Call for Papers. Our work aligns more closely with \"Probabilistic methods,\" a specific area specified in the Call for Papers. We believe our research introduces a fresh perspective and valuable insights to the NeurIPS community. While we recognize potential concerns about its interdisciplinary nature, we firmly believe that our contributions can foster novel ideas and collaborations across research threads, ultimately enhancing the scientific landscape.\n\n\n> `No reference is made to the central works of Driscoll & Healy and McEwen & Wiaux and the follow-up articles of their groups.`\n\nRest assured, we fully recognize the significance of these foundational works and in our final submission version, we will provide a comprehensive review of the key works by Driscoll & Healy and McEwen & Wiaux. Additionally, we will include references to the relevant follow-up articles to ensure we properly acknowledge their contributions.\n\n\n> `The use of the term \"near optimal\" in the title and throughout the article is somewhat overstated.`\n\nWe want to highlight the common practice of describing \"near optimal\" performance up to a logarithmic factor, which is widely recognized in algorithmic literature.\nTo provide some context, many prominent works in the field have adopted this perspective. For instance, the paper by Lin, Tianyi, Chi Jin, and Michael I. Jordan titled \"Near-optimal algorithms for minimax optimization\" (COLT 2020) employs the $\\tilde{O}$ notation to discuss near-optimal solutions. Similarly, the study by Candès and Terence Tao titled \"The power of convex relaxation: Near-optimal matrix completion\" (IEEE Transactions on Information Theory 2010) employs this notation to describe the performance of their proposed approach.\nWe believe that the use of \"near optimal\" to indicate optimality up to a logarithmic factor provides a more nuanced and accurate representation of our method's performance and aligns with the broader trend in algorithmic analysis.\n\n\n"}, {"review_id": "hSmvpo2SKS", "reviewer": "Reviewer_L1Mq", "summary": "### Result ###\n\nThe paper studies the problem of recovering a function from a finite number of noisy observations, for the class of \"square-integrable functions on the unit sphere\" (denoted by $L^2(\\mathbb{S}^{d-1})$ when the sphere in question is of degree $d$).\n\nThe paper shows, by developing an algorithm, that the number of samples required for recovery up to an $\\epsilon$-multiplicative error is proportional to $\\beta_{q, d}$, where $q$ is the degree of spherical harmonics desired and $\\beta_{q, d}$ is the dimension of the space of degree $q$ spherical harmonics on the sphere of dimension $d-1$. This sample complexity is optimal, as shown by the paper via lower bounds. \n\n---------------\n\n### Broader Contributions ###\n\nTo achieve the result above, the paper makes the following conceptual and technical contributions. \n\nFirst, the problem at hand is a regression problem with potentially infinite-dimensional continuous cost functions. Inspired by the approach of Avron, Kapralov, Musco, Musco, Velingker, and Zandieh (AKMMVZ19), which discretizes this problem via the leverage scores of the regression matrix. \n\nSecond, the authors utilize connections between spherical harmonics and zonal harmonics to enable the implementation of leverage score sampling of the regression matrix in question. This leverage score sampling result is novel. \n\n", "questions": "### Questions ### \n\n1. Would the authors be able to situate their results in the context of machine learning? For instance, in lines 20-21, the potential applications mentioned seem to be leaning towards physics/astronomy; it would be quite interesting to see such applications in machine learning (e.g., published in prior ICML/NeurIPS conferences) and also some that might be somewhat more recent. \n\n2. I'm actually slightly confused about the main result (Theorem 1 in the full version of the submission; lines 53-56 in the main_full.pdf) for the following reason: In \"Randomized algorithms for matrices and data\" by Mahoney (2011), we see that the (essentially optimal) sample complexity of least-squares regression is $d/\\epsilon^2$, obtained by leverage score sampling. Here, $d$ can be thought of as the \"intrinsic dimension\" of the projection matrix in question; I'm surprised the result in the paper (Theorem 1 in the full submission) seems to have a better dependence of $\\epsilon^{-1}$ on the sample complexity despite (seemingly) being a more general result. Would the authors be able to comment on this? \n\n--------------------------\n\n### Suggestions ### \n\n1.  I think it would really help with readability if the authors could present \"specific cases\" of their statements wherever possible. As an example, the Definition 3, Equation 8 is essentially the definition of the sensitivities function (introduced by Langberg and Schulman; for a recent paper with this fact explicitly stated, see, for example, \"Sharper Bounds for $\\ell_p$ Sensitivity Sampling\" by Woodruff and Yasuda); similarly, the minimum characterization of the leverage score function is essentially the one seen in Lemma 2 of \"Uniform Sampling for Matrix Approximation\" by Cohen, Lee, Musco, Musco, Peng, and Sidford). It would make it easier for readers to see these (possibly more familiar) statements first followed by the generalized versions stated in the paper. \n\n\n", "limitations": "N/A ", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "### Result's Strength ### \n\nI think the optimality of the stated result of the paper makes the paper mathematically strong.  \n\n----------------------------------------\n\n# During the Rebuttal Phase # \n\nThe author's explanations, repeated readings of the paper, and Reviewer BcP8's review and questions have helped me understand and appreciate the paper better. I'm therefore raising the score from 4 to 6 and confidence from 2 to 3. \n", "weaknesses": "### Difficulty in reading. ### \n\nI found the paper quite difficult to read though I could tell it was written well. I think this is simply because of the deeply technical nature of the problem and my unfamiliarity with the topic (at least to the level of generality of this paper) within the context of machine learning. \n\nI am not sure if I have any concrete feedback to this end but due to this particular reason I feel this paper might be a much better fit (in terms of reaching a wide audience who'd actually understand and appreciate the results) at a mathematical/optimization journal. That said, I acknowledge, based on my unfamiliarity with the topic, that my judgement could be completely misplaced. \n\nI'll be grateful to the authors if they could situate their work more in the context of machine learning. I'll also be happy to keep reading the submission and improving my understanding through the rebuttal period.\n\n----------------------\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We greatly appreciate the reviewer's detailed and insightful comments. \n\n> `Would the authors be able to situate their results in the context of machine learning? For instance, in lines 20-21, the potential applications mentioned seem to be leaning towards physics/astronomy; it would be quite interesting to see such applications in machine learning (e.g., published in prior ICML/NeurIPS conferences) and also some that might be somewhat more recent.`\n\nOur main focus is recovering (unknown) functions defined in the sphere, a critical task in scenarios where rotational invariance is a fundamental property. In real-life machine learning applications, this property becomes very important as a foundational requirement for modeling 3D point clouds. Notable examples that appeared in top-tier machine learning conferences include molecular/atom systems, where understanding the underlying functions within a spherical context can significantly enhance predictive modeling and simulation accuracy [1,2,3,4]. Other examples are in the field of computer vision, specifically in the recognition, classification, and reconstruction of 3D objects [5,6,7,8].\n\n[1] Eickenberg, Michael, et al. \"Solid harmonic wavelet scattering: Predicting quantum molecular energy from invariant descriptors of 3D electronic densities.\" NeurIPS 2017\n\n[2] Frank, Thorben, et al. \"So3krates: Equivariant attention for interactions on arbitrary length-scales in molecular systems.\" NeurIPS 2022\n\n[3] Zitnick, Larry, et al. \"Spherical channels for modeling atomic interactions.\" NeurIPS 2022\n\n[4] Liu, Yi, et al. \"Spherical message passing for 3d molecular graphs.\" ICLR 2022\n\n[5] Gardner, James, et al. \"Rotation-Equivariant Conditional Spherical Neural Fields for Learning a Natural Illumination Prior.\" NeurIPS 2022\n\n[6] Melnyk, Pavlo, et al. \"Steerable 3D spherical neurons.\" ICML 2022\n\n[7] Gerken, Jan, et al. \"Equivariance versus augmentation for spherical images.\" ICML 2022\n\n[8] Shakerinava, Mehran, and Siamak Ravanbakhsh. \"Equivariant networks for pixelized spheres.\" ICML 2021\n\n\n> `I'm actually slightly confused about the main result (Theorem 1 in the full version of the submission; lines 53-56 in the main_full.pdf) for the following reason: In \"Randomized algorithms for matrices and data\" by Mahoney (2011), we see that the (essentially optimal) sample complexity of least-squares regression is d/eps^2, obtained by leverage score sampling. Here, d can be thought of as the \"intrinsic dimension\" of the projection matrix in question; I'm surprised the result in the paper (Theorem 1 in the full submission) seems to have a better dependence of eps^{-1} on the sample complexity despite (seemingly) being a more general result. Would the authors be able to comment on this?`\n \nIn our proof of Theorem 1, we invoked Theorem 6.3 from the paper [Chen, Price'19]. Your intuition is correct that in order to achieve a subspace embedding guarantee for the design matrix in the regression problem with an approximation factor $\\epsilon$ (i.e. to preserve all singular values of the gram matrix to within a factor of $1 \\pm \\epsilon$), it is generally necessary to have about $d \\log d / \\epsilon^2$ leverage score samples. The quadratic dependence on $1/\\epsilon$ arises from the birthday's paradox, and the $\\log d$ factor comes from the coupon-collector problem. The approximate regression results in Mahoney (2011) are derived from subspace embedding guarantee with error factor $\\epsilon$.\n\nHowever, when approximately solving the regression problem, it is possible use fewer samples since we do not need to guarantee subspace embedding with an error parameter $\\epsilon$. Instead, it suffices to have a subspace embedding that preserves the eigenvalues of the design matrix up to some constant factor, like $1 \\pm 1/2$. For more precise details, please refer to the definition of \"well-balanced sample\" in Definition 2.1 of [Chen, Price'19] which is a sufficient condition for approximately solving linear regression.\n\n[Chen, Price'19]: \"Active regression via linear-sample sparsification.\" COLT 2019.\n\n\n> `I think it would really help with readability if the authors could present \"specific cases\" of their statements wherever possible. ... It would make it easier for readers to see these (possibly more familiar) statements first followed by the generalized versions stated in the paper.`\n\nYour suggestion regarding enhancing readability is highly appreciated. To make it easier for readers to understand the concepts, we'll include references that they're familiar with. We'll follow your advice by adding specific cases to all our definitions. We agree with you that this will help readers start with recognizable statements before moving on to the generalized versions in our paper.\n\n"}]

Author rebuttal:

Meta-review: This paper uses randomized linear algebra (specifically leverage scores sampling) to solve near-optimally the natural problem of fitting a function on the sphere with spherical harmonics using a limited number of samples. Most reviewers felt the work studies and interesting problems, contains a strong theoretical result (not the most novel, but novel enough) and is reasonably well presented, We hope the authors take into account suggestions for how to improve clarity and accessibility to the broader NeurIPS audience when preparing their final version of the paper.
Decision: Accept (poster)