The steerable graph Laplacian
and its application to filtering image data-sets

arXiv:1802.01894v1 [] 6 Feb 2018

Boris Landa and Yoel Shkolnisky
February 7, 2018

Boris Landa
Department of Applied Mathematics, School of Mathematical Sciences
Tel-Aviv University
sboris20@gmail.com

Yoel Shkolnisky
Department of Applied Mathematics, School of Mathematical Sciences
Tel-Aviv University
yoelsh@post.tau.ac.il

Please address manuscript correspondence to Boris Landa, sboris20@gmail.com, (972) 549427603.

1

Abstract
In recent years, improvements in various scientific image acquisition techniques gave rise to the
need for adaptive processing methods aimed for large data-sets corrupted by noise and deformations.
In this work, we consider data-sets of images sampled from an underlying low-dimensional manifold
(i.e. an image-valued manifold), where the images are obtained through arbitrary planar rotations.
To derive the mathematical framework for processing such data-sets, we introduce a graph Laplacianlike operator, termed steerable graph Laplacian (sGL), which extends the standard graph Laplacian
(GL) by accounting for all (infinitely-many) planar rotations of all images. As it turns out, a properly
normalized sGL converges to the Laplace-Beltrami operator on the low-dimensional manifold, with
an improved convergence rate compared to the GL. Moreover, the sGL admits eigenfunctions of the
form of Fourier modes multiplied by eigenvectors of certain matrices. For image data-sets corrupted by
noise, we employ a subset of these eigenfunctions to “filter” the data-set, essentially using all images
and their rotations simultaneously. We demonstrate our filtering framework by de-noising simulated
single-particle cryo-EM image data-sets.

1

Introduction

Developing efficient and accurate processing methods for scientific image data-sets is a central research
task, which poses many theoretical and computational challenges. In this work, motivated by certain
experimental imaging and tomography problems [10, 32], we put our focus on the task of reducing the
noise in a large data-set of images, where the in-plane rotation of each image is arbitrary.
To accomplish such a task, it is generally required to have prior knowledge, or a model assumption, on
the data-set at hand. One popular approach is to assume that the data lies on a low-dimensional linear
subspace, whose parameters can then be estimated by the ubiquitous Principal Components Analysis
(PCA). In our setting, where images admit arbitrary planar rotations, it is reasonable to incorporate
all rotations of all images into the PCA procedure, resulting in what is known as “steerable PCA”
(sPCA) [18, 34, 35].
In practice, however, experimental data-sets typically admit more complicated non-linear structures.
Therefore, we adopt the more flexible notion that the images were sampled from a low-dimensional
manifold M embedded in a high-dimensional Euclidean space, an assumption that lies at the heart

of many effective machine learning, dimensionality reduction, and signal processing techniques (see for
example [25, 31, 4, 7, 23]).
When processing and analyzing manifold data, a fundamental object of interest is the Laplace-Beltrami
operator ∆M [24], which encodes the geometry and topology of M. Essentially, the Laplace-Beltrami

operator is a second-order differential operator generalizing the classical Laplacian, and can be therefore
considered as accounting for the smoothness of functions on M. In this context, it is a common approach

to leverage the Laplace-Beltrami operator and its discrete counterpart, the graph Laplacian [4], to process
surfaces, images, and general manifold data [20, 30, 33, 19, 9, 16, 21]. Incorporating the Laplacian in data
processing algorithms typically follows one of two approaches. The first is based on solving an inverse
problem which includes a regularization term involving the Laplacian of the estimated data coordinates,
and the second is based on directly using the Laplacian or its eigenfunctions for filtering the data-set.
We mention that the eigenfunctions of the Laplace-Beltrami operator, which we refer to as “manifold
harmonics”, are analogous to classical Fourier modes as they constitute a basis on M favorable for
2

expanding smooth functions. Here, we focus on the the second approach, namely on filtering the dataset by the manifold harmonics. In particular, we assume (as mentioned above) that each point in the
data-set is a high-dimensional point that lies on some manifold M with low intrinsic dimension. Thus,
the coordinate functions of M (which are functions defined on the manifold) can be expanded by the

manifold harmonics of M and filtered by truncating the expansion (see for example [20, 33] for similar
approaches in the contexts of image processing and surface fairing).

As the manifold M is unknown a priori, we do not have access to its Laplace-Beltrami operator
directly. Consequently, it must approximated from the data, which can be achieved through the graph
Laplacian [4, 7]. Specifically, given points {x1 , . . . , xN } ⊂ RD , we consider the fully connected graph
Laplacian, denoted as L ∈ RN ×N and given by
L = D − W,

Wi,j

n

o
= exp − kxi − xj k /ε ,
2

Di,i =

N
X

Wi,j ,

(1)

j=1

where W is known as the affinity matrix (using the Gaussian kernel parametrized by ε), and D is a
diagonal matrix with {Di,i } on its diagonal. Then, as was shown in [7, 26, 4, 5], the normalized graph
Laplacian L̃ = D −1 L converges to the negative-defined Laplace-Beltrami operator ∆M when ε → 0 and
N → ∞. In particular, it was shown in [26] that for a smooth function f : M → R
N

4X
1
L̃i,j f (xj ) = ∆M f (xi ) + O( 1/2 1/2+d/4 ) + O(ε),
ε
N ε

(2)

j=1

where d is the intrinsic dimension of M. Therefore, it is evident that for a fixed parameter ε, the

error in the approximation of ∆M depends directly on the intrinsic dimension d and inversely on the
number of data points N . If d is large, then we need a large number of samples to achieve high accuracy.
In our scenario, as images admit arbitrary planar rotations, the number of images required to use the
approximation (2) may be prohibitively large as images which differ only by an in-plane rotation may not
be encoded as similar by the affinity matrix W (since the Euclidean distance between them may be large).
To overcome this obstacle, we construct the steerable graph Laplacian, which is conceptually similar to
the standard graph Laplacian, except that it also accounts for all rotations of the images in the data-set.
We then propose to employ the eigenfunctions of this operator to filter our image data-set, allowing for
an efficient procedure for mitigating noise.
The contributions of this paper are as follows. First, we introduce and analyze the steerable graph
Laplacian operator, characterize its eigen-decomposition (together with a general family of operators),
and show that it can be diagonalized by Fourier modes multiplied by eigenvectors of certain matrices.
Second, we introduce the normalized steerable graph Laplacian, and demonstrate that it is more accurate
than the standard graph Laplacian in approximating the Laplace-Beltrami operator, in the sense that it
admits a smaller variance error term. Essentially, the improved variance error term can be obtained by
replacing d in equation (2) with d− 1. Third, we propose to employ the eigenfunctions of the (normalized)
steerable graph Laplacian for filtering image data-sets, where the explicit appearance of Fourier modes
in the form of the eigenfunctions allows for a particularly efficient filtering procedure. To motivate and
justify our approach, we provide a bound on the error incurred by approximating an embedded manifold

3

by a truncated expansion of its manifold harmonics. We also analyze our approach in the presence of
white Gaussian noise, and argue that in a certain sense our method is robust to the noise, and moreover,
allows us to reduce the amount of noise inversely to the number of images in the data-set.
The paper is organized as follows. Section 2.1 lays down the setting and provides the basic notation and
assumptions. Then, Section 2.2 defines the steerable graph Laplacian and derives some of its properties,
including its eigen-decomposition. Section 2.3 presents the normalized steerable graph Laplacian and
derives its convergence rate to the Laplace-Beltrami operator, while providing its eigen-decomposition
similarly to the preceding section. Section 2.4 numerically corroborates the convergence rate of the
normalized steerable graph Laplacian by a simple toy example, and Section 2.5 proposes and analyzes
a filtering scheme for image data-sets based on the eigenfunctions of the (normalized) steerable graph
Laplacian. Section 3 summarizes all relevant algorithms and presents the computational complexities
involved. Section 4 provides an analysis of our approach in the presence of white Gaussian noise, followed
by Section 5 which demonstrates our method for de-noising a simulated cryo-EM image data-set. Lastly,
Section 6 provides some concluding remarks and possible future research directions.

2
2.1

Setting and main results
The setting

Suppose that we have N points {x1 , . . . , xN } ⊂ CD sampled from a probability distribution p(x), which
is restricted to a smooth and compact d-dimensional submanifold M without boundary. Furthermore,

each point x ∈ M is associated with an image through a correspondence between points in the ambient

space CD and images. Specifically, each point x ∈ CD corresponds to an image I(r, θ) ∈ L2 (D), where D
is the unit disk, by

I(r, θ) =

ℓm
M
X
X

xm,ℓ ψm,ℓ (r, θ),

ψm,ℓ (r, θ) = Rm,ℓ (r)eımθ ,

(3)

m=−M ℓ=1

where xm,ℓ is the (m, ℓ)’th coordinate of x, and {ψm,ℓ } is an orthogonal basis of L2 (D) whose radial part
is {Rm,ℓ }ℓ (orthogonal on [0, 1) w.r.t the measure rdr). In other words, the points xi sampled from the

manifold M are the expansion coefficients of some underlying images in the basis {ψm,ℓ }. We mention

that the points xi do not correspond to the pixels of the images directly since such a representation does
not allow for a natural incorporation of planar rotations. We shall refer to m ∈ Z as the angular index,
and to ℓ ∈ N+ as the radial index, where {ℓm } of (3) are the numbers of radial indices taking part in the
P
expansion for each angular index m, satisfying M
m=−M ℓm = D. Therefore, the data-set can be organized
as the N × D matrix

z

 x1,(−M,1)

 x
X = 
2,(−M,1)

..


.


xN,(−M,1)

m=−M

}|
{

. . . x1,(−M,ℓ−M )

. . . x2,(−M,ℓ−M ) 
 ...
..

.

. . . xN,(−M,ℓ−M )
4


}|
{

. . . x1,(M,ℓM ) 

. . . x2,(M,ℓM ) 
,

..

.


. . . xN,(M,ℓM )

m=M

z


x1,(M,1)

 x2,(M,1)

..

.


xN,(M,1)

(4)

where xi,(m,ℓ) denotes the (m, ℓ)’th coordinate (with angular frequency m and radial frequency ℓ) of the
i’th data-point xi .
Representing image data-sets via their expansion coefficients obviously does not impose any restrictions, as any image data-set can be first expanded in basis functions of the form of ψm,ℓ (see Remark 1
below), and all subsequent analysis can be carried out in the domain of the resulting expansion coefficients

xi,(m,ℓ) .

Basis functions of the form of ψm,ℓ , which are separable in polar coordinates into radial functions
Rm,ℓ (r) multiplied by Fourier modes eımθ , are called “steerable” [14, 22], as they allow for simple and
efficient rotations. In particular, every ψm,ℓ can be rotated by multiplying it with a complex constant
ψm,ℓ (r, θ + ϕ) = eımϕ ψm,ℓ (r, θ)

(5)

and thus, we can describe image rotation by modulation of the expansion coefficients, with each coefficient xm,ℓ transformed into xm,ℓ eımϕ . Consequently, we endow the ambient space CD with the rotation
operation R : CD × [0, 2π) → CD , defined as
ımϕ
xϕ
.
m,ℓ = xm,ℓ e

R(x, ϕ) , xϕ ,

(6)

Therefore, if x is the coefficients vector of the image I, then xϕ is the coefficients vector of the image I
rotated by ϕ, obtained by modulating each coefficient appropriately.
Lastly, we assume that the manifold M is rotationally-invariant, that is, it is closed under R, such

that for every x ∈ M and ϕ ∈ [0, 2π) we have that R(x, ϕ) ∈ M. A key observation here, is that this
property enables us to generate new data-points on M by rotating existing images.

Our goal is to derive adaptive processing methods for our image data-set, allowing for filtering and

de-noising, while making use of the rotation-invariance of M to provide accurate and efficient algorithms.
Remark 1. Examples for bases of the form of ψm,ℓ include the 2D Prolate Spheroidal Wave Functions
(PSWFs) [29, 17], the Fourier-Bessel functions [35], and data-adaptive steerable principal components [18,
34], all of which allow approximating image data-sets provided by their samples. We note that the choice
of the particular basis may depend on the application and specific model assumptions.
Remark 2. It is important to mention that our framework can also support functions/images sampled on
a polar grid, and as a special case – 1D periodic signals. That is, in place of eq. (3), every point x ∈ CD

can be defined via the correspondence

I(rℓ , θ) =

M
X

xm,ℓ eımθ ,

(7)

m=−M

where ℓ = 1, . . . , ℓ̂ enumerates over the different radii of I. In the case that ℓ̂ = 1, each point x corresponds
exactly to a 1D periodic signal. Then, if images/functions sampled on a polar grid are provided, xm,ℓ can
be computed efficiently by the FFT of the (equally-spaced) angular samples of I for each radius.

5

2.2

The steerable graph Laplacian for image-manifolds

To derive a natural basis on the manifold M, we employ graph Laplacian operators which encode the

geometry and topology of M. To this end, since the manifold M is rotationally-invariant, we propose to
form a graph Laplacian over the points {x1 , . . . , xN } and all of their (infinitely many) rotations.

We start by defining an appropriate function space for constructing our operators. Consider the

domain Γ = {1, ..., N } × S1 , where S1 is unit circle (parametrized by an angle ϑ ∈ [0, 2π)), and functions

2 1
f : Γ → C of the form f (i, ϑ) = fi (ϑ), with {fi }N
i=1 ∈ L (S ). The space of the functions f is defined as

H = L2 (Γ), which is a Hilbert space endowed with the inner product
hg, f iH =

N Z
X

2π

0

i=1

gi∗ (ϑ)fi (ϑ)dϑ

(8)

for any f, g ∈ H, where (·)∗ denotes complex-conjugation. Loosely speaking, every f ∈ H can be
considered as a column vector of periodic functions, namely f = [f1 (ϑ), . . . , fN (ϑ)]T , assigning a value to
every index i ∈ {1, ..., N } and an angle ϑ.

In order to capture pair-wise similarities between different points and rotations in our data-set, we

define the steerable affinity operator W : H → H as
{W f } (i, ϑ) =

N Z
X
j=1


Wi,j (ϑ, ϕ) = exp − xϑi − xϕ
j

2π

Wi,j (ϑ, ϕ)fj (ϕ)dϕ,

0

2


/ε ,

(9)

where f ∈ H, 1 ≤ i, j ≤ N , {ϑ, ϕ} ∈ [0, 2π), ε is a tunable parameter, and xϑi stands for the rotation of
xi by an angle ϑ (via (6)). Therefore, W can be considered as describing the affinity between any two

rotations of any two points in our data-set. Note that since {ψm,ℓ } (of (3)) are orthonormal, the distance
xϑi − xϕ
j

2

agrees with the natural distance (in L2 (D)) between the images corresponding to xi and xj ,

after rotating them by ϑ and ϕ, respectively.

Before we proceed to the define the steerable graph Laplacian over H, we mention that we lift any

complex-valued matrix A ∈ CN ×N to act over H by

{Af } (i, ϑ) =

N
X

Ai,j fj (ϑ),

(10)

j=1

for any f ∈ H. Then, we define the (un-normalized) steerable graph Laplacian L : H → H by
Lf = Df − W f,

Di,i =

N Z
X
j=1

2π

Wi,j (0, α)dα,

(11)

0

where D is a diagonal matrix with {Di,i }N
i=1 on its diagonal. If we implicitly augment our data-set to
include all planar rotations of all images, then the steerable graph Laplacian can be viewed as the standard
graph Laplacian (equation (1)) constructed from the (infinitely-many) data points of the augmented dataset.

6

Similarly to the standard graph Laplacian, we show in Appendix B that L admits the quadratic form
Z
N Z
1 X 2π 2π
hf, Lf iH =
Wi,j (ϑ, ϕ) |fi (ϑ) − fj (ϕ)|2 dϑdϕ,
2
0
0

(12)

i,j=1

which is analogous to the quadratic form of the standard graph Laplacian (see [4]) in the sense that it
accounts for the regularity of the function f over the domain Γ w.r.t the pair-wise similarities Wi,j (ϑ, ϕ)
(measured between different data-points and rotations). In other words, the quantity hf, Lf iH penalizes

large differences |fi (ϑ) − fj (ϕ)| particularly when Wi,j (ϑ, ϕ) is large, i.e. when the images corresponding
to xi and xj , rotated by ϑ and ϕ, respectively, are similar. Therefore, hf, Lf iH is expected to be small
for functions f which are smooth (in a certain sense) over Γ with the geometry induced by W .

As we expect the operator L to encode certain geometrical aspects of our data-set, as in the case of
the standard graph Laplacian (see [4, 7]), it is beneficial to investigate its eigen-decomposition. In this
context, it is important to mention that a naive evaluation of W (and consequently L) by discretizing all
rotation angles is generally computationally prohibitive, and moreover, is less accurate then considering
the continuum of all rotation angles. To obtain the eigen-decomposition of L, we demonstrate in Appendix C that the steerable graph Laplacian L belongs to a family of operators, which we term Linear and
Rotationally-Invariant (LRI), that admit eigenfunctions with a convenient analytic form. In particular,
we show that LRI operators can be diagonalized by tensor products between Fourier modes and vectors
in CN , where the vectors can be computed efficiently by diagonalizing a certain sequence of matrices. In
the case of the steerable graph Laplacian L, this inherently stems from the fact that Wi,j (ϑ, ϕ) is only a
function of ϕ − ϑ (following immediately from (6) and (9)), and therefore can be expanded in a Fourier
series as

∞
1 X
(m)
Wi,j (ϑ, ϕ) = Wi,j (0, ϕ − ϑ) =
Ŵ e−ım(ϕ−ϑ) ,
2π m=−∞ i,j

(m)
Ŵi,j

=

Z

2π

Wi,j (0, α)eımα dα,

(13)

0
(m)

where 1 ≤ i, j ≤ N . We define the matrix Ŵ (m) whose (i, j)’th entry is Ŵi,j , and observe from (13)
n
o∞
provides a complete characterization of the steerable affinity
that the sequence of matrices Ŵ (m)
m=−∞
n
o∞
also plays a key
operator W (and consequently L). Therefore, the sequence of matrices Ŵ (m)
m=−∞

role in the evaluation of the eigen-decomposition of L, as detailed by the following theorem.

Theorem 1. The steerable graph Laplacian L admits a sequence of non-negative eigenvalues {λm,1 , . . . , λm,N }∞
m=−∞ ,
and a sequence of eigenfunctions {Φm,1 , . . . , Φm,N }∞
m=−∞ which are orthogonal and complete in H and

are given by

Φm,k = vm,k · eımϑ ,

(14)

where vm,k and λm,k are the k’th eigenvector and eigenvalue, respectively, of the matrix
Sm = D − Ŵ (m) ,
and D and Ŵ (m) are given by (11) and (13), respectively.
The the proof is provided in Appendix D.
7

(15)

We point out that Di,i =
the matrices

Ŵ (m) ,

PN

(0)
j=1 Ŵi,j ,

hence all quantities involving Sm can be computed directly from

which in turn can be approximated (to arbitrary precision) by
(m)

Ŵi,j ≈

K−1
2π X
Wi,j (0, 2πk/K)eı2πmk/K
K

(16)

k=0

for a sufficiently large integer K, and evaluated rapidly using the FFT.
Analogously to the separation of variables of the basis functions ψm,ℓ of (3), the basis functions Φm,k
of (14) adopt a similar separation into products of vectors vm,k ∈ CN and Fourier modes eımϑ . As such,

we consider Φm,k as “steerable” over H, and hence the term steerable in “steerable graph Laplacian”.

2.3

Normalized steerable graph Laplacian and the Laplace-Beltrami operator

In the previous section, we constructed and analyzed the steerable graph Laplacian L, which can be
considered as a generalization of the standard graph Laplacian. In particular, the steerable graph Laplacian inherits many of the favorable properties of the graph Laplacian. Based upon the construction in
Section 2.2, in what follows we consider a certain normalized variant of L which not only provides us with
steerable basis functions adapted to our data-set, but moreover, is shown to approximate the continuous
(negative-defined) Laplace Beltrami operator ∆M .
We start by defining the normalized steerable graph Laplacian L̃ : H → H, similarly to the normalized

variant of the standard graph Laplacian (see [7]), as

L̃ = D −1 L,

(17)

where D −1 is the inverse of the matrix D from (11). Explicitly, we have that L̃f = f − D −1 W f for every
f ∈ H. It then turns out that the normalized steerable graph Laplacian L̃ converges to the negativedefined Laplace-Beltrami operator ∆M [24] when ε → 0 and N → ∞, while improving on the convergence
rate of the standard (normalized) graph Laplacian (equation (2)), as reported by the next theorem.
Theorem 2. Suppose that

P

m6=0

Pℓ m

2
ℓ=1 |xm,ℓ |

> 0 for all x ∈ M (up to a set of measure zero), and

let {x1 , . . . , xN } ∈ M be i.i.d with probability distribution p(x) = 1/V ol {M}, i.e. uniform sampling

distribution. If f : M → R is a smooth function, and we define g ∈ H s.t. g(i, ϑ) = f (xϑi ) (where xϑi is
given by (6)), then with high probability we have that

4n o
1
L̃g (i, ϑ) = ∆M f (xϑi ) + O( 1/2 1/2+(d−1)/4 ) + O(ε).
ε
N ε

(18)

The proof is provided in Appendix E. Comparing (18) with (2), it is evident that both graph Laplacians
converge to ∆M with the same bias error term of O(ε). However, the steerable graph Laplacian admits
a smaller variance error term (second term from the right in (18)), which depends on d − 1 instead of

d. The intuition behind this improved bound is that the steerable graph Laplacian takes all rotations of
all images into consideration, and so it analytically accounts for one of the intrinsic dimensions of M,

that is, the dimension corresponding to the rotation R (see (6)). A numerical example demonstrating the
improved convergence rate due to Theorem 2 can be found in Section 2.4.

8

Remark 3. The condition

P

m6=0

Pℓ m

2
ℓ=1 |xm,ℓ |

> 0 in Theorem 2 essentially requires that the images

associated with the points of M are not radially-symmetric (i.e. have a non-constant angular part). This
is because the coordinates xm,ℓ of x corresponding to the angular index m = 0 contribute only to the

radial part of the image (see equation (3)). Of-course, if the images are all radially-symmetric, then the
steerable graph Laplacian would not provide any improvement over the convergence rate of the standard
graph Laplacian.
Next, we evaluate the eigenfunctions and eigenvalues of the normalized steerable graph Laplacian L̃,
where analogously to Theorem 1, the next theorem relates the eigenfunctions and eigenvalues of L̃ to the
matrices Ŵ (m) of (13).
Theorem
3. The
a sequenceoof non-negative eigenvalues
n
o∞normalized steerable graph Laplacian L̃ admits
n
∞
, and a sequence of eigenfunctions Φ̃m,1 , . . . , Φ̃m,N
which are complete
λ̃m,1 , . . . , λ̃m,N
m=−∞

in H and are given by

m=−∞

Φ̃m,k = ṽm,k · eımϑ ,

(19)

where ṽm,k and λ̃m,k are the k’th eigenvector and eigenvalue, respectively, of the matrix
S̃m = I − D −1 Ŵ (m) ,

(20)

I is the N × N identity matrix, and D and Ŵ (m) are given by (11) and (13), respectively.
The proof is provided innAppendix
o F.
Let us denote the basis Φ̃m,k of (19) by Φ̃. Due to the convergence of L̃ to the Laplace-Beltrami

operator ∆M , we consider Φ̃ as a basis adapted to our data-set through the n
geometry
o and topology of
M, and hence a favorable basis for expanding and filtering our data-set. Since Φ̃m,k are also steerable,

we shall refer to them (with a slight abuse of notation) as steerable manifold harmonics.

2.4

Toy example

At this point, we wish to demonstrate our setting as well as the improved convergence rate of the steerable
graph Laplacian by the following example. Consider images of the form
I(r, θ) = x0,1 R0,1 (r) + x1,1 R1,1 (r)eıθ ,

(21)

which is a special case of (3), where R0,1 , R1,1 are arbitrary radial functions, and M = 1, ℓ−1 = 0, ℓ0 =
1, ℓ1 = 1. Additionally, we take the unit sphere S2 (d = 2) in R3 , and embed it in C2 by mapping every
point p = [px , py , pz ] ∈ S2 (px , py , pz are the x, y, z coordinates) to the point x = [x0,1 , x1,1 ] ∈ M via
x0,1 = pz ,

x1,1 = px + ıpy .

Note that the rotation operation R of (6) in this case is
R(x, ϕ) =

"

1

0

0 eıϕ
9

#
#"
x0,1
x1,1

,

(22)

(23)

which is equivalent to rotating the point p ∈ R3 (corresponding to x) in the xy-plane as
 
px
 

− sin(ϕ) cos(ϕ) 0 py  .
0
0
1
pz


cos(ϕ)

sin(ϕ) 0

(24)

Hence, all rotations of all images sampled from the sphere remain on the sphere, and therefore M is
rotationally-invariant (as defined in Section 2.1).

In order to demonstrate numerically the convergence rate of the (normalized) steerable graph Laplacian
to the Laplace-Beltrami operator (as asserted by Theorem 2), we chose a test function f : M → R
f (x) = Re {x1,1 } + x0,1 ,

(25)

and a testing point x0 = [0, 1] (corresponding to p = [1, 0, 0] on S2 ), for which ∆M f (x0 ) = −2 (see

example in [26]). We then uniformly sampled N = 2, 000 points {x1 , . . . , xN } from M and approximated

∆M f by applying the steerable graph Laplacian L̃. Specifically, ∆M f (x0 ) was approximated from (18)
and (17) by defining g(i, ϑ) = f (xϑi ) for i = 0, 1, . . . , N and computing




"
#
PN R 2π
ϕ
)dϕ
W
(0,
ϕ)f
(x
4
0,j
4
j=0
j
0
−1

D0,j
W0,j (0, ϕ)f (xϕ
f (x0 ) −
L̃g (0, 0) =
f (x0 ) −
PN R 2π
j )dϕ =
ε
ε
ε
0
j=0 0 W0,j (0, ϕ)dϕ
j=0
"
PN PK−1
2πk/K #
)
4
j=0
k=0 W0,j (0, 2πk/K)f (xj
,
(26)
f (x0 ) −
≈
PN PK−1
ε
j=0
k=0 W0,j (0, 2πk/K)

4n

o

N Z
X

2π

where W is given by (9), D is given in (11), and we replaced integration with summation using a sufficiently
large integer K. Note that f (x0 ) = 1, and by (6) we have that
2πk/K

f (xj

o
n
) = Re xj,(1,1) eı2πk/K + xj,(0,1) ,

(27)

where xj,(m,ℓ) is the (m, ℓ)’th coordinate of the j’th point. Figure 1 depicts the errors of estimating
∆M f (x0 ) using the steerable graph Laplacian (equation (26)) versus the standard graph Laplacian (equations (1) and (2)), for K = 256 and different values of ε. The slope of the log-error in the variancedominated region (obtained by a linear curve fit and averaged over 1, 000 experiments) is −0.97 for the

standard graph Laplacian, and −0.74 for the steerable graph Laplacian, agreeing with equation (2) and

Theorem 2, which predict slopes of −1 and −0.75, respectively, when substituting d = 2. Moreover, the

errors due to the steerable and standard graph Laplacians coincide in the region where the errors are
dominated by the bias error term, also in agreement with Theorem 2.

Additionally, we computed the eigenvalues of L̃ as described in Section 2.3, and compared them with
the eigenvalues of the standard (normalized) graph Laplacian. The results can be seen in Figure 2. It
is evident that the eigenvalues in both cases agree with the well-known multiplicities of the spherical
harmonics (the eigenfunctions of the Laplacian on the unit sphere). However, is clear that the eigenvalues
of L̃ admit smaller fluctuations compared to the eigenvalues of the standard (normalized) graph Laplacian,
owing to the improved convergence rate of L̃ to the Laplace-Beltrami operator.

10

2

Standard Graph-Laplacian
Steerable Graph-Laplacian

1

0

-1

-2

-3

-4

-5
-6

-5

-4

-3

-2

-1

0

1

Figure 1: Errors in approximating ∆M f (x0 ) by the standard graph Laplacian (equations (1) and (2))
and by the steerable graph Laplacian (equation (26)) as a function of ε (in log scale). The region to the
left of the dashed vertical line is dominated by the variance error term, whereas the region to the right of
the dashed vertical line is dominated by the bias error term.

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0
0

10

20

30

40

50

60

0

10

20

30

40

50

60

Figure 2: Eigenvalues of the steerable (left) and standard (right) normalized graph Laplacians, for ε =
1/4 and 2, 000 data points sampled uniformly from the sphere. For the steerable graph Laplacian, the
eigenvalues were sorted in ascending order and enumerated over (m, k) using a single joint index.

11

2.5

Filtering image data-sets by the steerable manifold harmonics

Next, we propose to expand our data-set of images and all of their rotations by a carefully-chosen subset of
the steerable manifold harmonics (the eigenfunctions of the steerable graph Laplacian L̃, see Theorem 3).
Consider the function Fm,ℓ ∈ H given by
Fm,ℓ (i, ϕ) = xϕ
i,(m,ℓ) ,

(28)

where xϕ
i,(m,ℓ) stands for the (m, ℓ)’th coordinate of the i’th data-point rotated by ϕ (via (6)). In essence,
the function Fm,ℓ describes the (m, ℓ)’th coordinate of all points in the data-set and all of their rotations.
As Fm,ℓ ∈ H, it can be expanded in the basis Φ̃, and we can write
∞
N
X
X

Fm,ℓ =

A(m′ ,k),(m,ℓ) Φ̃m′ ,k

(29)

m =−∞ k=1
′

for all (m, ℓ) pairs, where A(m′ ,k),(m,ℓ) are some associated expansion coefficients. We propose to “filter”
the functions Fm,ℓ for each pair (m, ℓ) by considering a truncated expansion of the form of (29), with
expansion coefficients obtained by solving
k

′

Fm,ℓ −

M
X

2

′

m
X

′

−→ min,

A(m′ ,k),(m,ℓ) Φ̃m′ ,k

m =−M k=1
′

A

(30)

H
′

where A is a matrix of expansion coefficients with rows indexed by (m , k) and columns indexed by (m, ℓ).

′
As for the numbers of chosen basis functions km′ and M , we propose the following natural truncation rule based on a cut-off frequency λc ∈ R+ :

where

n

λ̃m,k

o

o
n
km′ = max k : λ̃m′ ,k < λc ,

(31)

are the eigenvalues (sorted in non-decreasing order w.r.t k) of the normalized steerable
′

graph Laplacian L̃. Then, M is simply the largest m

′

s.t. km′ > 0. Fundamentally, this truncation

rule can be viewed as the analogue of the classical truncation of Fourier expansions. Figure 3 illustrates
′

a typical configuration of index-pairs (m , k) resulting from the truncation rule of (31).
We motivate the above-mentioned approach (series expansion and truncation rule) as follows. It is
well known that for smooth and compact manifolds the Laplace-Beltrami operator ∆M admits a sequence
∞
of eigenvalues {µk }∞
k=0 and eigenfunctions {φk }k=0 , which are orthogonal and complete in the class of

square-integrable functions on M, denoted by L2 (M). Therefore, every function f ∈ L2 (M) can be
expanded as

f (x) =

∞
X

ak φk (x),

ak =

k=1

Z

M

f (x)φ∗k (x)dx.

(32)

In this context, it is possible to consider the coordinates of M in the ambient space, i.e. xm,ℓ for every
x ∈ M, as smooth functions over M, which can be approximated by truncating the above-mentioned

expansion. In particular, we provide the following proposition, which bounds the error in approximating
12

the coordinates of M using a truncated series of manifold harmonics.
∞
Proposition 4. Let {φk }∞
k=1 and {µk }k=1 be the eigenfunctions and eigenvalues (sorted in non-decreasing

order), respectively, of the negative-defined Laplace-Beltrami operator ∆M . Then, we have that
1
Vol {M}

Z

M

n
X
k=1

2

ck φk (x) − x

2

dx ≤

d
µn+1

,

ck =

Z

M

xφ∗k (x)dx,

(33)

where Vol {M} is the volume of M, d is the intrinsic dimension of M, and ck ∈ CD are the expansion

coefficients of x ∈ CD (i.e. of every coordinate function of the embedded manifold) w.r.t φk .

Proof. The proof follows immediately from combining equation (3.1) in [1] and Proposition 3.1 in [21].
It is important to note that by the properties of ∆M we have that µn →∞ when n → ∞ [24], and

therefore we can get an arbitrarily small approximation error for the coordinates of M using a sufficiently

large number of manifold harmonics. As we have shown in Section (2.3) that L̃ approximates the Laplace-

Beltrami
operator
n
o
n ∆M ,owe follow the common practice and use the eigenfunctions and eigenvalues of L̃,
i.e. Φ̃m,k and λ̃m,k , as discrete proxies for {φk } and {µk } in (33).

Next, we proceed to derive a simple and efficient solution to problem (30). By our construction of the

Hilbert space H, one can write (30) explicitly as
N Z
X
i=1

2π
0

k

′

xϕ
i,(m,ℓ) −

M
X

m =−M k=1
′

′

2

′

m
X

A(m′ ,k),(m,ℓ) Φ̃m′ ,k (i, ϕ) dϕ −→ min,
A

(34)

which is interpreted as performing regression over the entire data-set of images
and all oof their planar
n

′
′
rotations using the functions Φ̃m,k restricted to k ∈ 1, . . . , km′ , m ∈ −M , . . . , M . Recall that

by (6), we have that

ımϕ
,
xϕ
i,(m,ℓ) = xi,(m,ℓ) e

(35)

where xi,(m,ℓ) stands for the (m, ℓ)’th coordinate of the i’th data-point. It turns out that (34) can be
significantly simplified by substituting (35) into (34) together with the steerable form of Φ̃m,ℓ (i.e. (19)),
while making use of the orthogonality of the Fourier modes {eımϕ }∞
m=−∞ over [0, 2π). It then immediately

follows that the matrix of coefficients A in (34) is block-diagonal, where the blocks can be obtained by
solving ordinary least-squares problems. In particular, we have that

A(m′ ,k),(m,ℓ) =


B (m) , m = m′ ,
k,ℓ
0,

′

m 6= m ,

(36)

where B (m) is the m’th block on the diagonal of A, obtained by solving the least-squares system
X (m) − Ve (m) B (m)

13

2

−→ min,

F B (m)

(37)

where k·kF stands for the Frobenius norm, and X (m) and Ve (m) are given by
X (m)




x1,(m,1) . . . x1,(m,ℓm )


..
..
..
,
=
.
.
.


xN,(m,1) . . . xN,(m,ℓm )



|


Ve (m) = ṽm,1 · · ·
|

|




ṽm,km  ,
|

(38)

′

with ṽm,k given by (19) and (20). We mention that km′ changes with the angular index m , and in
particular, is typically smaller for higher angular frequencies (larger |m|). Therefore, the size of the

blocks B (m) reduces with |m|, as illustrated by Figure 4. Once the coefficients matrices B (m) were
obtained by solving (37), we define

X̂ (m) , Ve (m) B (m)

(39)

as the filtered data-set corresponding to the angular index m.

140
200
120

400
600

100

800
80
1000
60

1200
1400

40
1600
20
1800
2000
0

10

20

100 200 300 400 500

Figure 3: A typical configuration of index pairs
′
(m , k) obtained by the truncation rule of (31), for
non-negative angular indices m.

Figure 4: A typical structure of the matrix
′
A for non-negative angular indices m and m .
Black-coloured regions correspond non-zero entries, where each block corresponds to a different
′
angular index m = m .

Lastly, a favorable interpretation of this procedure can be derived as follows. If we denote by Q(m)
a matrix whose columns are orthonormal and span the columns of Ve (m) , then X̂ (m) can be written
equivalently as

i∗
h
X̂ (m) = Q(m) Q(m) X (m) = C (m) X (m) ,

14

(40)

where (·)∗ denotes complex-conjugate and transpose, and we defined the N × N “filtering” matrices
i∗
h
C (m) , Q(m) Q(m) ,

(41)

which are applied to our data-set for every angular index separately. Essentially, C (m) acts as a “lowpass filter”, in the sense that it retains only the contribution of steerable manifold harmonics with low
frequencies (i.e. eigenvalues below the threshold λc ). In this context, the cut-off frequency λc controls
the rank of C (m) , which is equal to km , and the degree to which C (m) suppresses oscillations in the data.

3

Algorithms summary and computational cost

We outline the algorithms for evaluating the steerable manifold harmonics and employing them for filtering
image data-sets in Algorithms 1 and 2, respectively. We note that a modification to Algorithm 1 is
proposed in Section 4 for improving its robustness to noise. We added this modification to the algorithm
(step 3) under the label “Implicit debiasing (optional)”.
Algorithm 1 Evaluating the steerable manifold harmonics
Required: A data-set of N points {x1 , . . . , xN } ⊂ CD , where xi,(m,ℓ) is the (m, ℓ)’th coordinate of xi
(see Section 2.1).
1: Choose a numerical-integration parameter K (see (16)), and a Gaussian kernel parameter ε.
2: For every 1 ≤ i, j ≤ N , k ∈ {0, . . . , K − 1}, compute the affinities
)
(
ℓm
M
X
X
2
ı2πmk/K
k
xi,(m,ℓ) − xj,(m,ℓ)e
/ε .
(42)
Wi,j = exp −
m=−M ℓ=1

k = 0 for 1 ≤ i ≤ N and k = 0, . . . , K − 1.
Implicit debiasing (optional): Set Wi,i
4: For every angular index m = −M, . . . , M and 1 ≤ i, j ≤ N , evaluate

3:

m
Ŵi,j
=

K−1
2π X k ı2πmk/K
Wi,j e
,
K

Di =

k=0

5:

N
X

0
Ŵi,j
.

(43)

j=1

For every angular index m = −M, . . . , M form the matrix
S̃m = I − D −1 Ŵ (m) ,

(44)

N
where D is a diagonal matrix with {Di }N
i=1 on its diagonal, and return its eigenvectors {ṽm,k }k=1 and
n
oN
eigenvalues λ̃m,k
.
k=1

We now turn our attention to the computational complexity of Algorithms 1 and 2. We begin with
k , which can be evaluated efficiently by
Algorithm 1. The first step is to compute all affinity measures Wi,j

the FFT if we notice that
ℓm
M
X
X

m=−M ℓ=1

xi,(m,ℓ) − xj,(m,ℓ)eı2πmk/K

2

= kxi k22 + kxj k22 − 2 Re

15

(

M
X

m=−M

(m)
ci,j e−ı2πmk/K

)

,

(45)

Algorithm 2 Rotationally-invariant data-set filtering
Required:
(a) A data-set of N points {x1 , . . . , xN } ⊂ CD , where xi,(m,ℓ) is the (m, ℓ)’th coordinate of xi (see
Section 2.1).
n oM
n
oM
(b) Eigenvectors {ṽm,1 , . . . , ṽm,N }M
of
and
eigenvalues
λ̃
,
.
.
.
,
λ̃
S̃m
m,1
m,N
m=−M
m=−M

from Algorithm 1.

1:
2:

3:

m=−M

Choose a cut-off frequency λc .
For m = −M, . . . , M do
o
n
(a) Compute km = max k : λ̃m,k < λc , and form the matrices X (m) and Ve (m) of (38).
(b) Estimate the coefficients matrix B (m) by solving the least squares system of (37).
(c) Compute X̂ (m) = Ve (m) B (m) .

The filtered data-set is given by X̂ = [X̂ (−M ) · · · X̂ (M ) ].

where we defined
(m)

ci,j =

ℓm
X

xi,(m,ℓ) x∗j,(m,ℓ) .

(46)

ℓ=1

(m)

Note that computing ci,j for all i, j, m takes O(N 2 D) operations. Therefore, if we denote
M̄ = max {M, K},

(47)

the computational complexity of this step is O(N 2 D + N 2 M̄ log M̄ ) when using the FFT to compute (45).
m (of step 4) by the FFT takes O(N 2 M̄ log M̄ ) operations. Lastly,
In a similar fashion, computing Ŵi,j
n oM
forming the matrices S̃m
requires O(M N 2 ) operations, and evaluating its eigenvectors and eigenm=−M

values takes O(M N 3 ) operations. Overall, the computational complexity of Algorithm 1 is therefore

O M N 3 + N 2 (D + M̄ log M̄ ) .

(48)

In practice, it is often the case that only a small fraction of pairs of indices (i, j) contributes significantly to
k , since only images which are similar up to a planar rotation admit a non-negligible affinity (assuming
Wi,j
k , allowing for
that ε is sufficiently small). Hence, it is often possible to zero-out the small values of Wi,j

cheaper sparse-matrix computations. Additionally, computing the eigen-decomposition in step 5 for large
data-sets (large N ) may be accomplished more efficiently using randomized methods [15, 2].
2
−→ min, where X (m)
As of Algorithm 2, in part (b) of step 2 we need to solve X (m) − Ve (m) B (m)
F B (m)

is of dimension N × ℓm , Ve (m) is N × km , and B (m) is km × ℓm . For each angular index m, this amounts to
solving ℓm least-squares problems (one for each column of B (m) ), each of size of size N ×km. Assuming that
2 N +N k ℓ +k 2 ℓ )
N ≥ km and using the QR factorization to solve least-squares, this part requires O(km
m m
m m
2 N ) operations to compute the QR decomposition
operations for every angular index m, as we need O(km

of V (m) (which needs to be computed only once), O(km ℓm N ) operations to apply Q of the QR to X (m) ,
2 ℓ ) operations to solve the resulting ℓ triangular systems. Then, since part (c) of step 2 takes
and O(km
m
m

16

O(N km ℓm ) operations for each m, it follows that Algorithm 2 requires
O

M
X

2
2
km
N + N km ℓm + km
ℓm

m=−M

!


= O M N k̄2 + D k̄(N + k̄)

(49)

operations, where k̄ = maxm {km }. As typically k̄ << N , the computational cost of Algorithm 2 is

negligible compared to that of Algorithm 1.

4

Analysis under Gaussian noise

Next, we analyze our method under Gaussian white noise, and argue that in a certain sense the steerable
graph Laplacian (following a certain modification) is robust to noise. Moreover, we argue that the filtering
procedure (described in Section 2.5) allows us to reduce the amount of noise in the filtered data-set
proportionally to the number of images N .
Suppose that we are given a noisy data-set
yi,(m,ℓ) = xi,(m,ℓ) + ηi,(m,ℓ) ,

(50)


where xi,(m,ℓ) is the (m, ℓ)’th coordinate of the i’th clean data point, and ηi,(m,ℓ) are independent and
normally distributed complex-valued noise variables with mean zero and variance σ 2 .

We start by considering the effect of noise on the construction of L̃ (of (17)). Clearly, the noise
changes the pair-wise distances computed in Wi,j (ϑ, ϕ), where we note that from Theorem 3 and (13) it is
−1
sufficient to consider the effect of noise only on Di,i
Wi,j (0, α), for all α ∈ [0, 2π) and i, j = 1, . . . , N . To

α
α
α , y , yα , · · · , yα
D
this end, consider the set of points Yi = y1 , · · · , yi−1
i i+1
N ⊂ C , where all points except

the i’th were replaced with their rotations by an angle α (via (6)). We have that
yjα = xαj + ηjα ,

α
yj,(m,ℓ)
= xj,(m,ℓ) eımα + ηj,(m,ℓ) eımα ,


α , η , ηα , · · · , ηα
for j 6= i, j = 1, . . . , N , and it is evident that the set of noise points η1α , · · · , ηi−1
i i+1
N

(51)
are

still i.i.d Gaussian. Then, Theorem 2.1 (and specifically equation (1)) in [11], when applied to the set Yiα ,

asserts that if we denote γ = Dσ 2 and vary D and σ 2 such that γ remains constant, then
n
Wi,j (0, α) = exp − yi − yjα

2

o
n 
/ε −→ exp − xi − xαj
D→∞

2

 o
+ 2γ /ε

(52)

in probability, for all j 6= i. Essentially, this result is due to the concentration of measure of high-

dimensional Gaussian random vectors, and in particular the fact that {ηiα } are uncorrelated and are

concentrated around the surface of a sphere in CD . Therefore, in the regime of high dimensionality and
small noise-variance, the effect of the noise on the pair-wise distances (between different data-points and
rotations) is only an additive constant bias term. Note that even though the noise variance tends to
zero, the overall noise magnitude γ = Dσ 2 is kept constant and may be large, corresponding to a low
signal-to-noise ratio (SNR). We further mention that this constant-bias effect is not restricted to Gaussian

white noise, as it occurs also when the noise admits a general covariance matrix Σ, and even when the
noise takes other certain non-Gaussian distributions (see [11] for specific details and conditions).
17

Next, in order to correct for the bias in the distances, we follow [12] and zero-out the diagonal of W ,
that is, we update
Wi,j (0, α) ←

(

Wi,j (0, α),

i 6= j,

0,

i = j.

(53)

Then, we expect D −1 to correct (implicitly) for the bias in W , since
Z
N
X

Di,i =

2π

Wi,j (0, α)dα =

D→∞

2π

j6=i,j=1 0

j6=i,j=1 0

−→ e−2γ/ε

Z
N
X

Z
N
X

2π

j6=i,j=1 0

n
exp − xi − xαj

2

n
exp − yi − yjα

o
/ε dα,

2

o
/ε dα
(54)

in probability, and thus


2
α
exp − xi − xj /ε
−1


Wi,j (0, α) −→
Di,i
R 2π
′ 2
D→∞ PN
α
/ε dα′
j6=i,j=1 0 exp − xi − xj

(55)

in probability, which is equivalent to its clean counterpart for i 6= j (after zeroing-out the diagonal).

Lastly, we argue that zeroing-out the diagonal of W does not change the point-wise convergence rate
of the clean steerable graph Laplacian (as reported by Theorem 2), as it results in an error which is
negligible compared to the leading error terms (see the end of Section E.2 in the proof of Theorem 2,
and [26] for an analogous argument in the case of the standard graph Laplacian). In summary, the above
discussion suggests that when the dimension D is large, the noise variance σ 2 is small, and the overall

noise magnitude γ = Dσ 2 is fixed (and may be large compared to the magnitude of the signal), the

steerable graph Laplacian constructed from the noisy data after implicit de-biasing (by zeroing-out the
diagonal of W ) is expected to be close to its clean counterpart.
Next, we consider the eigenvectors and eigenvalues computed from the clean (normalized) steerable
graph Laplacian, and analyze the affect of the filtering procedure (described in Section 2.5) on the noise
in the data-set. From (40), the de-noised data for angular frequency m is given by
i∗
h
e (m) = C (m) Y (m) = Q(m) Q(m) Y (m) ,
X

(56)
(m)

where Y (m) is the matrix of noisy data points corresponding to angular index m, i.e. Yi,ℓ

= yi,(m,ℓ) . Since

Q(m)

consists of km orthonormal column vectors independent of the noise, and recalling that γ = Dσ 2 ,
PM
D = m=−M ℓm , we have that
1
e
E X̂ − X
N

2

F

1
=
N

M
X

m=−M

E X̂

(m)

e (m)
−X

2

F

=

σ2

PM

m=−M

N

km ℓm

≤ max {km }
m

γ
,
N

(57)

e =
where X̂ = [X̂ (−M ) · · · X̂ (M ) ] (X̂ (m) is defined in (40)) represents the clean filtered data-set, and X
e (−M ) · · · X
e (M ) ] represents the noisy filtered data-set. Hence, larger data-sets are expected to provide
[X
18

e reduces proportionally to 1/N .
improved de-noising results, as the noise in the filtered data-set X

At this point, it is worthwhile to point out that the error bound of (57) is significantly better then what

we would expect from using the standard graph Laplacian and its eigenvectors (to filter the coordinates
of the data-set). Fundamentally, this is due to the block diagonal structure of the coefficients matrix A
(see Section 2.5), and more specifically, the fact that we only need to use eigenfunctions with angular
index m to expand data coordinates with the same angular index, in contrast to using all eigenfunctions.
In particular, since the limiting operators of the steerable and standard graph Laplacians are the same
(the Laplace-Beltrami operator), we expect the truncation rule of (31) to provide a similar number of
eigenfunctions/eigenvectors from both methods. Then, if we were to use the eigenvectors of the standard
P
graph Laplacian to filter our data-set, we would be required to use all ∼ M
m=−M km eigenvectors, and
PM
by a computation equivalent to (57) we would expect an error of ∼ m=−M km Nγ , which is considerably
larger than maxm {km } Nγ . In conclusion, as the steerable graph Laplacian is more informative than the

standard graph Laplacian, in the sense that it provides us with the angular part of each eigenfunction,

it allows us to be more precise when filtering our data-set by matching the angular frequencies of the
eigenfunctions to those of the data, thereby reducing the computational complexity and improving the
de-noising performance considerably. This feature of the steerable graph Laplacian stands on its own,
and is separate from the improved convergence rate to the Laplace-Beltrami operator (Theorem 2), which
improves the accuracy of the eigenfunctions and eigenvalues compared to those of the standard graph
Laplacian.
Lastly, we mention that the error term (57) can be viewed as a variance error term in a classical
bias-variance trade-off, as we can write, conditioned on the clean data-set X, that
=X̂

1
e
E X −X
N

2
F

=X̂

z h}| i{ z h}| i{
2
1
e +E X
e −X
e
= E X −E X
N
F
2
2
1
1
e .
= E X − X̂ + E X̂ − X
F
F
|N
{z
} |N
{z
}
Bias

(58)

Variance

Consequently, the overall error cannot get arbitrarily small, as there exists a bias term when approximating
the clean data points by finitely many eigenfunctions (see Proposition 4). Therefore, in practice, the
optimal de-noising results for a given data-set and noise variance would be attained as an optimum in a
bias-variance trade-off, where a large cut-off frequency λc would result in larger {km } values and a larger
variance error (as noise is mapped to more expansion coefficients), and a smaller cut-off frequency λc
would result in smaller {km } values and a larger bias error.
Remark 4. While the discussion in this section suggests that our method is robust to noise in the highdimensional regime, it is not to say that reducing the dimensionality of a given data-set (with a given and
fixed noise variance σ 2 ) would degrade the accuracy of the quantities computed by our method. On the
contrary, a close examination of the results in [11] reveals that the errors in pair-wise distances computed
√
from noisy data points are dominated by Dσ 2 , meaning that projecting the data onto a lower-dimensional

subspace (while retaining a sufficient approximation accuracy w.r.t the clean data) is encouraged – as it
improves the accuracy of the pair-wise affinities on one hand, and reduces the overall noise magnitude
γ = Dσ 2 on the other.
19

5

Example: De-noising cryo-EM projection images

In this section, we demonstrate how we can use our framework to de-noise single-particle cryo-electron
microscopy (cryo-EM) image data-sets.

5.1

Cryo-EM

In single-particle cryo-EM [13, 6], one is interested in reconstructing a three-dimensional model of a
macromolecule (such as a protein) from a set of two-dimensional images taken by an electron microscope.
The procedure begins by embedding many copies of the macromolecule in a thin layer of ice, where due to
the experimental set-up, the different copies are frozen at random unknown orientations. Then, an electron
microscope acquires two-dimensional images of the these macromolecules (more precisely, it samples the
Radon transform of the density function of the macromolecule). Consequently, it can be shown that the
set of all projection images lies on a three-dimensional manifold diffeomorphic to the group SO(3). Thus,
the manifold model assumption discussed in this work is natural for describing cryo-EM data-sets. Note
that due to the experimental set-up in cryo-EM, the in-plane rotation of each copy of the macromolecule
is arbitrary, and therefore, so are the planar rotations of the two-dimensional images. Additionally, the
images acquired in cryo-EM experiments are very noisy, with a typical SNR (Signal-to-Noise Ratio) of
1/10 and lower. Simulated clean and noisy cryo-EM images of the 70S ribosome subunit can be seen in
Figure 6 (top two rows).

5.2

De-noising recipe

Given a collection of cryo-EM projection images {I1 , . . . , IN } sampled on a Cartesian grid, we start by
performing steerable principal components analysis (sPCA), as described in [18]. This procedure provides
us with steerable basis functions (the steerable principal components) {ψm,ℓ } of the form of (3), which are
optimal for expanding the images in the data-set and all of their rotations. For each basis function ψm,ℓ ,

the steerable PCA also returns its associated eigenvalue νm,k , which encodes the contribution of ψm,ℓ to
the expansion (analogously to the eigenvalues of the covariance matrix in standard PCA). Therefore, we
have that
Ii ≈

ℓm
M
X
X

yi,(m,ℓ) ψm,ℓ ,

(59)

m=−M ℓ=1

where yi,(m,ℓ) is the (m, ℓ)’th expansion coefficient of the i’th image (provided by sPCA, see [18] for
appropriate error bounds associated with (59)). Expanding the image data-set using such basis functions
allows us to apply our filtering scheme in the domain of the expansion coefficients, as required by our
algorithms. We note that for images corrupted by additive white Gaussian noise, the noise variance σ 2 is
estimated from the corners of the images (where no molecule is expected to be present), and the number
of basis functions used in the expansion, governed by M and {ℓm }, is determined by estimating which

eigenvalues νm,k are above the noise level (i.e. exceed the Baik-Ben Arous-Péché transition point [3], see
also [35, 34]) via

(



ℓm = max ℓ : νm,ℓ > σ 2 1 +

20

r

nm
N

2 )

,

(60)

Data-set Y

Steerable PCA [18]

Steerable GraphLaplacian (Alg. 1)

Input images

Eigenvectorsn{ṽm,ko
}
λ̃m,k

Basis {ψm,ℓ }

P
Iˆi = m,ℓ x̂i,(m,ℓ) ψm,ℓ

Eigenvalues

Filtered data-set X̂

Rotationally-invariant
filtering (Alg. 2)

Expansion in basis
Filtered images

Figure 5: Schemtaic view of the de-noising procedure. We start by applying the steerable
 PCA [18] to the
input images, obtaining basis functions {ψm,ℓ } and associated expansion coefficients yi,(m,ℓ) (organized
into the matrix Y following the layout of (4)), where the truncation of the expansion is due to (60). Then,
the data-set Y is used to construct the steerable graph Laplacian, whose eigenfunctions (the steerable
manifold harmonics) are obtained via Algorithm 1 with implicit debiasing (step 3). Using the
 steerable
manifold hamronics, we filter the data-set via Algorithm 2, and use the filtered coefficients x̂i,(m,ℓ) in
conjunction with the basis functions {ψm,ℓ } to get back the filtered images.
where {nm } can be found in [18] (nm is the size of the m’th block in the block-diagonal covariance matrix

associated with steerable PCA), and assuming that {νm,k }k are sorted in a non-increasing order for every

m. Correspondingly, M in (59) is simply the largest |m| s.t. ℓm > 0.

Using the above setting, the task of de-noising the images {Ii } is reduced to the task of de-noising the

sPCA coefficients yi,(m,ℓ) . We then estimate the steerable manifold harmonics {ṽm,k } (as described by

Algorithm 1) from the data-set {yi }N
i=1 , and follow by employing {ṽm,k } for filtering the data-set according

to Algorithm 2. After obtaining the de-noised expansion coefficients x̂i,(m,ℓ) , we can plug them back
in the expansion (3) to get de-noised images. The procedure is summarized in Figure 5. Note that since
{Ii } are real-valued images, their expansion coefficients satisfy the symmetry
∗
yi,(−m,ℓ) = yi,(m,ℓ)
.

(61)

Therefore, it is sufficient to de-noise only the coefficients with non-negative angular frequencies.

5.3

Experimental results

We demonstrate the de-noising performance of our approach using simulated images of the 70S ribosome,
of size 128 × 128 pixels, after applying a filter to all images corresponding to a typical Contrast Transfer

Function (CTF) [13] of the electron microscope. As described previously, we first map all images to their

sPCA coefficients via [18] (with T = 10 and half-Nyquist bandlimit), and then proceed according to our
filtering scheme (Algorithms 1 and 2). We mention that throughout our experiments the choice K = 256
was found satisfactory, and that ε and λc were chosen automatically for every experimental set-up (determined by the number of images N and noise variance σ 2 ) as described in Appendix A. In every experiment,

21

we compare the de-noised images resulting from our method to the images obtained directly from the

sPCA coefficients (i.e. images computed from the coefficients yi,(m,ℓ) ), and to images obtained after
applying a shrinkage to the sPCA coefficients via yi,(m,ℓ) wm,ℓ , where the weights {wm,1 , . . . , wm,ℓm }M
m=−M ,

which were computed as described in [35], correspond to the asymptotically-optimal Wiener filter [28].
Essentially, this is the optimal filter for the expansion coefficients in the sense of minimizing the mean
squared error.
First, we demonstrate our method on 10, 000 projection images at signal-to-noise ratio of 1/20. The
de-noised images can be seen in Figure 6, where it is visually evident that the final de-noised images using
our method contain many more details compared to sPCA Wiener filtering, which results in somewhat
blurred images due to the aggressive shrinkage of sPCA coefficients. In terms of performance measures,
our method (which we term “sMH filtering”, where sMH stands for “steerable manifold harmonics”)
results in an average peak-SNR (pSNR) of 25.37dB, where the sPCA Wiener filter provided 21.17dB
pSNR, and sPCA alone resulted in 17.64dB pSNR.
It is therefore evident that Wiener filtering of sPCA coefficients is far from optimal in terms of denoising and image recovery, as it essentially applies a single linear operator on the individual images,
which is only optimal when the data resides on a linear subspace. However, in the case of cryo-EM, as the
data resides on a manifold, it is reasonable to apply non-linear methods which account for the geometry
and topology of the manifold. In this respect, our method is able to make use of all images and their
rotations simultaneously to accurately estimate the structure of the manifold, and thereby provides an
improved de-noising of the image data-set.
Next, Figure 7 demonstrates the performance of our method for SNR = 1/10 and different values of
N (data-set size). As anticipated, our method is able to exploit larger data-sets for improved de-noising,
whereas the sPCA Wiener filtering offers only a mild gain beyond 2, 000 images. The reason for that is
that the Wiener filtering is applied to each image separately, and therefore reaches saturation once the
estimation of the sPCA from the noisy data is sufficiently accurate (approaches the sPCA of the clean
data). Note that the pSNR from the projection onto the sPCA components (without shrinkage) reduces
with N , because more basis functions ψm,ℓ are chosen (according to (60)) as N increases, even if their
contribution to expanding the data-set is negligible. Therefore, the dimension D increases, and with it

also the overall noise magnitude Dσ 2 . It is important to mention that even though the variance error
term in (57) behaves like 1/N , the improvement in the pSNR of our method is not expected to follow

this trend, since the overall error also includes a bias error term (see (58)), such that the minimal error
for every value of N is attained as a different optimum in the bias-variance trade-off.
Lastly, we evaluated the de-noising quality for N = 10, 000 images and varying amounts of noise. The
results are displayed in Figure 8, where we can see that our method outperforms the sPCA Wiener filter
considerably in a wide range of SNRs. We remark that as the SNR decreases the asymptotics considered
in Section 4 become less valid, thus at some point the steerable graph Laplacian becomes too noisy, and
the performance gain of our method drops. This phenomenon is mostly evident for SNRs below −14dB,

and our method eventually under-performs the sPCA Wiener filter at −20dB SNR.

22

Figure 6: Images after de-noising, for N = 10, 000 and SNR = 1/20. Different colummns correspond to
different images from the dataset (different in-plane rotations and viewing directions of the molecule),
while different rows correspond to (from top to bottom): clean images, noisy images, sPCA (pSNR =
17.64dB), sPCA Wiener filter (pSNR = 21.17dB), and sMH filtering (this paper, pSNR = 25.37dB).

28

sMH filtering (this paper)
sPCA Wiener filter
sPCA

27
26
25
24
23
22
21
20
19
0

2000

4000

6000

8000

10000

Figure 7: pSNR of de-noised images for SNR = 1/10 and different data-set sizes N .

23

30

sMH filtering (this paper)
sPCA Wiener filter
sPCA

28
26
24
22
20
18
16
14
12
-6

-8

-10

-12

-14

-16

-18

-20

Figure 8: pSNR of de-noised images for N = 10, 000 and various SNR levels.

6

Conclusions and discussion

In this work, we introduced the steerable graph Laplacian, which generalizes the standard graph Laplacian
by incorporating all planar rotations of all images in the data-set. We demonstrated that the (normalized)
steerable graph Laplacian is both more accurate and more informative than the standard graph Laplacian,
in the sense that it allows for an improved approximation of the Laplace-Beltrami operator on one hand,
and admits eigenfunctions with a closed-form analytic expression of their angular part (i.e. the angular
Fourier modes) on the other. This closed-form expression is essentially what allows for the efficient
filtering procedure of the data coordinates (see Section 2.5), as we only need to estimate a block-diagonal
coefficients matrix. Then, we have shown that under a suitable modification, the (normalized) steerable
graph Laplacian is robust to noise in the regime of high dimensionality due to the concentration of
measure of Gaussian noise. Moreover, we have seen that the proposed filtering procedure reduces the
noise proportionally to the number of images in the data-set, which was corroborated by the experiments
of de-noising cryo-EM projection images, where we demonstrated that our method can provide excellent
de-noising results on highly noisy image data-sets.
It is interesting to point-out that the steerable graph Laplacian, while utilized for filtering image
data-sets, can be employed for many other purposes. One application immediately coming to mind is the
filtering of data-sets consisting of periodic signals (see last remark in Section 2.1). However, and more
importantly, the steerable graph Laplacian can replace the standard graph Laplacian in all applications
where the domain is known to be rotationally invariant (by our definition in Section 2.1). For instance,
it can be used for regularization over general signal/data recovery inverse problems, or for dimensionality
reduction in the framework of Diffusion Maps [7] and Laplacian Eigenmaps [4]. In this context, we
mention the method of Vector Diffusion Maps (VDM) [27], which allows for diffusion-based dimensionality
reduction for manifold data in the presence of nuisance parameters (such as planar rotations). However,

24

as VDM computes the group ratios only between pairs of data points (e.g. optimal rotational alignments),
it may be of interest to compare it to the steerable graph Laplacian, which essentially considers all planar
rotations of all images. Lastly, we note that a possible future research direction is the extension of our
techniques to other group actions (such as SO(3)).

Acknowledgements
We would like to thank Amit Singer for useful comments and suggestions regarding this work. This
research was supported by THE ISRAEL SCIENCE FOUNDATION grant No. 578/14, by Award Number
R01GM090200 from the NIGMS, and by the European Research Council (ERC) under the European
Unions Horizon 2020 research and innovation programme (grant agreement 723991 - CRYOMATH).

Appendix A

Choosing ε and λc

In this section, we provide some guidelines on how to choose the parameters ε and λc , and propose a
method for determining them automatically.
In the clean setting, it is reasonable to choose the optimal value of ε by minimizing the bias-variance
error related to the convergence of the steerable graph Laplacian (i.e. equation (18)). However, in general
this cannot be achieved without prior knowledge on geometrical quantities of M, such as its curvature
(see [26] and the related discussion therein). Nonetheless, methods for automatic picking of ε for the

standard graph Laplacian were proposed (see for example [8]). Unfortunately, such methods cannot be
applied to noisy data in a straightforward manner, due to the effects of noise on the pair-wise distances,
and we are not aware of any work detailing the automatic choice of ε for this case. We note however, that
when the distances are computed from the noisy data-points {y1 , . . . , yN }, and if N is sufficiently large

such that the noise is the dominant factor in determining ε, then it is reasonable to choose
ε∝

√ 2
Dσ ,

(62)

as this expression dominates the standard deviation of the errors in the pair-wise (squared) distances
(see [11] and analysis therein), and so can be used to set a region in which it is likely to find the optimal
ε.
In order to find the optimal ε and cut-off frequency λc for a given data-set, we propose to use a
cross-validation procedure. In particular, we choose ε and λc such that the log-likelihood of a noisy set
of points, given a de-noised set (where the sets are disjoint), is maximized. In more detail, suppose that
{y1 , . . . , yN } is a collection of noisy data-points as in (50), i.e.
yi = xi + ηi ,

(63)

where ηi ∈ CD are i.i.d Gaussian vectors with mean zero and covariance σ 2 I, and xi are sampled from the

manifold M with uniform distribution. Then, the log-likelihood of obtaining the subset yN ′ +1 , . . . , yN ,

25

′

where N < N , is given (up to additive constants) by
N
X

i=N ′ +1



1
log
Vol {M}

Z

M

o 
n
2
2
exp − kyi − xk /2σ dx .

(64)

We then propose to approximate this log-likelihood by Monte-Carlo integration via the points

x̂1 (ε, λc ), . . . , x̂N ′ (ε, λc ) ,

(65)


obtained from the de-noising of y1 , . . . , yN ′ using the parameters (ε, λc ). That is, we define the empirical

log-likelihood of yN ′ +1 , . . . , yN (as a function of the parameters ε and λc ) by
N
X

J(ε, λc ) =

′

i=N +1

 ′

N
K−1
X
2π X
2πk/K

log
exp − yi − x̂j
(ε, λc )
K
j=1

k=0

2





/2σ 2 ,

(66)

2πk/K

where x̂j
(ε, λc ) stands for the rotation of x̂j (ε, λc ) by an angle of 2πk/K according to (6). We mention
that summing over the rotations improves the accuracy of the Monte-Carlo integration, as we account for
one of the dimensions of the manifold M in the integration. Finally, we choose the parameters ε and λc
by maximizing J, i.e.



εopt , λopt
= arg max {J(ε, λc )}.
c

(67)

(ε,λc )


Essentially, we expect the realizations of the noisy data-points yN ′ +1 , . . . , yN to be best explained

when the de-noised points x̂1 (ε, λc ), . . . , x̂N ′ (ε, λc ) (and their rotations) lay as close as possible to M,

and we seek the best parameters for achieving that. For example, if we take λc to be too small, then the

de-noising error will be dominated by a bias term as we did not take a sufficient number of components
to represent the features of the clean images accurately. Therefore we would expect the empirical loglikelihood to be small as the de-noised images will not be close to M. On the other hand, if we take λc
to be too large, then over-fitting will occur, in the sense that some features of the noise will be preserved
in the de-noised images. In that case, as the empirical likelihood is computed on a set independent from
the set of de-noised points, we would again expect the empirical log-likelihood to be small.

Appendix B

Quadratic form of L

In what follows, we derive the quadratic form of L appearing in equation (12). First, for f = [f1 (ϑ), . . . , fN (ϑ)]T
we have that
(Lf )(i, ϑ) = Di,i fi (ϑ) −

N Z
X
j=1

2π

Wi,j (ϑ, ϕ)fj (ϕ)dϕ,

(68)

0

and therefore
hf, Lf iH =

N
X
i=1

Di,i

Z

2π
0

2

|fi (ϑ)| dϑ −

N Z
X

i,j=1 0

26

2π

Z

0

2π

fi∗ (ϑ)Wi,j (ϑ, ϕ)fj (ϕ)dϑdϕ.

(69)

Then, if we notice that
Di,i =

N Z
X
j=1

2π

Wi,j (0, α)dα =
0

N Z
X
j=1

2π

Wi,j (ϑ, ϕ)dϕ

(70)

0

for every ϑ ∈ [0, 2π) (by changing the integration parameters), we have that
N
X

Di,i

i=1

Z

2π
0

|fi (ϑ)|2 dϑ =

N Z
X

i,j=1 0

2π

Z

2π
0

Wi,j (ϑ, ϕ) |fi (ϑ)|2 dϑdϕ =

N Z
X

i,j=1 0

2π

Z

2π
0

Wi,j (ϑ, ϕ) |fj (ϕ)|2 dϑdϕ,
(71)

due to the symmetry Wi,j (ϑ, ϕ) = Wj,i (ϕ, ϑ) (see the definition of W in (9)). Finally, using (71) we can
write
hf, Lf iH =
=

Z
N Z
h
i
1 X 2π 2π
Wi,j (ϑ, ϕ) |fi (ϑ)|2 + |fj (ϕ)|2 − fi (ϑ)fj∗ (ϕ) − fi∗ (ϑ)fj (ϕ) dϑdϕ
2
0
0
i,j=1

Z
N Z
1 X 2π 2π
Wi,j (ϑ, ϕ) |fi (ϑ) − fj (ϕ)|2 dϑdϕ.
2
0
0

(72)

i,j=1

Appendix C

Linear rotationally-invariant operators

We start with the following definition of linear rotationally-invariant operators over H.
Definition 5 (LRI operators). An operator G : H → H is linear rotationally-invariant (LRI) over H,
if

1. For any fixed (i, ϑ) ∈ Γ (where Γ defined in Section 2.2), the functional {Gf } (i, ϑ) is linear and
continuous in f .

2. G satisfies
{Gf } (i, ϑ − α) = {Gf α } (i, ϑ),

f α (i, ϑ) , fi (ϑ − α),

(73)

for all f ∈ H, α ∈ [0, 2π), and (i, ϑ) ∈ Γ.
In the first requirement of Definition 5, the continuity property essentially means that if f1 and f2 are
close (in H), then Gf1 (i, ϑ) and Gf2 (i, ϑ) are also close (in absolute value). As for the second requirement

(rotational-invariance), loosely speaking, it means that shifting the output of the operator cyclically by
an angle α is equivalent to shifting the input by α, hence the action of the operator itself does not depend
on the angle ϑ (of (i, ϑ) ∈ Γ). We mention that this property of linearity and rotational-invariance can be
viewed as analogous to that of Linear-Time Invariant (LTI) operators, native to classical signal processing.
Remark 5. Our definition of LRI operators is somewhat more restrictive than the name suggests (compared also to classical LTI operators) because of our requirement for continuity of every functional {G·} (i, ϑ).
We note that while this requirement can be removed, allowing for a broader class of operators, it is simpler
to handle and sufficient for our purposes.
The next lemma characterizes the form of LRI operators explicitly.
27

Lemma 6 (Explicit form of LRI operators). Let G be an LRI operator over H. Then, there exist
2 1
unique {Gi,j }N
i,j=1 ∈ L (S ), s.t. for any f ∈ H and (i, ϑ) ∈ Γ

{Gf } (i, ϑ) =

N Z
X
j=1

2π

Gi,j (ϕ − ϑ)fj (ϕ)dϕ.

0

(74)

Proof. By the Riesz representation theorem, if {G·} (i, ϑ) is a linear and continuous functional over H,

then there exists a unique g(i,ϑ) ∈ H such that

{Gf } (i, ϑ) = g(i,ϑ) , f

H

=

N Z
X
j=1

2π

0

∗
gi,j
(ϑ, ϕ)fj (ϕ)dϕ,

(75)

where g(i,ϑ) , [gi,1 (ϑ, ·), . . . , gi,N (ϑ, ·)]T . Additionally, from (73), we have that
α

{Gf } (i, ϑ) = {Gf } (i, ϑ + α − α) = {Gf } (i, ϑ + α) =

N Z
X
j=1

2π

0

∗
gi,j
(ϑ + α, ϕ + α)fj (ϕ )ϕ ,
′

′

′

(76)

′

when changing the integration parameter via ϕ = ϕ−α. Lastly, taking α = −ϑ and defining Gi,j (ϕ−ϑ) =

∗ (0, ϕ − ϑ) concludes the proof.
gi,j

The main contribution of Lemma 6 is to point out that every LRI operator can be characterized by
a finite number of functions {Gi,j }N
in a Fourier series. Therefore, G can be
i,j=1 which can be expanded
n
o∞
defined by
mapped to (and described by) a sequence of matrices Ĝ(m)
m=−∞

(m)
Ĝi,j

=

Z

2π

Gi,j (α)eımα dα.

(77)

0

It is important to notice that from (6) and (9), it immediately follows that the steerable affinity
operator W of (9) is LRI with Gi,j = Wi,j (0, ϑ − ϕ). However, we note that for any fixed (i, ϑ) ∈ Γ, it is

evident that Df (i, ϑ) (where D is the diagonal matrix defined in (11)) is not a continuous functional of

f (in H), as small perturbations in f may lead to arbitrarily large changes in Df (i, ϑ) as it depends on
point-wise values of f . Therefore, the steerable graph Laplacian L (from (11)) is not LRI. Nonetheless,
as we shall see next, we can still employ the properties of LRI operators to characterize a broader family
of operators, which includes the steerable affinity operator W and the steerable graph Laplacian L as
special cases, with eigenfunctions admitting a particularly convenient form. This is the subject of the
next proposition.
Proposition 7. Consider an operator H : H → H of the form
H = A + G,

(78)

where G is LRI and A ∈ CN ×N is a complex-valued matrix. If (λ, v) is an eigenvalue-eigenvector pair of
the matrix

Ĥ (m) = A + Ĝ(m) ,
28

(79)

where Ĝ(m) is from (77), then Φ = v · eımϑ is an eigenfunction of H with eigenvalue λ.
Proof. The proof follows directly from Lemma 6 and the Fourier expansion of Gi,j . Let us write
{HΦ} (i, ϑ) =

N
X

Ai,j Φ(j, ϑ) +

N Z
X
j=1

j=1

2π
0

Gi,j (ϕ − ϑ)Φ(j, ϕ)dϕ,

(80)

where we have used the explicit form of the LRI operator G given by (74). Then, if we expand Gi,j (·) in
a Fourier series as
∞
1 X
(m)
Gi,j (ϕ − ϑ) =
Ĝ e−ım(ϕ−ϑ) ,
2π m=−∞ i,j

(m)
Ĝi,j

=

′

′

Z

2π

Gi,j (α)eımα dα,

(81)

0

we have that
{HΦ} (i, ϑ) =

N
X
j=1

N
X
1
Ai,j Φ(j, ϑ) +
2π
j=1

∞
X

ım ϑ
Ĝm
i,j e

′

m =−∞

Z

2π

′

Φ(j, ϕ)e−ım ϕ dϕ.

(82)

0

Therefore, by substituting Φ(i, ϑ) = vi eımϑ , where vi stands for i’th element of v, we get
{HΦ} (i, ϑ) = eımϑ

where we have used the orthogonality of eımϑ

N
X

Ai,j vj + eımϑ

Ĝm
i,j vj ,

(83)

j=1

j=1

∞
m=−∞

N
X

over [0, 2π). Finally, it follows that

h
i
HΦ = (A + Ĝ(m) )v eımϑ = λveımϑ = λΦ,

(84)

since v is an eigenvector of A + Ĝ(m) with eigenvalue λ.
Therefore, even though the steerable graph Laplacian L is not strictly LRI according to Definition 5,
it still takes the form of the operators considered in Proposition
n 7, and
o∞consequently, we can derive its
(m)
of (13).
eigen-decomposition by making use of the sequence of matrices Ŵ
m=−∞

Appendix D

Proof of Theorem 1

Proof. First, we note that by (15) and (9), it follows directly that Ŵ (m) is Hermitian, which implies that
Sm is also Hermitian, and therefore can be diagonalized by a set of orthogonal eigenvectors. Next, as L is
of the form A + G as required by Proposition 7 (where A is a matrix and G is an LRI operator), we can
obtain a sequence of eigenfunctions and eigenvalues of L by diagonalizing the matrices Sm = D − Ŵ (m) for

every m ∈ Z. Then, the eigenvalues {λm,k }N
k=1 must be real-valued (since Sm is Hermitian), and moreover,

by the quadratic form of L (12) it follows that hf, Lf iH ≥ 0, which implies that L is semi-positive definite
and thus the eigenvalues {λm,k }N
k=1 are non-negative.

Lastly, the fact that {Φm,k }m,k are orthogonal and complete follows from the orthogonality and com
∞
N
pleteness of eımϑ m=−∞ over L2 (S1 ), and the orthogonality and completeness of {vm,k }N
k=1 over C
29

(since D − Ŵm is Hermitian) for every m ∈ Z. In particular, it easily follows that we can expand every
f ∈ H as

f (i, ϑ) =

∞
X

αim eımϑ =

N
∞ X
X

βm,j vi,(m,j) eımϑ ,

(85)

m=∞ j=1

m=∞

where vi,(m,j) stands for the i’th element of the vector vm,j (which is the j’th eigenvector of Sm =
D − Ŵ (m) ), which concludes the proof.

Appendix E
E.1

Proof of Theorem 2

The limit and bias terms

By (17) and (11), we can write


N Z 2π
n
o
X
4
4
−1

L̃g (i, ϑ) = f (xϑi ) −
Di,i
Wi,j (ϑ, ϕ)f (xϕ
j )dϕ
ε
ε
0
j=1




2
R 2π
ϕ
ϕ
1 PN
ϑ
exp
−
x
−
x
/ε
f
(x
)dϕ
i
j=1 0
j
j
N

4
ϑ
.


= 
f
(x
)
−
i

2
R 2π
ε
ϕ
1 PN
ϑ
/ε dϕ
j=1 0 exp − xi − xj
N

(86)

We begin by deriving the limit of (86) for N → ∞ and a fixed ε > 0, showing that it is essentially

the Laplace-Beltrami operator ∆M with an additional bias error term of O(ε). First, let us focus our

attention on the expression
1
Ci,N
(ϑ)


N Z
1 X 2π
,
exp − xϑi − xϕ
j
N
0
j=1

2


/ε f (xϕ
j )dϕ,

(87)

which is the numerator of the second term of (86) (inside the brackets). Before we proceed with the
evaluation of the expression in (87), we construct a convenient parametrization of the manifold M. To

this end, since our manifold M is rotationally-invariant, it would be beneficial to parametrize it by
a rotationally-invariant coordinate z coupled with a rotation angle β ∈ [0, 2π) (analogously to polar

coordinates in R2 ). In particular, in Section E.3 we construct a parametrization x 7→ (z, β) for every
′

′

x ∈ M , where M ⊂ M is a certain smooth neighbourhood of xϑi (defined explicitly in Section E.3),

such that

x = R(z, β) = z β ,

z ∈ N,

β ∈ [0, 2π),

(88)

′

and N ⊂ M is a smooth (d − 1)-dimensional submanifold.

Next, let us continue with the evaluation of (87), and define
Hiϑ (x)

,

Z

2π
0


exp − xϑi − xϕ

30

2


/ε f (xϕ )dϕ.

(89)

′

Now, using the rotationally-invariant parametrization x 7→ (z, β) for every x ∈ M , we can write
xϕ = R(x, ϕ) = R(z β , ϕ) = z ϕ+β ,

(90)

and thus
Hiϑ (x)

=

Z

2π
0

establishing that


exp − xϑi − z ϕ+β
Hiϑ

2


Z
ϕ+β
/ε f (z
)dϕ =

2π

0


exp − xϑi − z ϕ

2


/ε f (z ϕ )dϕ = Hiϑ (z),

(91)

′

is only a function of z ∈ N for all points x ∈ M . Now, since {xi } are i.i.d samples

from M, then by the law of large numbers

N
1 X ϑ
1
Hi (xj ) = lim
N →∞ N
N →∞ N

N
X

1
lim Ci,N
(ϑ) = lim

N →∞

j=1

j6=i,j=1

h
i Z
Hiϑ (xj ) = E Hiϑ (x) =

M

Hiϑ (x)p(x)dx.

(92)

′

By our construction of the rotationally-invariant parametrization, and in particular the set M (see
Section (E.3)), we have that xϑi − x
Z

M

2

′

′

′

> (δ )2 for all x ∈
/ M and some positive constant δ > 0. Hence

Hiϑ (x)p(x)dx

=

Z

M

′



′ 2
Hiϑ (x)p(x)dx + O e−(δ ) /ε ,

(93)



′ 2
and we mention that the exponential term O e−(δ ) /ε is negligible w.r.t to any polynomial asymptotic
expansion in ε, and is therefore omitted in subsequent analysis.

Continuing, we are interested in changing the integration variable from x to (z, β), as considered by
the following lemma.
′

′

Lemma 8 (Integration and volume form on M ). For any smooth h : M → R, we have that
Z

M

′

h(x)dx =

Z

z∈N

Z

2π

h(z β )V (z β )dzdβ,

(94)

β=0

′

where V (x) is associated with the volume form of M at x when integrating by (z, β), and is invariant to
β, i.e.

V (z β ) = V (z)

(95)

for all β ∈ [0, 2π).
The proof of Lemma 8 is provided in Section E.4. Hence, By Lemma 8, equations (91) – (93), and
the fact that p(x) = 1/Vol {M}, we have
1
lim Ci,N
(ϑ) =

N →∞

Z Z
N

2π
0

Hiϑ (z)V (z)p(z β )dzdβ =

2π
Vol {M}

Z

N

Hiϑ (z)V (z)dz.

(96)

Then, by substituting (89) into (96) we obtain
lim

N →∞

1
Ci,N
(ϑ)

2π
=
Vol {M}

Z Z
N

2π
0


exp − xϑi − z ϕ
31

2


/ε f (z ϕ )V (z)dzdϕ.

(97)

Eventually, if we change parametrization from (z, ϕ) back to x, via x = z ϕ , then by Lemma 8 we arrive
at
lim

N →∞

1
Ci,N
(ϑ)



Z
2
2π
ϑ
=
exp − xi − x /ε f (x)dx
Vol {M} M′


Z


′
2
2π
ϑ
−(δ )2 /ε
=
exp − xi − x /ε f (x)dx + O e
Vol {M} M


Z
2
2π
ϑ
exp − xi − x /ε f (x)dx,
=
Vol {M} M
2

ϑ
where we again
the fact that xi − x
 used
′ 2
resulting O e−(δ ) /ε term.

′

(98)

′

> (δ )2 for all x ∈
/ M (see Section (E.3)), and omitted the

In a similar fashion, if we consider the denominator of the second term in (86)
2
Ci,N
(ϑ)


N Z
1 X 2π
,
exp − xϑi − xϕ
j
N
0
j=1

2


/ε dϕ,

(99)

1 (ϑ) with f ≡ 1, we get that
and by repeating the calculations for Ci,N

lim

N →∞

2
Ci,N
(ϑ)

2π
=
Vol {M}

Z

M



2
ϑ
exp − xi − x /ε dx.

(100)

Lastly, if we substitute (98) and (100) into (86), we have that

lim

N →∞

4n
ε

o

L̃g (i, ϑ) =



4
f (xϑi ) −
ε

n
o

ϑ − x 2 /ε f (x)dx
exp
−
x
i
M

n
o
R
1
ϑ − x 2 /ε dx
exp
−
x
i
Vol{M} M

1
Vol{M}

R

= ∆M f (xϑi ) + O(ε),

(101)

where the last simplification (eq. (101)) is justified in [26].

E.2

The variance term

The variance error term in the convergence of the steerable graph Laplacian to the Laplace Beltrami
1 (ϑ), C 2 (ϑ) for finite N , and their
operator arises from the discrepancy between the values of Ci,N
i,N

limits (98), (100), respectively, when N → ∞. To prove the improved convergence rate of steerable graph
Laplacian, we follow the technique used in [26] which makes use of the Chernoff tail inequality. Such an

inequality provides a bound for the probability that a sum of random variables deviates from its mean by
a certain quantity.
Let us begin by defining
Gϑi (x)

,

Z

2π
0


exp − xϑi − xϕ

32

2


/ε dϕ.

(102)

Then, we are interested in evaluating the probabilities
p+ (N, α) =
p− (N, α) =
where we mention that the use of

PN

( PN

ϑ
j6=i Hi (xj )
P r PN
ϑ
j6=i Gi (xj )
( PN
ϑ
j6=i Hi (xj )
P r PN
ϑ
j6=i Gi (xj )

PN

ϑ
j6=i Hi (xj )/

)
 
E Hiϑ
−  ϑ > α ,
E Gi
)
 
E Hiϑ
−  ϑ  < −α ,
E Gi

ϑ
j6=i Gi (xj )

instead of

PN

(103)
(104)
PN

ϑ
j=1 Hi (xj )/

ϑ
j=1 Gi (xj )

(i.e.

without the diagonal) is justified at the end of the proof. We proceed by evaluating p+ (N, α), as p− (N, α)
can be obtained in a similar fashion. As was shown in [26], p+ (N, α) can be rewritten as

p+ (N, α) = P r


N
X


j6=i


 h i2 
,
Jiϑ (xj ) > α(N − 1) E Gϑi


(105)

where Jiϑ (xj ) are zero-mean i.i.d random variables (indexed by j), given by


h i h i
h i
h i
Jiϑ (xj ) , E Gϑi Hiϑ (xj ) − E Hiϑ Gϑi (xj ) + αE Gϑi E Gϑi − Gϑi (xj ) .

(106)

At this point, making use of Chernoff’s inequality gives
 4 

α2 (N − 1)2 E Gϑi
h
i
,
p+ (N, α) ≤ exp −
 2(N − 1)E J ϑ 2 + O(α) 
i



and it remains to evaluate the variance term E

h

Jiϑ

2 i

(107)

, which can be expressed as

    h i   
i  h i2  2 
h i h i h
2
2
2
ϑ ϑ
ϑ
ϑ
ϑ
ϑ
ϑ
+O(α). (108)
E Gϑi
−2E Gi E Hi E Hi Gi + E Hiϑ
E Hi
= E Gi
E Ji

 
 
Now, the integral expressions of E Hiϑ and E Gϑi (obtained in the previous section in equations (98)
and (100), respectively), admit asymptotic expansions via the following proposition.

˜
f be a smooth and compact d-dimensional
f → R be
Proposition 9. [26] Let M
submanifold, and let f˜ : M
f
a smooth function. Then, for any y ∈ M
(πε)

˜
−d/2

Z

n
o
i
εh
2
˜
exp −ky − xk2 /ε f˜(x)dx = f˜(y) +
E(y)f˜(y) + ∆M
ff (y) + O(ε ),
4
f
M

(109)

f at y.
where E(y) is a scalar function of the curvature of M
Then, it follows from Proposition 9 (see also [7, 5]) that



Z
h
i
2
2π
2π
ϑ
=
E
exp − xi − x /ε f (x)dx =
(πε)d/2 f (xϑi ) + O(ε) ,
Vol {M} M
Vol(M)


Z
h i
2
2π
2π
E Gϑi =
exp − xϑi − x /ε dx =
(πε)d/2 [1 + O(ε)] .
Vol {M} M
Vol(M)
h

Hiϑ

i

33

(110)
(111)

Thus, it remains to evaluate the second order moments E

h

Hiϑ

2 i

,E

h

Gϑi

2 i



, and E Hiϑ Gϑi , which is

the subject of the next lemma, whose proof is given in Section E.5.
r
PM
Pℓ m 2 2 2
Lemma 10. Let µ(x) ,
m=−M
ℓ=1 m xm,ℓ > 0 for every x ∈ N . Then


(πε)(d+1)/2 f 2 (xϑi )pN (xϑi )
+ O(ε) ,
=
E
2(d−1)/2
µ2 (xϑi )
  


2
(πε)(d+1)/2 pN (xϑi )
ϑ
=
E Gi
+ O(ε) ,
2(d−1)/2
µ2 (xϑi )

i (πε)(d+1)/2  f (xϑ )p (xϑ )
h
N
ϑ ϑ
i
i
E Hi Gi =
+ O(ε) ,
2(d−1)/2
µ2 (xϑi )


Hiϑ

2 

where pN (x) = 2πV (x)/ Vol {M} (and V (x) is from Lemma 8).
We mention that since we required (in Theorem 2) that

P

m6=0

Pℓ m

ℓ=1 |xm,ℓ |

(112)
(113)
(114)

2

> 0 for all x ∈ M (up to

a set of measure zero on M), then it is also the case that µ(x) > 0 for all x ∈ N with probability one.
Now, using Lemma 10, (110), (111), and substituting all quantities into (108), we get that
E



2
2π
(πε)(d+1)/2 f 2 (xϑi )pN (xϑi )
+
O(ε)
(πε)d (d−1)/2
Vol(M)
2
µ2 (xϑi )



2
(d+1)/2
f 2 (xϑi )pN (xϑi )
2π
d (πε)
−2
+ O(ε)
(πε)
Vol(M)
2(d−1)/2
µ2 (xϑi )

2


(πε)(d+1)/2 f 2 (xϑi )pN (xϑi )
2π
+
O(ε)
+ O(α)
(πε)d (d−1)/2
+
Vol(M)
2
µ2 (xϑi )

2
2π
(πε)(d+1)/2
=
(πε)d (d−1)/2 · O(ε) + O(α) = O(ε3d/2+3/2 ) + O(α).
Vol(M)
2

   
2
Jiϑ
=

Additionally, from (111) we have
E
and thus

h

Gϑi


p+ (N, α) ≤ exp −

i4

= O(ε2d ),


α2
.
O(ε−d/2+3/2 /N ) + O(α)

(115)

(116)

(117)

√
Henceforth, by taking α = O(ε−d/4+3/4 / N ) we can make p+ (N, α) arbitrarily small with exponential
decay. Additionally, we mention that p− (N, α) leads to the same asymptotic expression. Therefore, it
follows that with high probability
PN

ϑ
j6=i Hi (xj )
PN
ϑ
j6=i Gi (xj )

 
E Hiϑ
1
ε−d/4+3/4
−  ϑ  = |α| = O(
) = O( 1/2 (d−1)/4−1/2 ).
1/2
N
N ε
E Gi

34

(118)

Continuing, we can write (using (101))
!
PN
PN
 !
 
ϑ
ϑ
E Hiϑ
E Hiϑ
4
4
4
j6=i Hi (xj )
j6=i Hi (xj )
ϑ
ϑ
f (xi ) − PN
f (xi ) −  ϑ 
−  ϑ =
−
PN
ϑ (x )
ϑ (x )
ε
ε
ε
E
G
E Gi
G
G
j
j
i
j6=i i
j6=i i
!
PN
ϑ


4
1
j6=i Hi (xj )
f (xϑi ) − PN
− ∆M f (xϑi ) + O(ε) = O( 1/2 (d−1)/4+1/2 ),
=
ϑ
ε
N ε
j6=i Gi (xj )

(119)

which gives us that
4
ε

f (xϑi ) −

PN

ϑ
j6=i Hi (xj )
PN
ϑ
j6=i Gi (xj )

!

= ∆M f (xϑi ) + O(

1
N 1/2 ε(d−1)/4−1/2

) + O(ε).

(120)

The last step of the proof is to justify that removing the diagonal of the steerable affinity operator W
(i.e. computing all sums with j 6= i) does not change the convergence rate. Indeed, this is the case since
P

PN

N
ϑ
ϑ
j6=i Hi (xj )
j=1 Hi (xj )
− PN
=
PN
ϑ
ϑ
j=1 Gi (xj )
j6=i Gi (xj )
PN
PN
ϑ
ϑ
j=1 Hi (xj )
j6=i Hi (xj )
= PN
− PN
ϑ
ϑ
j=1 Gi (xj )
j=1 Gi (xj )

=O

Gϑ (xi )
PN i ϑ
j=1 Gi (xj )

!

PN

ϑ
j=1 Hi (xj )
PN
ϑ
j=1 Gi (xj )

,

−

Gϑi (xi )

1 + PN

PN

ϑ
j6=i Hi (xj )
PN
ϑ
j=1 Gi (xj )

ϑ
j6=i Gi (xj )

!

PN

ϑ
j=1 Gi (xj )
ϑ
j6=i Gi (xj )

PN

Hiϑ (xi )

= PN

Gϑi (xi )

ϑ
j=1 Gi (xj )

− PN

ϑ
j=1 Gi (xj )

PN

ϑ
j6=i Hi (xj )
PN
ϑ
j6=i Gi (xj )

(121)

where we used the fact that since f is smooth then it is also bounded on M, satisfying |f (x)| ≤ c, and
hence Hiϑ (xi ) ≤ c Gϑi (xi ) . Therefore, by using Proposition 9 (specifically, retaining the zero-order
element in the asymptotic expansion in (109)) it follows that
Gϑ (xi )
=
PN i ϑ
j=1 Gi (xj )

1
N

1 ϑ
1
ε1/2 /N
N Gi (xi )
) = O( (d−1)/2 ),
=
O(
P
N
d/2
ϑ
ε
Nε
j=1 Gi (xj )

(122)

which is negligible compared to variance error term of (120). Overall, we get that
4n o
4
L̃g (i, ϑ) =
ε
ε

PN

ϑ
j=1 Hi (xj )
f (xϑi ) − PN
ϑ
j=1 Gi (xj )

!

= ∆M f (xϑi ) + O(

1
N 1/2 ε(d−1)/4−1/2

) + O(ε),

(123)

which concludes the proof.

E.3

Construction of a rotationally-invariant parametrization

We construct a parametrization (z, β) of all points x in a certain neighbourhood of xϑi . This parametrization has favorable properties for our purposes, and is specific for every index i and rotation angle ϑ. The
parametrization is defined by the mapping x 7→ (z, β), given by
z(x) = R(x, α̂(x)),

β(x) = −α̂(x),

35

α̂(x) = arg min xα − xϑi
α∈[0,2π)

2
2

.

(124)

That is, z is the rotation of x (by (6)) which is closest to xϑi , and therefore this parametrization satisfies
x = R(z, β) = z β .

(125)

Note that this mapping is invariant to the intrinsic rotation, that is, different values of x which differ only
by a rotation will be mapped to the same z. Therefore, the parametrization (z, β) can be perceived as
a form of a polar parametrization, where coordinates are parametrized by a radius (the equivalent of z)
and a rotation angle (the equivalent of β).
Now, a solution to
xα − xϑi

2

−→ min

(126)

2 α∈[0,2π)

must exist, but may not be unique for all x ∈ M. We start by showing that it is guaranteed to be

unique in a sufficiently small neighbourhood of xϑi . To this end, note that due to our requirements in
P
P
2
Theorem (2), we have with probability one that |m|>0 ℓ xi,(m,ℓ) > 0 (i.e. the image corresponding
to xϑ is not radially symmetric), and since M is smooth, there exists a neighbourhood of xϑi for which
P i P
2
> 0. Let us consider a ball of radius δ > 0 around xϑi , denoted by Bδ (xϑi ), such
|m|>0
ℓ x(m,ℓ)
P
P
2
> 0 for all points x ∈ Bδ (xϑi ) ∩ M. Clearly, fixing some x ∈ Bδ (xϑi ) ∩ M,
that
|m|>0
ℓ x(m,ℓ) n
o
2
we have that minα∈[0,2π) xα − xϑi 2 ≤ δ2 , and moreover, the curve {xα }2π
α=0 has bounded curvature
′

and does not self intersect, which means that there exists a sufficiently small δ ≤ δ such that any

solution to xα − xϑi

2
−→
2 α∈[0,2π)

′

min must satisfy xα ∈ Bδ′ (xϑi ) ∩ M. Moreover, we choose δ such that

the solution is unique, which is justified by the fact that the problem xα − xϑi

2
−→
2 α∈[0,2π)

min restricted

′

to xα ∈ Bδ′ (xϑi ) ∩ M, must be convex for a sufficiently small δ (due to the smoothness of the curve
ϑ
{xα }2π
α=0 ). To conclude, for points x in a sufficiently small neighbourhood of xi , the mapping x 7→ (z, β)

given by (124) is unique.

Next, it is of interest to characterize the resulting set of feasible points of z, and we proceed by showing
that

o
n
N , z(x) : x ∈ Bδ′ (xϑi ) ∩ M

(127)

is a smooth and compact (d − 1)-dimensional submanifold. Now, from (126) it follows that α = α̂(x)

must be a solution of

Re



∂xα α
, x − xϑi
∂α



= 0,

(128)

which can be written explicitly via (6) as
Re



∂xα α
, x − xϑi
∂α



=

ℓm
M
X
X

m=−M ℓ=1

=

ℓm
M
X
X

m=−M ℓ=1

=−

o
n

Re ım · x∗(m,ℓ) e−ımα · x(m,ℓ) eımα − xϑi,(m,ℓ)

n
Re ım · x(m,ℓ)

ℓm
M
X
X

m=−M ℓ=1

2

− ım · x∗(m,ℓ) e−ımα · xϑi,(m,ℓ)

o
n
Re ımx∗(m,ℓ) xϑi,(m,ℓ) · e−ımα .

36

o

(129)

Next, since at α = α̂(x) we have that z = xα , then by (124) each z satisfies
Re

(*

∂z β
∂β

β=0

,z −

xϑi

+)

= 0,

(130)

which can be written explicitly via (129) as
ℓm
M
X
X

m=−M ℓ=1

n
o
∗
Re ım · z(m,ℓ)
· xϑi,(m,ℓ) = 0,

(131)

where z(m,ℓ) denotes the (m, ℓ)‘th coordinate of z. Essentially, equation (131) defines a linear subspace
of the ambient space, such that all feasible points for z (z = xα for α satisfying (126)) must lie in the
intersection between Bδ′ (xϑi ) ∩ M and this subspace. In particular, the submanifold N can be explicitly
defined through

N =

(

ℓm
M
X
X

x:

m=−M ℓ=1

o

n

)

Re ım · x∗(m,ℓ) · xϑi,(m,ℓ) = 0, x ∈ Bδ′ (xϑi ) ∩ M ,

(132)

which is a smooth and compact submanifold due to the smoothness and compactness of M, and is

of intrinsic dimension (d − 1) due to the additional linear constraint (note that this constraint is not
P
P
2
degenerate since |m|>0 ℓ x(m,ℓ) > 0 for x ∈ Bδ (xϑi ) ∩ M).
Lastly, we make the observation that all rotations of any point in the neighbourhood Bδ′ (xϑi ) ∩ M, i.e.

xϕ for all ϕ ∈ [0, 2π) and any x ∈ Bδ′ (xϑi ) ∩ M, share the same solution (the same xα for α from (126)) as

the point x. This allows us to extend the neighbourhood in which our parametrization is valid by taking
′
all rotations of all points in this neighbourhood, and in particular, we conclude that for all x ∈ M , where
o
n
′
M , xϕ : x ∈ Bδ′ (xϑi ) ∩ M, ϕ ∈ [0, 2π) ,

(133)
′

the parametrization x 7→ (z, β), where x = z β , is unique. Additionally, it is evident that M ⊂ M is also
a smooth and compact d-dimensional submanifold.

E.4

Proof of Lemma 8



Proof. Let N be parametrized locally by u = u1 , . . . , u(d−1) ∈ Rd−1 around a point z0 . That is, every

coordinate zm,ℓ of the manifold N is expressed as a function of u in the vicinity of the point z0 . Then,

M can be parametrized locally around x0 = z0β0 by [u, β] ∈ Rd , where β ∈ [0, 2π) is the rotation angle
′

from our rotationally-invariant parametrization (see section E.3), in the sense that every coordinate xm,ℓ
′

of M can be expressed as a function of [u, β] in the neighbourhood of the point x0 . Hence, the integral
′

of a function h(x) over M can be expressed through the parametrization x 7→ (z, β) by
Z

M

′

h(x)dx =

where
dV (x) =

Z

z∈N

q

Z

2π

h(z β )dV (z β ),


det gM′ (x) du1 . . . dud−1 dβ
37

(134)

β=0

(135)

′

is the volume form at the point x, gM′ (x) is the metric tensor on M , given by pull-back as
o
n
∗
′ (x)
,
gM′ (x) = Re JM
′ (x)J
M

(136)

and JM′ (x) is the Jacobian matrix


h
i
JM′ (x) = Ju Jβ ,

∂x(−M,1)
∂u1

∂x(−M,1)
∂u(d−1)

...



..

.

 ∂x(−M,ℓm )
 ∂u1


..
Ju = 
.

 ∂x(M,1)
 ∂u1


..

.


...

..
.

...

∂x(−M,ℓm )
∂u(d−1)

..
.
∂x(M,1)
∂u(d−1)

...

..
.

...

∂x(M,ℓm )
∂u1

∂x(M,ℓm )
∂u(d−1)

...











,










∂x(−M,1)
∂β





..



.


 ∂x
 (−M,ℓm ) 


∂β


.
.

.
Jβ = 
.

 ∂x(M,1) 


∂β




.
..





(137)

∂x(M,ℓm )
∂β

Note that since x = z β we have xm,ℓ = zm,ℓ eımβ (from (6)), and thus
∂z(m,ℓ) ımβ
∂x(m,ℓ)
=
·e ,
∂ui
∂ui

i = 1, . . . , d − 1;

∂x(m,ℓ)
= ım · zm,ℓ · eımβ .
∂β

(138)

Therefore, it is evident that the metric tensor
o
n
∗
′ (x)
=
gM′ (x) = Re JM
′ (x)J
M

"

Ren
{Ju∗ (x)Ju (x)}
{Ju∗ (x)Jβ (x)}
o Ren
o
Re Jβ∗ (x)Ju (x)
Re Jβ∗ (x)Jβ (x)

#

(139)

does not depend on β, i.e. gM′ (x) = gM′ (z), since the term eımβ cancels-out in all entries of gM′ (x).
Consequently, we have that
dV (x) =

q


det gM′ (z) du1 . . . dud−1 dβ = V (z)dzdβ,

V (z) =

q


det gM′ (z) ,

where we denoted

and it follows that

Z

M

E.5

′

h(x)dx =

Z

z∈N

dz = du1 . . . dud−1 ,

Z

2π

h(z β )V (z)dzdβ.

(140)

(141)

(142)

β=0

Proof of Lemma 10

h
h
2 i
2 i
and
We put our focus on evaluating the term E Hiϑ , as the other second-order terms E Gϑi
 ϑ ϑ
E Hi Gi can be obtained in a similar fashion. Thus, we are interested in evaluating the term
   Z
2
=
E Hiϑ

M


2
Hiϑ (x) p(x)dx,

38

(143)

recalling that
Hiϑ (x) =

Z

2π
0


exp − xϑi − xϕ

2


/ε f (xϕ )dϕ.

(144)

Using the construction of our rotationally-invariant parametrization (see Section E.3), and in particular
′

the submanifold M , we can write
E



2 

Hiϑ

=

Z

M′


2


′ 2
Hiϑ (x) p(x)dx + O e−(δ ) /ε ,

(145)



′ 2
′
for some constant δ > 0. Next, as the O e−(δ ) /ε term is negligible w.r.t to any polynomial asymptotic

expansion in ε, we omit it. Then, by Lemma 8 (i.e. change of integration variables x 7→ (z, β)) and (91),
we have

E



Hiϑ

2 

=

Z Z

2π



2

Hiϑ (z β )

2π
p(z )V (z)dzdβ =
Vol {M}
β

ZN 0
2
Hiϑ (z) pN (z)dz,
=

Z 
N

2
Hiϑ (z) V (z)dz
(146)

N

where we defined
pN (z) =

2π
V (z).
Vol {M}

(147)

We start by deriving an asymptotic expression for Hiϑ (z). Let us write
xϑi − z ϕ

2


2
xϑi − z + (z − z ϕ )
Eo
nD
2
+ kz − z ϕ k2 ,
= xϑi − z + 2 Re xϑi − z, z − z ϕ
=



and denote
δiϑ (z, x) , 2 Re
Therefore, we have that

Hiϑ (z) = exp − xϑi − z

2

/ε

Z

2π

0

nD

xϑi − z, z − x

Eo

.

n
o
n
o
exp −kz − z ϕ k2 /ε exp −δiϑ (z, z ϕ )/ε f (z ϕ )dϕ.

(148)

(149)

(150)

Next, we use Taylor expansion to write
exp

n

−δiϑ (z, z ϕ )/ε

o

δϑ (z, z ϕ )
=1− i
+O
ε

39



2 !
δiϑ (z, z ϕ )
,
ε2

(151)

which gives us that
Hiϑ (z)


= exp − xϑi − z

2

/ε

 Z

0

2π

Z

n
o
exp −kz − z ϕ k2 /ε f (z ϕ )dϕ

n
o
ϕ 2
exp −kz − z k /ε δiϑ (z, z ϕ )f (z ϕ )dϕ
0 Z 2π

n
oh
i2
1
ϕ 2
ϑ
ϕ
ϕ
+O 2
exp −kz − z k /ε δi (z, z ) |f (z )| dϕ .
ε 0

1
−
ε

2π

(152)

In what follows, we evaluate the terms in the square brackets of (152) one by one, where we mention
that Proposition 9 is the main workhorse for obtaining the asymptotic expansions of the integrals taking
part in our analysis. To this end, let us first define the set of points Cz = {z ϕ }2π
ϕ=0 , which is a smooth
curve in CD . We then change the integration parameter in (152) from the angle ϕ to the variable x = z ϕ

(which is equivalent to parametrizing by arc-length), and if we recall that f (x) is a smooth function, then
by Proposition 9 we get that
Z

0

where µ(z) =

r

2π

Z
n
o
1
exp −kz − xk2 /ε f (x)dx
µ(z) x∈Cz
√
πε
=
[f (z) + O(ε)] ,
µ(z)

n
o
exp −kz − z ϕ k2 /ε f (z ϕ )dϕ =

(153)

h n
oi
det Re Jβ∗ (z)Jβ (z)
(Jβ is defined in (137)) is associated with the change of the

integration variable, and is given explicitly by

v
u M ℓm
u X X
m2 |zm,ℓ |2 .
µ(z) = t

(154)

m=−M ℓ=1

Next, we evaluate the second term in the square brackets of (152). Since δiϑ (z, z ϕ ) is a smooth function
in ϕ, and using the previous change of variable x = z ϕ together with Proposition 9, we have
Z
n
o
n
o
1
exp −kz − z ϕ k2 /ε δiϑ (z, z ϕ )f (z ϕ )dϕ =
exp −kz − xk2 /ε δiϑ (z, x)f (x)dx
εµ(z) x∈Cz
0




√
n
o
πε
ε
ϑ
2
ϑ
ϑ
=
+ O(ε )
δ (z, z)f (z) +
E(y)δi (z, z)f (z) + ∆Cz δi (z, x)f (x)
εµ(z) i
4
x=z


√
n
o
πε ε
ϑ
2
+ O(ε )
(155)
∆C δ (z, x)f (x)
=
εµ(z) 4 z i
x=z

1
ε

Z

2π

since it is clear from (149) that δiϑ (z, z) = 0. Moreover, we have (see Lemma 3.3 in [24]) that
n
o
D
E
∆Cz δiϑ (z, x)f (x) |x=z = ∆Cz f (x)|x=z · δiϑ (z, z) − 2 ∇Cz δiϑ (z, x)|x=z , ∇Cz f (z) + f (z) · ∆Cz δiϑ (z, x)|x=z
= f (z) · ∆Cz δiϑ (z, x)|x=z ,

(156)

40

where we have used the fact that δiϑ (z, z) = 0, and moreover, that (using (149))
∇Cz δiϑ (z, x)|x=z

Eo
nD
ϑ
= −2 Re
= −2 Re xi − z, ∇Cz x|x=z

as xϑi − z is perpendicular to

∂z β
∂β

(*

xϑi

1 ∂z β
− z,
µ(z) ∂β

β=0

+)

=0

(157)

by our rotationally-invariant parametrization (equation (130) in
β=0

Section E.3). Therefore, we are left with
1
ε

Z

2π
0

√
√
n
o

πε 
πε
[q(z) + O(ε)] ,
exp −kz − z ϕ k2 /ε δiϑ (z, z ϕ )f (z ϕ )dϕ =
εq(z) + O(ε2 ) =
εµ(z)
µ(z)

(158)

where we defined the function q(z) as
q(z) =

Eo
nD
f (z)
f (z)
,
∆Cz δiϑ (z, x)|x=z = −
Re xϑi − z, ∆Cz x|x=z
4
2

(159)

with the second equality following from (149). The notation ∆Cz x|x=z denotes the Laplacian of each
coordinate of x, taken w.r.t to the curve Cz , and sampled at the point z. Consequently, q(z) is a smooth
function satisfying

q(z) = O( xϑi − z ),

q(xϑi ) = 0,

(160)

where we applied the Cauchy-Schwarz inequality to (159), together with the fact that k∆Cz x|x=z k is

bounded (since Cz is smooth).

Now, as for the last term in the square brackets of (152), we first mention that since f (x) is smooth,

it is bounded (on a compact domain) and therefore
O



1
ε2

Z

2π
0


n
oh
i2
ϕ 2
ϑ
ϕ
ϕ
exp −kz − z k /ε δi (z, z ) |f (z )| dϕ =
 Z 2π
n
oh
i2 
1
ϕ 2
ϑ
ϕ
exp −kz − z k /ε δi (z, z ) dϕ . (161)
O
ε2 0

Moreover, if we expand
ϕ
zm,ℓ

ımϕ

= zm,ℓ e

= zm,ℓ +

then we have that
zϕ − z =

∂z ϕ
∂ϕ

ϕ
∂zm,ℓ

ϕ=0

∂ϕ

ϕ=0

· ϕ + O(ϕ2 ),

(162)

· ϕ + O(ϕ2 ),

(163)

and it is evident that
Eo
nD
= 2 Re
δiϑ (z, z ϕ ) = −2 Re xϑi − z, z − z ϕ
2

= O(ϕ

xϑi

(*

∂z ϕ
xϑi − z,
∂ϕ

ϕ=0

+)

· ϕ + O(ϕ2 xϑi − z )

− z ),

(164)

where we used Cauchy-Schwarz inequality and again the fact that xϑi − z is perpendicular to

41

∂z ϕ
∂ϕ

due
ϕ=0

to our rotationally-invariant parametrization (eq. (130)). Eventually, we obtain that
h
i2
δiϑ (z, z ϕ ) = O(ϕ4 xϑi − z

2

).

(165)

Continuing, from (163) it is clear that
ϕ = O(kz − z ϕ k),
and therefore

h

i2
δiϑ (z, z ϕ ) = O(kz − z ϕ k4 xϑi − z

(166)
2

).

(167)

When plugging (167) back into (161) and changing the integration parameter from ϕ to x = z ϕ , we arrive
at
O



1
ε2

=O

Z


Z
n
oh
i2 
n
oh
i2 
1
2
ϕ 2
ϑ
ϕ
ϑ
exp −kz − z k /ε δi (z, z ) dϕ = O
exp −kz − xk /ε δi (z, x) dx
ε2 µ(z) x∈Cz
0
!
!
!
2
2
2 Z
n
o
xϑi − z
xϑi − z √
xϑi − z
2
4
2√
· ε πε = O
exp −kz − xk /ε kz − xk dx = O
πε ,
ε2 µ(z)
ε2 µ(z)
µ(z)
x∈Cz
2π

(168)

where we used the asymptotic expansion in Proposition 9 together with the fact that the function kz − xk4

(in z) and its Laplacian vanish at z = x (leaving only the O(ε2 ) term in the asymptotic expansion of
Proposition 9).
Altogether, when plugging (153), (158), and (168) into (152), we get

Hiϑ (z)

=

n
exp − xϑi − z
µ(z)

2

o
/ε √ 
πε f (z) + q(z) + O( xϑi − z

2


) + O(ε) ,

(169)

where q(z) was defined in (159). Therefore, we have
n

2 exp −2 xϑi − z
Hiϑ (z) =
µ2 (z)
n
exp −2 xϑi − z
=
µ2 (z)

2

2

/ε
/ε

o

o


πε f (z) + q(z) + O( xϑi − z


2

πε f (z) + 2f (z)q(z) + O(

2

2
) + O(ε)

xϑi

−z

2


) + O(ε) ,

(170)

where we used the fact (from (160)) that q(z) = O( xϑi − z ), and retained only the asymptotically
dominant factors inside the square brackets.
h
2 i
by plugging (170) into (143). We have
We are now ready to evaluate E Hiϑ
   Z 
2
2
Hiϑ (z) pN (z)dz
=
E Hiϑ
nN
o
Z exp −2 xϑ − z 2 /ε 
i
= πε
f 2 (z) + 2f (z)q(z) + O( xϑi − z
2 (z)
µ
N

2


) + O(ε) pN (z)dz.

(171)

Before we proceed with the asymptotic analysis, we mention that if µ(z) > 0 for all z ∈ N , then 1/µ2 (z)
42

is a smooth function. Additionally, the smoothness of pN (z) = 2πV (z)/ Vol {M} is guaranteed by the
smoothness of M and N (see the definition of V (x) in Section E.4). Then, we expand the square brackets
in (171) and evaluate (asymptotically) all resulting integrals by applying Proposition 9. We get
Z

N

Z

N



exp −2

xϑi

2

−z


exp −2 xϑi − z

2




 2 ϑ
ϑ
f 2 (z)pN (z)
(d−1)/2 f (xi )pN (xi )
/ε
+ O(ε) ,
dz = (πε/2)
µ2 (z)
µ2 (xϑi )



(172)



ϑ
ϑ
ϑ
f (z)q(z)pN (z)
(d−1)/2 f (xi )q(xi )pN (xi )
/ε
dz = (πε/2)
+ O(ε)
µ2 (z)
µ2 (xϑi )
= (πε/2)(d−1)/2 · O(ε),

(173)

since q(xϑi ) = 0 (see (160)), and
n
exp −2 xϑi − z

Z

2

µ2 (z)

N

where we used the fact that xϑi − z

2

/ε

o

O( xϑi − z

2

)pN (z)dz = (πε/2)(d−1)/2 · O(ε),

(174)

is smooth and vanishes at z = xϑi (using Proposition 9, we are left

only with the O(ε) term in the expansion).
Finally, by substituting all of the above asymptotic integral expansions into (171), it follows that



  
 2 ϑ
ϑ
2
(πε)(d+1)/2 f 2 (xϑi )pN (xϑi )
(d−1)/2 f (xi )pN (xi )
ϑ
= πε(πε/2)
+ O(ε) =
+ O(ε) .
E Hi
2(d−1)/2
µ2 (xϑi )
µ2 (xϑi )
h

(175)

2 i



and E Hiϑ Gϑi can be obtained in exactly the same manner, and we omit the derivation
h
2 i
it is sufficient to take f ≡ 1 throughout the
for the sake of brevity (note that to compute E Gϑi
Then, E

Gϑi

derivation).

Appendix F

Proof of Theorem 3

We mention that this proof follows very closely the proof of Theorem 1.
Proof. First, as


D

−1

W f (i, ϑ) =

N Z
X
j=1

2π

(Wi,j (ϑ, ϕ)/Di,i ) fj (ϕ)dϕ,

(176)

0

it is evident that D −1 W is also LRI (since D −1 W merely alters Wi,j by constant factors independent of
ϑ), and hence L = I − D −1 W is of the form A + G as required by Proposition 7. Therefore, we can obtain

a sequence of eigenfunctions and eigenvalues of L by diagonalizing the matrices S̃m = I − D −1 Ŵ (m) for
every m ∈ Z. However, it is important to mention that in contrast to Sm = D − Ŵ (m) , the matrix
S̃m = I − D −1 Ŵ (m) is not Hermitian. Nonetheless, if we make the observation that S̃m is similar to the
Hermitian matrix Sm = I − D −1/2 Ŵ (m) D −1/2 by
′

′

D 1/2 S̃m D −1/2 = Sm ,
43

(177)

then it follows that S̃m can be diagonalized with a set of eigenvectors complete in CN and the eigenvalues
′

of Sm , which are real-valued. Then, as it follows from Theorem 1 that the eigenvalues of Sm are non′

negative, we have that the eigenvalues of Sm = D −1/2 Sm D −1/2 must be also non-negative (since surely
v ∗ D −1/2 Sm D −1/2 v ≥ 0 for any v ∈ CN ), which lastly implies that the eigenvalues of S̃m are non-negative.

As to the fact that {Φm,k }m,k are complete in H, the same arguments as in the proof of Theorem 1

hold.

References
[1] Yonathan Aflalo, Haim Brezis, and Ron Kimmel. On the optimality of shape and data representation
in the spectral domain. SIAM Journal on Imaging Sciences, 8(2):1141–1160, 2015.
[2] Yariv Aizenbud and Amir Averbuch. Matrix decompositions using sub-gaussian random matrices.
arXiv preprint arXiv:1602.03360, 2016.
[3] Jinho Baik, Gérard Ben Arous, Sandrine Péché, et al. Phase transition of the largest eigenvalue for
nonnull complex sample covariance matrices. The Annals of Probability, 33(5):1643–1697, 2005.
[4] Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps for dimensionality reduction and data
representation. Neural computation, 15(6):1373–1396, 2003.
[5] Mikhail Belkin and Partha Niyogi. Towards a theoretical foundation for laplacian-based manifold
methods. Journal of Computer and System Sciences, 74(8):1289–1308, 2008.
[6] Yifan Cheng, Nikolaus Grigorieff, Pawel A Penczek, and Thomas Walz. A primer to single-particle
cryo-electron microscopy. Cell, 161(3):438–449, 2015.
[7] Ronald R Coifman and Stéphane Lafon. Diffusion maps. Applied and computational harmonic
analysis, 21(1):5–30, 2006.
[8] Ronald R Coifman, Yoel Shkolnisky, Fred J Sigworth, and Amit Singer. Graph laplacian tomography
from unknown random projections. IEEE Transactions on Image Processing, 17(10):1891–1899, 2008.
[9] Mathieu Desbrun, Mark Meyer, Peter Schröder, and Alan H Barr. Implicit fairing of irregular meshes
using diffusion and curvature flow. In Proceedings of the 26th annual conference on Computer graphics
and interactive techniques, pages 317–324. ACM Press/Addison-Wesley Publishing Co., 1999.
[10] Michael Eisenstein. The field that came in from the cold. Nature methods, 13(1):19, 2016.
[11] Noureddine El Karoui et al. On information plus noise kernel random matrices. The Annals of
Statistics, 38(5):3191–3216, 2010.
[12] Noureddine El Karoui, Hau-tieng Wu, et al. Graph connection laplacian methods can be made robust
to noise. The Annals of Statistics, 44(1):346–372, 2016.
[13] Joachim Frank. Three-Dimensional Electron Microscopy of Macromolecular Assemblies: Visualization of Biological Molecules in Their Native State. Oxford, 2006.
44

[14] William T. Freeman and Edward H Adelson. The design and use of steerable filters. IEEE Transactions on Pattern Analysis & Machine Intelligence, 13(9):891–906, 1991.
[15] Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM review, 53(2):217–
288, 2011.
[16] Matthias Hein and Markus Maier. Manifold denoising. In NIPS, volume 19, pages 561–568, 2006.
[17] Boris Landa and Yoel Shkolnisky. Approximation scheme for essentially bandlimited and spaceconcentrated functions on a disk. Applied and Computational Harmonic Analysis, 2016.
[18] Boris Landa and Yoel Shkolnisky. Steerable principal components for space-frequency localized
images. SIAM Journal on Imaging Sciences, 10(2):508–534, 2017.
[19] Xianming Liu, Deming Zhai, Debin Zhao, Guangtao Zhai, and Wen Gao. Progressive image denoising
through hybrid graph laplacian regularization: a unified framework. IEEE Transactions on image
processing, 23(4):1491–1503, 2014.
[20] François G Meyer and Xilin Shen. Perturbation of the eigenvectors of the graph laplacian: Application
to image denoising. Applied and Computational Harmonic Analysis, 36(2):326–334, 2014.
[21] Stanley Osher, Zuoqiang Shi, and Wei Zhu. Low dimensional manifold model for image processing.
SIAM Journal on Imaging Sciences, 10(4):1669–1690, 2017.
[22] Pietro Perona. Deformable kernels for early vision. Pattern Analysis and Machine Intelligence, IEEE
Transactions on, 17(5):488–499, 1995.
[23] Gabriel Peyré. Manifold models for signals and images. Computer Vision and Image Understanding,
113(2):249–260, 2009.
[24] Steven Rosenberg. The Laplacian on a Riemannian manifold: an introduction to analysis on manifolds. Number 31. Cambridge University Press, 1997.
[25] Sam T Roweis and Lawrence K Saul. Nonlinear dimensionality reduction by locally linear embedding.
science, 290(5500):2323–2326, 2000.
[26] Amit Singer. From graph to manifold laplacian: The convergence rate. Applied and Computational
Harmonic Analysis, 21(1):128–134, 2006.
[27] Amit Singer and H-T Wu. Vector diffusion maps and the connection laplacian. Communications on
pure and applied mathematics, 65(8):1067–1144, 2012.
[28] Amit Singer and H-T Wu. Two-dimensional tomography from noisy projections taken at unknown
random directions. SIAM journal on imaging sciences, 6(1):136–175, 2013.
[29] David Slepian. Prolate spheroidal wave functions, Fourier analysis and uncertainty - IV: extensions to many dimensions; generalized prolate spheroidal functions. Bell System Technical Journal,
43(6):3009–3057, 1964.
45

[30] Gabriel Taubin. A signal processing approach to fair surface design. In Proceedings of the 22nd
annual conference on Computer graphics and interactive techniques, pages 351–358. ACM, 1995.
[31] Joshua B Tenenbaum, Vin De Silva, and John C Langford. A global geometric framework for
nonlinear dimensionality reduction. science, 290(5500):2319–2323, 2000.
[32] Rebecca F Thompson, Matt Walker, C Alistair Siebert, Stephen P Muench, and Neil A Ranson. An
introduction to sample preparation and imaging by cryo-electron microscopy for structural biology.
Methods, 100:3–15, 2016.
[33] Bruno Vallet and Bruno Lévy. Spectral geometry processing with manifold harmonics. In Computer
Graphics Forum, volume 27, pages 251–260. Wiley Online Library, 2008.
[34] Zhizhen Zhao, Yoel Shkolnisky, and Amit Singer. Fast steerable principal component analysis. IEEE
Transactions on Computational Imaging, 2(1):1–12, 2016.
[35] Zhizhen Zhao and Amit Singer. Fourier–Bessel rotational invariant eigenimages. JOSA A, 30(5):871–
877, 2013.

46

