Fast expansion into harmonics on the disk: a steerable basis with fast radial convolutions

1. Introduction

1.1. Motivation

Decomposing a function into its Fourier series can be viewed as representing a function in the eigenfunctions of the Laplacian on the torus $T := [0, 2 π]$ where $0$ and $2 π$ are identified. Indeed,

- Δ e^{ı k x} = k^{2} e^{ı k x} .

The eigenfunctions of the Laplacian (harmonics) on the disk ${x \in R^{2} : | x | < 1}$ that satisfy the Dirichlet boundary conditions can be written in polar coordinates $(r, θ) \in [0, 1) \times [0, 2 π)$ as

ψ_{n k} (r, θ) = c_{n k} J_{n} (λ_{n k} r) e^{ı n θ},

(1)

where $c_{n k}$ is a normalization constant, $J_{n}$ is the $n$ -th order Bessel function of the first kind, and $λ_{n k}$ is the $k$ -th smallest positive root of $J_{n}$ . The indices run over $(n, k) \in Z \times Z_{> 0}$ . The functions $ψ_{n k}$ satisfy

- Δ ψ_{n k} = λ_{n k}^{2} ψ_{n k} .

(2)

In this paper, we present a fast and accurate transform of digitized $L \times L$ images into this eigenfunction basis, often referred to as the Fourier-Bessel basis. For computational purposes, this basis is convenient for a number of reasons:

[label=()]
Orthonormal: these eigenfunctions are an orthonormal basis for square integrable functions on the disk.
Ordered by frequency:
the basis functions are ordered by eigenvalues, which can be interpreted as frequencies due to the connection with the Laplacian and Fourier series described above; therefore, low-pass filtering can be performed by retaining basis coefficients up to a given threshold.
Steerable: functions expanded in the basis can be rotated by applying a diagonal transform corresponding to phase modulation of the coefficients.
Fast radial convolutions: we show that the convolution with radial functions can be computed by applying a diagonal transform to the coefficients.

Figure 1. Illustration of method for $L \times L$ images with $L = 256$ . Original image (a), low-pass filter of the original image using a decreasing number of basis functions (b–d), radial function (e), convolution of original image with radial function (f).

Our method involves $O (L^{2} log L)$ operations and has precise accuracy guarantees. Python code that implements our method is publicly available online¹¹1An implementation is available at https://github.com/nmarshallf/fle_2d.. To the best of our knowledge, existing methods for computing the expansion coefficients in a steerable basis either require $O (L^{3})$ operations [19, 20] or lack accuracy guarantees [34, 35].

Steerable bases have been utilized in numerous image processing problems including image alignment [25], image classification [35] and image denoising [34]

, including applications to machine learning

[5, 7, 33] and data-driven science, such as applications to cryo-electron microscopy (cryo-EM) [6, 23]

, and computer vision

[24], among other areas.

There are many possible choices of steerable bases, for instance prolate spheroidal wave functions which are bandlimited functions optimally supported on the disk [19, 20, 28], or Zernike polynomials which are widely used in optics [31]. The harmonics on the disk (which satisfy Dirichlet boundary conditions) [34, 35] are one natural choice due to their orthogonality, ordering by frequency, and fast radial convolution.

We illustrate the frequency ordering property of the Laplacian eigenbasis by performing a low-pass filter by projecting onto the span of eigenfunctions whose eigenvalues are below a sequence of bandlimits that decrease the number of basis functions successively by factors of four, starting from $39593$ coefficients. Since the basis is orthonormal this is equivalent to setting coefficients above the bandlimit equal to zero; see Fig. 1 (a–d). Further, we demonstrate the radial convolution property by illustrating the convolution with a point spread function, which is a function used in computational microscopy [32]; see Fig. 1 (e–f). The image used for this example is a tomographic projection of a 3-D density map representing a bio-molecule (E. coli 70S ribosome) [26].

1.2. Notation

Suppose that $f$ is an $L \times L$ image representing a function $~ f : [- 1, 1]^{2} \to R$ that is supported on the unit disk ${x \in R^{2} : | x | < 1}$ . More precisely, we assume that

f (j_{1}, j_{2}) = ~ f (x_{j_{1}, j_{2}}), where x_{j_{1}, j_{2}} := \frac{1}{⌊ L / 2 ⌋} (j_{1} - ⌊ L / 2 ⌋, j_{2} - ⌊ L / 2 ⌋),

for $j_{1}, j_{2} = 1, \dots, L$ . Let $p := L^{2}$ be the number of pixels in the image $f$ , and let

x_{1}, \dots, x_{p} be an enumeration of {x_{j_{1}, j_{2}} : j_{1}, j_{2} = 1, \dots, L} .

We can consider the $L \times L$ image $f$

as a vector

f \in R^{p}

whose

j

-th entry is defined by

f_{j} = ~ f (x_{j}) .

Similarly, for a fixed bandlimit parameter $λ$ , let

ψ_{1}, \dots, ψ_{m} be an enumeration of {ψ_{n k} : λ_{n k} \leq λ},

where $ψ_{n k}$ is defined in (1) above. With this notation, coefficients of a function in the basis of eigenfunctions up to bandlimit $λ$ is a vector in $R^{m}$ . By Weyl’s law, see for example [17], we have $m \sim λ^{2} / 4$ .

In the following, we consider a linear transformation

B

which maps coefficients in

R^{m}

to images in

R^{p}

, and its adjoint

B^{*}

which maps images in

R^{p}

to coefficients in

R^{m}

1.3. Main result

We consider the linear transform $B : R^{m} \to R^{p}$ which maps coefficients to images by

(B a)_{j} = m \sum i = 1 a_{i} ψ_{i} (x_{j}) h,

(3)

and its adjoint transform $B^{*} : R^{p} \to R^{m}$ which maps images to coefficients by

(B^{*} f)_{i} = p \sum j = 1 f_{j} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯ ¯ ψ_{i} (x_{j}) h,

(4)

where $h = 1 / ⌊ (\sqrt{p} - 1) / 2 ⌋$ is a normalization constant, chosen such that $B^{*} B$ is close to the identity. In this paper, we present a fast and accurate method to apply $B$ and $B^{*}$ . In particular, the presented method can apply the operators $B$ and $B^{*}$ to vectors in $O (p log p)$ operations for any fixed relative error $ε$ . We again emphasize that this is in contrast to previous results, since to the best of our knowledge, existing methods for computing the expansion coefficients in a steerable basis either require $O (p^{3 / 2})$ operations [19, 20] or lack accuracy guarantees [34, 35]. The application of the operators $B$ and $B^{*}$ can be used in an iterative method to determine least squares optimal expansion coefficients for a given image. Alternatively, applying $B^{*}$ to $f$

can be viewed as estimating the continuous inner products that define the coefficients by using quadrature points on a grid (and potentially quadrature weights to provide an endpoint correction).

We report the relative errors of our fast method for applying $B$ and $B^{*}$ compared to matrix multiplication for images with up to $p = 160^{2}$ pixels with varying values of $ε$ , and report timings for images with up to $p = 512^{2}$ pixels. Moreover, we present a numerical example involving rotations, radial convolutions, and deconvolutions.

1.4. Organization

The remainder of the paper is organized as follows. In §2 we describe the analytical apparatus underlying the method. In §3 we describe the computational method. In §4 we justify why the discretization of the analytical apparatus performed by the computational method achieves the specified precision. In §5 we present numerical results. In §6 we discuss implications of the method and potential extensions.

2. Analytical apparatus

2.1. Notation

The eigenfunctions of the Laplacian on the unit disk (that satisfy Dirichlet boundary conditions) defined in (1) can be extended to $R^{2}$ as functions supported on the unit disk by

ψ_{n k} (r, θ) = c_{n k} J_{n} (λ_{n k} r) e^{ı n θ} χ_{[0, 1)} (r),

(5)

for $(n, k) \in Z \times Z_{> 0}$ , where $χ_{[0, 1)}$ denotes an indicator function for $[0, 1)$ . For the sake of completeness, we note that the normalization constants $c_{n k}$ which ensure that $∥ ψ_{n k} ∥_{L^{2}} = 1$ are defined by

\begin{matrix} c_{0 k} = \frac{1}{π^{1 / 2} J_{1} (λ_{0 k})}, and c_{n k} = \frac{2^{1 / 2}}{π^{1 / 2} J_{n + 1} (λ_{n k})}, for n \neq 0. \end{matrix}

(6)

We use the convention that the Fourier transform

ˆ f : R^{2} \to C

of an integrable function

f : R^{2} \to C

is defined by

ˆ f (ξ) = \frac{1}{2 π} \int_{R^{2}} f (x) e^{- ı x \cdot ξ} d x,

(7)

where $x \cdot ξ$ denotes the Euclidean inner product. We define the convolution of two functions $f, g : R^{2} \to C$ by

(f * g) (x) = \int_{R^{2}} f (x - y) g (y) d y .

Furthermore, we will make use of the identity

J_{n} (r) = \frac{1}{2 π} \int_{0}^{2 π} e^{ı r sin θ} e^{- ı n θ} d θ,

(8)

see for example [29, Eq. 9.19]. To verify this identity, recall that $J_{n} (r)$ is the solution to Bessel’s equation

r^{2} \partial_{r r} y + r \partial_{r} y + (r^{2} - n^{2}) y = 0,

(9)

see [9, Eq. 10.2.1]. Observe that

(r^{2} \partial_{r r} + r \partial r + r^{2}) e^{ı r sin θ} = - \partial_{θ θ} e^{ı r sin θ} .

Since taking the second derivative of a function multiplies its $n$ -th Fourier coefficient by $- n^{2}$ , it follows that

(r \partial_{r r} + r \partial_{r} + (r^{2} - n^{2})) \frac{1}{2 π} \int_{0}^{2 π} e^{ı r sin θ} e^{- ı n θ} d θ = 0,

(10)

which establishes that (8) is a solution to (9). It is straightforward to verify that (8) also satisfies the same scaling relations, see [9, Eq. 10.2.2], as $J_{n} (r)$ when $r \to 0$ .

2.2. Fourier transform of eigenfunctions

The analytic foundation for the presented fast method is the following expression for the Fourier transform of the functions $ψ_{n k}$ defined in (5).

Lemma 2.1.

The Fourier transform ${ˆ ψ}_{n k}$ can be expressed by

{ˆ ψ}_{n k} (ξ) = (- ı)^{n} e^{ı n ϕ} \int_{0}^{1} c_{n k} J_{n} (λ_{n k} r) J_{n} (ρ r) r d r,

(11)

where $(ρ, ϕ)$ are polar coordinates for $ξ$ .

We prove this lemma for completeness.

Proof of Lemma 2.1.

By the definition of the Fourier transform (7) we have

{ˆ ψ}_{n k} (ξ) = \frac{1}{2 π} \int_{R^{2}} ψ_{n k} (x) e^{- ı x \cdot ξ} d x .

Changing to polar coordinates $ξ = (ρ cos ϕ, ρ sin ϕ)$ and $x = (r cos θ, r sin θ)$ gives

{ˆ ψ}_{n k} (ξ) = \frac{1}{2 π} \int_{0}^{2 π} \int_{0}^{1} c_{n k} J_{n} (λ_{n k} r) e^{ı n θ} e^{- ı r ρ cos (θ - ϕ)} r d r d θ,

where we used the fact that $x \cdot ξ = r ρ cos (θ - ϕ)$ . Changing variables $θ \mapsto - θ + ϕ - π / 2$ and taking the integral over $θ$ gives

{ˆ ψ}_{n k} (ξ) = (- ı)^{n} e^{ı n ϕ} \int_{0}^{1} c_{n k} J_{n} (λ_{n k} r) J_{n} (ρ r) r d r,

as desired. ∎

2.3. Coefficients from eigenfunction Fourier transform

Next, we observe how the coefficients of a function in the eigenfunction basis can be computed by an application of Lemma 2.1. In the following, we will write the arguments of Fourier transforms of functions in polar coordinates $(ρ, ϕ)$ . We have the following result:

Lemma 2.2.

Suppose that $I \subset Z \times Z_{> 0}$ is a finite index set, and set

f = \sum (n, k) \in I α_{n k} ψ_{n k},

(12)

where $α_{n k} \in C$ are coefficients. Define $β_{n} : [0, \infty) \to R$ by

β_{n} (ρ) := ı^{n} \int_{0}^{2 π} ˆ f (ρ, ϕ) e^{- ı n ϕ} d ϕ .

(13)

It then holds that

α_{n k} = c_{n k} β_{n} (λ_{n k}) .

(14)

The proof is a direct consequence of Lemma 2.1.

Proof of Lemma 2.2.

Observe that (11) implies

ı^{n} \int_{0}^{2 π} {ˆ ψ}_{n^{'} k^{'}} (ρ, ϕ) e^{- ı n ϕ} d ϕ = 2 π δ (n - n^{'}) \int_{0}^{1} c_{n^{'} k^{'}} J_{n} (λ_{n^{'} k^{'}} r) J_{n} (ρ r) r d r,

(15)

where $δ (x) = 1$ if $x = 0$ and $δ (x) = 0$ otherwise. Evaluating (15) at radius $ρ = λ_{n k}$ gives

\begin{matrix} ı^{n} \int_{0}^{2 π} {ˆ ψ}_{n^{'} k^{'}} (λ_{n k}, ϕ) e^{- ı n ϕ} d ϕ = 2 π δ (n - n^{'}) \int_{0}^{1} c_{n^{'} k^{'}} J_{n} (λ_{n^{'} k^{'}} r) J_{n} (λ_{n k} r) r d r = 2 π δ (n - n^{'}) \int_{0}^{1} c_{n k^{'}} J_{n} (λ_{n k^{'}} r) J_{n} (λ_{n k} r) r d r = \frac{1}{c_{n k}} δ (n - n^{'}) δ (k - k^{'}), \end{matrix}

where the final equality follows from the orthogonality of the eigenfunctions $ψ_{n k^{'}}$ , (which is a consequence of the fact that the Laplacian is self-adjoint). By the definition of $β_{n}$ in (13), this implies that

\begin{matrix} β_{n} (λ_{n k}) & = \sum (n^{'}, k^{'}) \in I α_{n^{'} k^{'}} ı^{n} \int_{0}^{2 π} {ˆ ψ}_{n^{'} k^{'}} (ρ, ϕ) e^{- ı n ϕ} d ϕ = \sum (n^{'}, k^{'}) \in I \frac{α_{n^{'} k^{'}}}{c_{n^{'} k^{'}}} δ (n - n^{'}) δ (k - k^{'}) = \frac{α_{n k}}{c_{n k}}, \end{matrix}

which concludes the proof. ∎

Remark 2.1 (Special property of Bessel functions).

We emphasize that the integral expression (8) of the Bessel function is crucial for the fast method of this paper. The possibility of extending the approach to create other fast transforms defined on domains in $R^{2}$ therefore hinges on identifying equally useful integral expressions for the corresponding transforms.

2.4. Convolution with radial functions

Let $g (x) = g (| x |)$ be a radial function. In this section, we observe how the convolution with $g$ can be computed via a diagonal transform of the coefficients. More precisely, we compute the projection of the convolution with $g$ onto the span of any finite basis of the eigenfunctions $ψ_{n k}$ .

Lemma 2.3.

Let $f$ be a function with coefficients $α_{n k}$ as in (12), and $g (x) = g (| x |)$ be a radial function. We have

P_{I} (f * g) = \sum (n, k) \in I α_{n k} ˆ g (λ_{n k}) ψ_{n k},

where $P_{I}$ denotes the orthogonal projection onto the span of ${ψ_{n k}}_{(n, k) \in I}$ .

The proof is a direct application of Lemma 2.2.

Proof of Lemma 2.3.

We use the notation $g (x) = g (| x |)$ and $ˆ g (ξ) = ˆ g (| ξ |)$ . Since the functions $ψ_{n k}$ are an orthonormal basis, in order to compute the orthogonal projection $P_{I}$ , it suffices to determine the coefficients of $f * g$ with respect to $ψ_{n k}$ for $(n, k) \in I$ . Since $ˆ (f * g) (ρ, ϕ) = ˆ f (ρ, ϕ) ˆ g (ρ)$ , and $ˆ g$ is radial, we have

\begin{matrix} ı^{n} \int_{0}^{2 π} ˆ (f * g) (λ_{n k}, ϕ) e^{- ı n ϕ} d ϕ & = ı^{n} \int_{0}^{2 π} ˆ f (λ_{n k}, ϕ) ˆ g (λ_{n k}) e^{- ı n ϕ} d ϕ = \frac{α_{n k}}{c_{n k}} ˆ g (λ_{n k}), \end{matrix}

where the final equality follows from (14). An application of Lemma 2.2 then completes the proof. ∎

3. Computational method

3.1. Notation

Recall that $p$ is the number of pixels in the given image, and $m$ is the number of functions in our expansion. Throughout, we write all computational complexities in terms of $p$ with the understanding that $m \leq p$ . This assumption is not restrictive since it does not make sense to expand an image using more functions than the number of pixels.

Remark 3.1 (Precision of calculations).

The computational complexities also have a (logarithmic) dependence on the reciprocal of the desired precision $ε$ ; for simplicity, we treat $ε$ as a fixed constant, say, $ε = 10^{- 7}$ , and do not include it in complexity statements.

Remark 3.2 (Continuous to discrete).

Throughout this section, for clarity of exposition, we use the same notation used in the continuous setting described in §2 to represent the corresponding discrete objects. We justify why the discrete sums of the algorithm will approximate the continuous integrals to the desired precision in §4.

3.2. Overview

In the following, we describe how to apply the operators $B$ and $B^{*}$ defined above in §1.3 in $O (p log p)$ operations. For the purpose of exposition, we start by describing a simplified method before presenting the full method. The section is organized as follows:

In §3.3 we describe a simplified method to apply $B$ and $B^{*}$ in $O (p^{3 / 2} log p)$ operations. The simplified method is a direct application of the lemmas from the previous section.
In §3.4 we provide an informal description of how to modify the simplified method to create a fast method to apply $B$ and $B^{*}$ in $O (p log p)$
operations. The main additional ingredient is fast interpolation from Chebyshev nodes.
In §3.5 we give a detailed description of the fast method to apply $B$ and $B^{*}$ in $O (p log p)$ operations.

3.3. Simplified method

In this section, for the purpose of exposition, we present a simplified method that applies $B$ and $B^{*}$ in $O (p^{3 / 2} log p)$ operations. Recall that $f_{1}, \dots, f_{p} \in R$ are the pixel values of the image, $x_{1}, \dots, x_{p} \in [- 1, 1]^{2}$ are the pixel coordinates, and $ψ_{1}, \dots, ψ_{m}$ are the eigenfunctions below the given bandlimit, see §1.2. Let $λ_{1}, \dots, λ_{m}$ denote the Bessel function roots (square root of the eigenvalues) corresponding to $ψ_{1}, \dots, ψ_{m}$ . Further, suppose that $(n_{i}, k_{i})$ are the indices such that $ψ_{i} = ψ_{n_{i} k_{i}}$ and $λ_{i} = λ_{n_{i} k_{i}}$ .

Fix a positive integer $s = ⌊ max {7.08 \sqrt{p}, {log}_{2} ε^{- 1} + {log}_{2} p} ⌋$ (see §4.3 for a justification of this parameter choice). We first describe how to apply $B^{*} : R^{p} \to R^{m}$ which is defined by

(B^{*} f)_{i} = p \sum j = 1 f_{j} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯ ¯ ψ_{i} (x_{j}) h,

where $h$ is a normalization parameter. The simplified method has three steps, which make use of the non-uniform fast Fourier transform (NUFFT) [11, 14, 22].

[label=Step 0., parsep=2ex,wide=0pt]
Using the type-2 2-D NUFFT, in $O (p^{3 / 2} log p)$ operations, compute:

${^f}_{i ℓ} := p \sum j = 1 f_{j} e^{- ı x_{j} \cdot ξ_{i ℓ}} where ξ_{i ℓ} := λ_{i} (cos 2 π ℓ / s, sin 2 π ℓ / s),$

for $(i, ℓ) \in {1, \dots, m} \times {0, \dots, s - 1}$ to precision $ε$ .
Using the FFT, in $O (p^{3 / 2} log p)$ operations, compute:

$β_{n} (λ_{i}) := s \sum ℓ = 1 {^f}_{i ℓ} e^{- ı n 2 π ℓ / s},$

for $(i, n) \in {1, \dots, m} \times {0, \dots, s - 1}$ .
By Lemma 2.2 it follows that

$(B^{*} f)_{i} = β_{n_{i}} (λ_{i}) c_{i} h,$

where $n_{i}$ is the integer such that $ψ_{i} = ψ_{n_{i} k_{i}}$ and $c_{i} = c_{n_{i} k_{i}}$ .

In the following remarks, we discuss how reversing these steps provides an immediate way to apply $B$ and discuss limitations of the simplified method.

Remark 3.3 (Applying $B$ ).

Note that each step of the algorithm consists of applying a linear transform whose adjoint can be applied in a similar number of operations.

Indeed, the adjoint of Step 1 (which uses type-2 NUFFT) is a type-1 2-D NUFFT [1], and the adjoint of Step 2 (which uses a standard FFT) is an inverse FFT. Step 3 (interpolation) can be computed in a variety of ways (including by using the NUFFT or FFT and sparse matrices, see Remark 3.5) which each have similarly fast adjoints [11, 14, 22]. Therefore, in order to apply $B : R^{m} \to R^{p}$ , which is defined by

(B a)_{j} = m \sum i = 1 a_{i} ψ_{i} (x_{j}) h,

we can just apply the adjoint of each step in the reverse order.

Remark 3.4 (Limitations of simplified method).

The problem with the simplified method is that the first step is prohibitively expensive since we want a method that works in $O (p log p)$ operations. The key to overcome this limitation is to observe that $β_{n_{i}} (λ_{i})$ are samples from analytic functions, and to add the additional ingredient of fast interpolation from Chebyshev nodes.

3.4. Summary of fast method

In this section, we describe how the method of the previous section can be improved from $O (p^{3 / 2} log p)$ to $O (p log p)$ by using fast interpolation from Chebyshev nodes. Recall that Lemma 2.2 defines

β_{n} (ρ) := ı^{n} \int_{0}^{2 π} ˆ f (ρ, ϕ) e^{- ı n ϕ} d ϕ,

and states that the coefficients can be computed by evaluating the appropriate $β_{n}$ at $λ_{i}$ . The problem with the simplified method is the first step: computing $ˆ f (ρ, θ)$ for $ρ = λ_{1}, \dots, λ_{m}$ and a number of points $θ$ so that we are sampling at the Nyquist rate is a set of $O (p^{3 / 2})$ points, which is already prohibitive. Fortunately, there is a simple solution to this issue: for each $n$ , the functions $β_{n} (ρ)$ are analytic functions of $ρ$ which can be interpolated to determine the coefficients, if they are tabulated at appropriate points; see §4.2 for a discussion of why $O (\sqrt{p})$ points are sufficient. Here is an informal summary of the fast method:

Compute the Fourier transform of $f$ at

$ξ_{i ℓ} := t_{i} (cos 2 π ϕ_{ℓ}, sin 2 π ϕ_{ℓ}),$

for $O (\sqrt{p})$ Chebyshev nodes $t_{i}$ and $O (\sqrt{p})$ angles $ϕ_{ℓ}$ .
Compute $β_{n} (t_{i})$ for the $O (\sqrt{p})$ Chebyshev nodes $t_{i}$ and $O (\sqrt{p})$ frequencies $n$ .
For each of the $O (\sqrt{p})$ frequencies, use fast interpolation from the $O (\sqrt{p})$ Chebyshev nodes $t_{i}$ to the $O (\sqrt{p})$ Bessel function roots associated with each frequency $n$ . We illustrate the interpolation step in Fig. 2.

Figure 2. We visualize the interpolation step for a $64 \times 64$ input image. For $n = 40, 100$ we plot $β_{n} (ρ)$ (black line), interpolation source nodes (black dots), and target points (orange crosses).

3.5. Detailed description of fast method

Recall that $ψ_{1}, \dots, ψ_{m}$ are the eigenfunctions below our bandlimit, and $λ_{1}, \dots, λ_{m}$ are the associated Bessel function roots (which are the square root of the associated eigenvalues). Assume that $λ_{1}, \dots, λ_{m}$ are listed in ascending order so that they are contained in the interval $[λ_{1}, λ_{m}]$ . Fix $q = ⌈ max {2.4 \sqrt{p}, {log}_{2} ε^{- 1}} ⌉$ , where $ε$ is the desired relative error, see §4.2 for a justification of this parameter choice. Recall that Chebyshev nodes of the first kind for the interval $[λ_{1}, λ_{m}]$ are defined by

t_{i} = (λ_{m} - λ_{1}) \frac{1 - cos (π (i - 1) / q)}{2} + λ_{1},

for $i = 1, \dots, q$ . The fast method follows the same strategy as the simplified method. The key difference is that the fast method replaces the final step with fast interpolation, and therefore evaluates the Fourier coefficients on a smaller number of nodes (the Chebyshev nodes).

[label=Step 0., parsep=2ex, wide=0pt]
Using the NUFFT, in $O (p log p)$ operations, compute:

${^f}_{i ℓ} := p \sum j = 1 f_{j} e^{- ı x_{j} \cdot ξ_{i ℓ}} where ξ_{i ℓ} := t_{i} (cos 2 π ℓ / s, sin 2 π ℓ / s),$

for $(i, ℓ) \in {1, \dots, q} \times {0, \dots, s - 1}$ .
Using the FFT, in $O (p log p)$ operations, compute:

$β_{n} (t_{i}) := s \sum ℓ = 1 {^f}_{i ℓ} e^{- ı n 2 π ℓ / s},$

for $(i, n) \in {1, \dots, m} \times {0, \dots, s - 1}$ .
Using fast interpolation from Chebyshev nodes, in $O (p log p)$ operations, compute:

$(B^{*} f)_{i} = β_{n_{i}} (λ_{i}),$

for $i \in {1, \dots, m}$ where $n_{i}$ is the integer such that $ψ_{i} = ψ_{n_{i} k_{i}}$ .

As detailed in Remark 3.3, the fact that each step is a linear operator whose adjoint can be applied in a similar number of operations provides a similar fast method to apply $B$ .

Remark 3.5 (Methods for fast interpolation from Chebyshev nodes).

Theoretically, the most straightforward way to perform fast interpolation from Chebyshev nodes is to consider the function values at the Chebyshev nodes as function values associated with equally spaced points on a half circle (by a change of variables). Then the fast interpolation (called spectral interpolation) can be computed by using the NUFFT in $O (p log p)$ operations or from arbitrary nodes using the Fast Multipole Method, see [10]. Practically speaking, since the NUFFT run-time constant is high, choosing a fixed number of source points centered around the target points (say $20$ source points) and then applying a precomputed sparse (barycentric interpolation [3]) matrix may be more practical; sparse interpolation can be used in combination with standard spectral interpolation (using the discrete cosine transform) to double or quadruple the number of Chebyshev nodes before the sparse interpolation step.

4. Continuous to discrete

We conclude the description of the fast algorithm by discussing additional analytic details that justify transitioning from the continuous identities in §2 to the discrete implementations of §3.

4.1. Maximum bandlimit

In this section, we derive a bound for the bandlimit parameter, which controls the number of basis functions used in the expansion. The number of basis functions should not exceed the number of pixels of the image corresponding to points within the unit disk. In combination with Weyl’s law, this leads to a bound on the maximum bandlimit.

In more detail, from Weyl’s law, the number of Dirichlet eigenvalues of the Laplacian in the disk is bounded by $λ$ is

(16)

see [8]. Similarly, the inscribed circle in an $L \times L$ image has radius $⌊ \frac{L - 1}{2} ⌋$ . The number of pixels contained inside this circle is therefore

# {(i, j) \in Z \times Z : i^{2} + j^{2} \leq {⌊ \frac{L - 1}{2} ⌋}^{2}} = π {⌊ \frac{L - 1}{2} ⌋}^{2} + O (L^{2 / 3}),

(17)

see for example [13, 16, 30]. Equating (16) with (17) results in

λ = 2 \sqrt{π} ⌊ \frac{L - 1}{2} ⌋ + 1 + O (L^{- 1 / 3}) .

(18)

Practically speaking, it can be advantageous to expand the image in fewer basis functions than described by (18), but, in any case, $λ = O (L)$ is a natural choice and different motivations exist in the literature [34].

Remark 4.1 ( FFT Bandlimit heuristic).

One heuristic for setting the bandlimit is based on the fast Fourier transform (FFT). For a centered FFT on a signal of length

L

, the maximum frequency is

π^{2} (L / 2)^{2}

which corresponds to a bandlimit of

λ = π L / 2

. Note that

π L / 2 \approx 1.57 L < 1.77 L \approx \sqrt{π} ⌊ (L - 1) / 2 ⌋ / 2,

(19)

so this FFT bandlimit heuristic does indeed produce a reasonable bandlimit below the bound (18) derived from Weyl’s law. We use this bandlimit for our numerical experiments. The computational complexity and accuracy guarantees of the method presented in this paper hold for any bandlimit $λ = O (L)$ . However, the fact that the fast method performs interpolation in Fourier space inside a disk bounded by the maximum bandlimit provides additional motivation for this FFT-based heuristic since it will ensure that the disk will be contained within the square in frequency space used by the two-dimensional FFT.

4.2. Number of radial nodes

The following lemma shows that $O (\sqrt{p})$ Chebyshev nodes are sufficient for accurate interpolation.

Lemma 4.1.

When running the algorithm on an image with $p$ pixels at the maximum bandlimit, choosing the number of Chebyshev nodes $q$ to satisfy

q = ⌊ max {2.4 \sqrt{p}, {log}_{2} ε^{- 1}} ⌋,

(20)

achieves relative error less than $ε$ in the interpolation step of §3.5.

It will be clear from the proof that the constant $2.4$ in the statement of the lemma is an overestimate, see Remark 4.2 for a discussion of how this constant can be improved.

Proof of Lemma 4.1.

When interpolating a smooth differentiable function $h$ defined on the interval $[a, b]$ using an interpolating polynomial $p_{q}$ at $q$ Chebyshev nodes, the residual term $R (x) = h (x) - p_{q} (x)$ can be written as

| R (x) | \leq \frac{2 C_{q}}{q!} {(\frac{b - a}{4})}^{q},

where $C_{q} := {max}_{x \in [a, b]} | h^{(q)} (x) |$ ; indeed, this can be deduced from [9, § 3.3]. If we apply this result with $[a, b] = [λ_{1}, λ_{m}]$ , the residual satisfies

| R (x) | \leq \frac{2 C_{q}}{q!} {(\frac{λ_{m} - λ_{1}}{4})}^{q} \leq \frac{2 C_{q}}{q!} {(\frac{\sqrt{π p}}{4})}^{q}

where the final inequality follows from the bound $λ_{m} \leq \sqrt{π p}$ ; see § 4.1. In order to apply this bound to $h (ρ) = β_{n} (ρ)$ , we estimate

C_{q} := max ρ \in [λ_{1}, λ_{m}] | β_{n}^{(q)} (ρ) |,

where $β_{n}^{(q)}$ denotes the $q$ th derivative of $β_{n}$ defined in (13). Since $f$ is compactly supported, its Fourier transform $^f$ is analytic and the $C_{q}$ are therefore all finite. More precisely, the term $C_{q}$ can be bounded uniformly by using

| β_{n}^{(q)} (ρ) | = ∣ ∣ ∣ ∣ \frac{d^{q}}{d ρ^{q}} (\frac{}{1} 2 π \int_{0}^{2 π} ˆ f (ρ, ϕ) e^{- i n ϕ} d ϕ) ∣ ∣ ∣ ∣ \leq (\frac{1}{2 π} \int_{0}^{2 π} ∣ ∣ ∣ \frac{d^{q}}{% d ρ^{q}} ˆ f (ρ, ϕ) ∣ ∣ ∣ d ϕ) .

We write out

\begin{matrix} ∣ ∣ ∣ \frac{d^{q}}{d ρ^{q}} ˆ f (ρ, ϕ) ∣ ∣ ∣ & = ∣ ∣ ∣ \frac{1}{2 π} \frac{d^{q}}{d ρ^{q}} \int_{0}^{2 π} \int_{0}^{\infty} f (r, θ) e^{- ı r ρ cos (θ - ϕ)} r d r d θ ∣ ∣ ∣ = ∣ ∣ ∣ \frac{1}{2 π} \int_{0}^{2 π} \int_{0}^{\infty} {(- ı r cos (θ - ϕ))}^{q} f (r, θ) e^{- ı r ρ cos (θ - ϕ)} r d r d θ ∣ ∣ ∣ = ∣ ∣ ∣ \frac{1}{2 π} \int_{0}^{2 π} \int_{0}^{1} {(- ı r cos (θ - ϕ))}^{q} f (r, θ) e^{- ı r ρ cos (θ - ϕ)} r d r d θ ∣ ∣ ∣ \leq \frac{1}{2 π} \int_{0}^{2 π} \int_{0}^{1} | f (r, θ) | r d r d θ = \frac{1}{2 π} ∥ f ∥_{L^{1}} \leq \frac{1}{2 \sqrt{} π} ∥ f ∥_{L^{2}}, \end{matrix}

where the third equality uses the fact that $f$ is supported on the disk of radius $1$ and the last inequality uses Cauchy-Schwarz applied to $f$ multiplied by the indicator function of the unit disk. It follows that

Therefore, in order to achieve relative error $| R (x) | / ∥ f ∥_{L^{2}} = ε$ , it suffices to set $q$ such that

ε \leq \frac{1}{\sqrt{π} q!} {(\frac{\sqrt{π p}}{4})}^{q} .

(21)

It follows from Stirling’s approximation [9, 5.11.3] that

\frac{| R (x) |}{∥ f ∥_{L^{2}}} \leq {(\frac{\sqrt{π p} e}{4 q})}^{q}

(22)

Setting $\sqrt{π p} e / 4 q = 1 / 2$ and solving for $q$ gives

q = \frac{\sqrt{π} e \sqrt{p}}{2} \approx 2.4 \sqrt{p} ⟹ q \geq max {2.4 \sqrt{p}, {log}_{2} ε^{- 1}),

is sufficient to achieve relative error less than $ε$ . ∎

4.3. Number of angular nodes

The following lemma shows that $s = O (\sqrt{p})$ angular nodes is sufficient to achieve relative error $ε$ .

Lemma 4.2.

When running the algorithm on an image with $p$ pixels at the maximum bandwidth, choosing the number of equispaced angular nodes $s$ to satisfy

s = ⌊ max {7.08 \sqrt{p}, {log}_{2} ε^{- 1} + {log}_{2} p} ⌋ .

(23)

achieves relative error less than $ε$ in the FFT step of §3.5.

As above, we emphasize that the constant $7.08$ in the statement of this result is an overestimate. See 4.2 for discussion about how this constant can be improved.

Proof of Lemma 4.2.

We will require a classical result, see for example [18]: suppose that $f : [0, 2 π] \to R$ is a smooth periodic function on the torus $[0, 2 π]$ where $0$ and $2 π$ are identified. Then

∣ ∣ ∣ ∣ \int_{0}^{2 π} f (θ) d θ - \frac{1}{n} n - 1 \sum j = 0 f (θ_{j}) ∣ ∣ ∣ ∣ < 2 \frac{∥ f^{(n)} ∥_{2}}{n^{n}},

for all $n \geq 2$ , where $θ_{j} = 2 π j / n$ .

Indeed, we start by writing $f (θ) = \sum_{k = - \infty}^{\infty} a_{k} e^{ı k θ}$ and use the fact that $\frac{1}{n} \sum_{j = 0}^{n - 1} e^{ı k θ_{j}}$ is equal to $1$ if $k$ is a multiple of $n$ and is $0$ otherwise to obtain:

∣ ∣ ∣ ∣ \int_{0}^{2 π} f (θ) d θ - \frac{1}{n} n - 1 \sum j = 0 f (θ_{j}) ∣ ∣ ∣ ∣ = ∣ ∣ ∣ ∣ a_{0} - \frac{1}{n} n - 1 \sum j = 0 \infty \sum k = - \infty a_{k} e^{ı k θ_{j}} ∣ ∣ ∣ ∣ = ∣ ∣ ∣ ∣ \sum t \neq 0 a_{n t} ∣ ∣ ∣ ∣ .

Multiplying and dividing terms by $(ı n t)^{n}$ , and using Cauchy-Schwarz:

\begin{matrix} ∣ ∣ ∣ ∣ \sum t \neq 0 a_{n t} \frac{(ı n t)^{n}}{(ı n t)^{n}} ∣ ∣ ∣ ∣ & \leq {⎛ ⎝ \sum t \neq 0 | (ı n t)^{n} a_{n t} |^{2} ⎞ ⎠}_{n t}^{1 / 2} {(\sum t > 0 \frac{2}{(n t)^{2 n}})}^{1 / 2} \leq {(\infty \sum ℓ = - \infty | (ı ℓ)^{n} a_{ℓ} |^{2})}_{ℓ}^{1 / 2} {(\sum t > 0 \frac{2}{(n t)^{2 n}})}^{1 / 2} < 2 \frac{∥ f^{(n)} ∥_{L^{2}}}{n^{n}}, \end{matrix}

where the last inequality follows from Parseval’s identity and the identity $\sum_{t > 0} \frac{1}{t^{4}} = \frac{π^{4}}{90}$ , together with $\sqrt{π^{4} / 45} < 2$ .

In order to determine the required number of angular nodes we write the image $f$ and its Fourier transform $ˆ f$ as

f (x) = p \sum j = 1 f_{j} δ_{x_{j}} (x), ˆ f (ξ) = p \sum j = 1 f_{j} \frac{1}{2 π} e^{i ξ \cdot x_{j}} .

By the triangle inequality, it suffices to estimate integrals of functions of the form

g (ϕ) = e^{i r ρ cos (θ - ϕ) - i n ϕ}

to precision $\frac{ε}{p}$ for all $r \in [0, 1)$ , $θ \in [0, 2 π)$ , $ρ \in [λ_{1}, λ_{m}]$ , and $| n | \leq n_{m}$ . We have the estimate

| g^{(s)} | \leq (λ_{m} + n_{m})^{s} .

Thus, in order to integrate these functions, we need to choose $s$ such that

{(\frac{λ_{m} + n_{m}}{s})}^{s} \leq \frac{ε}{p} .

(24)

It follows that choosing

s = max {2 (λ_{m} + n_{m}), {log}_{2} ε^{- 1} + {log}_{2} p},

achieves relative error $ε$ . To complete the proof we note that $λ_{m} \leq λ$ , where $λ$ is the maximum bandlimit from § 4.1. Also by [9, 10.21.40] we have

λ_{n 1} = n + 1.8575 n^{1 / 3} + O (n^{- 1 / 3}),

which implies that the maximum angular frequency $n_{m} \leq λ$ . We conclude that

s = max {4 λ, log ε^{- 1} + {log}_{2} p},

is sufficient to achieve error $ε$ ; using the bound $λ \leq 1.77 \sqrt{p}$ from (19) completes the proof. ∎

Remark 4.2 (Improving estimates for number of radial and angular nodes).

While Lemmas 4.1 and 4.2 show that the number of radial nodes $q$ and angular nodes $s$ are $O (\sqrt{p})$ , the constants in the lemmas are not optimal. For practical purposes, choosing the minimal number of nodes possible to achieve the desired error is advantageous to improve the run time constant of the algorithm, and it is clear from the proofs how the estimates can be refined. For Lemma 4.1 we set $Q = ⌈ 2.4 \sqrt{p} ⌉$ , and motivated by (21) compute

ε_{radial} (q) = \frac{1}{\sqrt{π} q!} {(\frac{\sqrt{π p}}{} 4)}^{q},

for $q = 1, \dots, Q$ and choose the smallest value $q^{*}$ of $q$ such that $ε_{radial} (q^{*}) \leq ε$ . Similarly, for Lemma 4.2, we set $S = ⌈ 7.08 \sqrt{p} ⌉$ , and motivated by (24) compute

ε_{angular} (s) = p {(\frac{λ_{m} + n_{m}}{s})}^{s},

for $s = 1, \dots, S$ and choose the smallest value $s^{*}$ of $s$ such that $ε_{angular} (s) \leq ε$ . Then, it follows that $2.4 \sqrt{p}$ and $7.08 \sqrt{p}$ can be replaced by $q^{*}$ and $s^{*}$ , in the statements of Lemmas 4.1 and 4.2, respectively. This procedure improves the estimate of the required number of angular and radial nodes by a constant factor.

5. Numerical results

5.1. Accuracy results

In this section, we report numerical results for the accuracy of the fast algorithm compared to matrix multiplication. Recall that $B : R^{m} \to R^{p}$ maps coefficients to images by

(B a)_{j} = m \sum i = 1 a_{i} ψ_{i} (x_{j}) h,

and its adjoint transform $B^{*} : R^{p} \to R^{m}$ maps images to coefficients by

(B^{*} f)_{i} = p \sum j = 1 f_{j} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯ ¯ ψ_{i} (x_{j}) h,

see § 1.3. By defining the $m \times p$ matrix $B$ by

B (i, j) = ψ_{i} (x_{j}) h,

we can apply $B$ and $B^{*}$ by dense matrix multiplication to test the accuracy of our fast method. Since the size of the matrix scales like $L^{4}$ for $L \times L$ images, constructing these matrices quickly becomes prohibitive so the comparison is only given up to $L = 160$ , see Table 1, where

{err}_{a} = ∥ a_{fast} - a_{dense} ∥