Multitaper estimation on arbitrary domains

1 Introduction

Estimating the frequency content of a stochastic process is crucial to many data processing tasks. For example, in several inverse problems, such as denoising, access to a good noise model is necessary for accurate reconstruction. With other tasks, such as system identification, this frequency structure is itself the object of interest.

The frequency content of a stationary process $X$ over $Z^{d}$ is characterized by its spectral density $S$ , defined as the Fourier series of its autocovariance function [30]. Spectral estimation is the task of recovering $S$ given one or more realizations of $X$ over a finite subset $Ω$ of $Z^{d}$ . The challenge in the estimation problem is two-fold: (i) stochastic fluctuations in $X$

introduce variance into estimates of

S

, and (ii) spatial constraints restrict us to using only samples from the domain

Ω

. When many realizations are available, the first point may be addressed through ensemble averaging. For several applications, however, such as geosciences, cosmology, and electron microscope imaging, only a single realization is available, requiring accurate single-shot estimators.

One approach to reducing error in single-shot spectral estimation is to divide the domain $Ω$ into disjoint subsets, compute spectral density estimates for each, and average the result. In one dimension, with $Ω = {0, \dots, N - 1}$ for some $N > 0$ , this is typically done by partitioning $Ω$ into $K$ blocks of $N / K$ samples each, computing a spectral estimate, such as a periodogram, for each subset, and averaging. The $K$ parameter controls the trade-off between the two challenges raised previously: while a higher value of $K$ mitigates the stochastic fluctuations (reducing variance), each periodogram is computed from only $N / K$ samples, reducing their resolution (increasing bias) [5]. One problem with this approach, however, is that it introduces artifacts due to the boundaries imposed by the partitioning. The multitaper estimator, introduced by Thomson [42], refines this method by instead computing periodograms over all of $Ω$ , but with the samples first multiplied by a set of $K$ sequences called tapers. Optimal trade-off between bias and variance is obtained for specific sequences known as discrete prolate spheroidal sequences, or Slepian tapers [40].

Originally defined for $d = 1$ , multitaper estimators naturally generalize to higher dimensions [38]

. The Slepian tapers, however, are defined as solutions to an ill-posed eigenvalue problem. Special geometries, including rectangular grids

[20, Chapter 2] and domains involving other symmetries [15, 20, 37, 36], may be analyzed by means of commuting differential operators, yielding well-posed problems that allow for accurate calculation of the tapers. There are other stable algorithms for calculating Slepian tapers, but these are defined only for specific domains [24, 25, 21, 23] or for certain modified tapers [32].

For general domains $Ω$ , however, all available methods inherit the instability of the underlying eigenvalue problem [37]. As a result, high-precision multitaper estimators are not applicable to domains such as the complement of the disk (see Figure 0(c)). This geometry arises in cryo-electron microscopy (cryo-EM), where noisy projection images contain signal on a central disk and clean samples of the noise process are found only outside of that disk. In geosciences, a similar problem occurs when a physical quantity needs to be sampled on a subregion $Ω$ of the earth [36, 37, 31].

One solution to this problem is to partition this irregular domain into smaller rectangular regions and apply the tensor Slepian multitaper estimator to those. This strategy, however, lacks the simplicity Thomson’s multitaper estimator, and may typically only be implemented suboptimally.

In this work, we show that the multitaper estimator only depends on the linear span of the tapers used. Consequently, we may calculate it with any set of tapers that span the same subspace as the Slepian tapers. The ill-posed problem of computing the leading $K$ eigenvectors is therefore replaced with computing a basis for that eigenspace. This problem is well-posed, because of the large spectral gap between the first $K$ eigenvalues and the rest of the spectrum, which we validate.

Furthermore, we provide a mean squared error bound for the multitaper estimator on arbitrary domains. This generalizes previous results of Abreu and Romero [3] on performance bounds for the one-dimensional case.

The proposed method is evaluated on synthetic data, where it is shown to perform comparably to the Slepian multitaper method for rectangular domains. However, we also show that the method performs equally well on the complement of a disk. We also evaluate the method on cryo-EM images, both synthesized and from experimental datasets, where we also achieve good performance compared to previous approaches.

Section 2 introduces the spectral estimation problem for general domains. The multitaper spectral estimator is introduced in Section 3. Section 4 provides a bound for the mean squared error of the multitaper estimator for arbitrary domains. We show that the estimator only depends on the linear span of the tapers in Section 5 and use this to provide an implementation using randomized tapers. Finally, Section 6 illustrates the performance of the proposed estimator using numerical experiments.

Notation

For a vector

x \in R^{d}

, we let

| x | := (\sum_{k = 1}^{d} {| x_{k} |}_{k}^{2})^{1 / 2}

be its Euclidean norm. The

p

-norms are denoted

{| x |}_{p}

. For two non-negative functions

f, g : X \to [0, + \infty)

, we write

f ≲ g

is there exists a constant

C \geq 0

such that

f \leq C g

. We also write

f ≍ g

, if

f ≲ g

and

g ≲ f

. For a function

f \in ℓ^{1} (Z^{d})

, we denote its

ℓ^{1}

-norm

\sum_{q} | f [q] |

| f |

We say that a function $f : R^{d} \to C$ is 1-periodic if $f (x + k) = f (x)$ for all $k \in Z^{d}$ . The convolution of two periodic functions $f, g$ is defined as

f * g (x) = \int_{[- 1 / 2, 1 / 2]^{2}} f (y) g (x - y) d y .

(1)

The $L^{\infty}$ -norm of a measurable, 1-periodic function $f : R^{d} \to C$ is $∥f∥∞:=\operatornamewithlimitsesssupx∈[−1/2,1/2]d|f(x)|$ , and its $L^{1}$ -norm is $∥ f ∥_{1} = \int_{[- 1 / 2, 1 / 2]^{d}} | f (x) | d x$ . For clarity, we will sometimes write $∥ f ∥_{L^{1} ([- 1 / 2, 1 / 2]^{d})}$ .

The $C^{2}$ -norm of a twice continuously differentiable function $f$ is

∥ f ∥_{C^{2}} := max {∥ f ∥_{\infty}, ∥ \partial_{x_{j}} f ∥_{\infty}, ∥ \partial_{x_{j}} \partial_{x_{j^{'}}} f ∥_{\infty} : j, j^{'} = 1, \dots, d} .

2 Spectral estimation on irregular domains

Let us consider a real-valued Gaussian, zero-mean, stationary process $X$ defined on an infinite grid $Z^{d}$ . Our goal is to estimate the covariance matrix

Cov [q, q^{'}] = E {X [q] X [q^{'}]} q, q^{'} \in Z^{d} .

Since $X$ is stationary, $Cov [q, q^{'}]$ only depends on the difference $q - q^{'}$ . We therefore rewrite the matrix in terms of the autocovariance function $r$ , giving

Cov [q, q^{'}] = r [q - q^{'}] q, q^{'} \in Z^{d} .

The covariance information is equivalently encoded in the spectral density

S (ξ) = \sum q \in Z^{d} r [q] e^{2 π i q ξ}, ξ \in [- 1 / 2, 1 / 2]^{d} .

Given one or several realizations of $X$ on a subset $Ω \subseteq Z^{d}$ of cardinality $n_{Ω}$ , we wish to estimate $S$ . This is known as the spectral estimation problem [30].

(a)

(b)

Figure 1: (a, b) Two sample single-particle cryo-EM images from the EMPIAR-10028 dataset [43]. The images are projections of the Plasmodium falciparum 80S ribosome. (c) The domain used to estimate the noise spectral density.

An instance of this problem arises in imaging. Suppose that an image is defined on a grid

Q_{N} = {0, \dots, N - 1}^{2}

of pixels, and $X [q]$ represents additive noise at the pixel $q \in Q_{N}$ . An accurate estimate of the noise distribution is then needed for denoising and other tasks. For example, in single-particle cryo-EM, molecules are imaged by freezing them in a thin layer of ice and recording their tomographic projections using an electron microscope [13]. To reduce specimen damage, electron dose is kept low, resulting in exceptionally noisy images, as shown in Figures 0(a) and 0(b). The projection of the molecule is expected to lie in the center of the image, so pure noise samples are only found outside of a central disk. A reasonable domain for noise estimation is therefore the complement of that disk

Ω := {q \in Q_{N} : | q - (N / 2, N / 2) | > R},

(2)

as illustrated in Figure 0(c).

Another instance is found in geosciences, where $X$ is a physical quantity on an approximately flat portion of the earth, and measurements are only available on a subregion $Ω$ , corresponding, for example, to a continent [37]. A more sophisticated model replaces the Cartesian grid for a grid on the sphere [36, 31].

For both of these applications, we only have access to a single realization of the process whose spectral density we wish to estimate. In cryo-EM, this is due to changing experimental conditions between projection images, such as non-stationary optical parameters or variation in ice thickness, while in geosciences, only one realization (i.e., one planet) exists. We therefore require single-shot estimators which provide accurate estimates from a single realization of $X$ . The aim of the present work is to study the performance one such method, the multitaper estimator, and provide a high-precision implementation.

3 Multitaper estimators

The inherent variability of $X$ induces error into any estimate of its spectral density $S$

. With access to only a single realization, ensemble averaging cannot be used to reduce the error. Instead, we must impose some constraint on the estimate. This often involves post-processing of the spectral density estimate by smoothing or fitting parametric models

[9, 30].

The multitaper spectral estimator provides another solution to this challenge [42]. Initially introduced for one-dimensional signals, we present it here for signals of arbitrary dimension. Let $m \in ℓ^{2} (Z^{d})$ be a function of unit $ℓ^{2}$ -norm supported on $Ω$ , that is, satisfying

m [q] = 0, if q \notin Ω .

(3)

We now define the tapered periodogram of $X$ with taper $m$ as

{ˆ S}_{m} (ξ) := {∣ ∣ ∣ ∣ \sum q \in Ω m [q] X [q] e^{2 π i q ξ} ∣ ∣ ∣ ∣}^{2}, ξ \in [- 1 / 2, 1 / 2]^{d} .

(4)

Given an orthonormal family $m_{1}, \dots, m_{K} \in ℓ^{2} (Z^{d})$ of tapers supported on $Ω$ , we define the corresponding multitaper estimator as the average

{ˆ S}^{m t} (ξ) := \frac{1}{K} K \sum k = 1 ˆ S_{m_{k}} (ξ) ξ \in [- 1 / 2, 1 / 2]^{d} .

(5)

The accuracy of ${ˆ S}^{m t}$ in estimating $S$ depends on the choice of tapers $m_{1}, \dots, m_{K}$ .

In the traditional multitaper estimator, the tapers are defined using the following spectral concentration problem:

		$maximize \int_{[- W / 2, W / 2]^{d}} {∣ ∣ ∣ ∣ \sum q \in Ω m [q] e^{2 π i q ξ} ∣ ∣ ∣ ∣}^{2} d ξ,$		(6)
		$subject to: \sum q \in Ω {\| m [q] \|}^{2} = 1, and supp (m) \subseteq Ω,$		(6)

where

W = (K / n_{Ω})^{1 / d},

(7)

and $n_{Ω}$ is the number of elements in $Ω$ . The parameter $W$ is known as the bandwidth of the estimator. It is chosen so that

K = ⌈ n_{Ω} \times W^{d} ⌉,

(8)

and controls the amount of smoothness imposed on the estimate ${ˆ S}^{m t}$ .

The Slepian tapers $m_{1}, \dots, m_{K}$ are defined as the set of mutually orthogonal solutions to (6). Specifically, the first taper $m_{1}$ is the solution to (6), while $m_{2}$ is the solution to (6) with the constraint that $m_{2} ⊥ m_{1}$ , and so on. In general, $m_{j}$ is the solution to (6) subject to $m_{j} ⊥ m_{1}, \dots, m_{j} ⊥ m_{j - 1}$ .

Alternatively, (6) may be formulated as the maximization of the quadratic form corresponding to the $d$ -Toeplitz matrix

T^{Ω, W} [q, q^{'}] = W^{d} 1_{Ω} (q) s i n c (W (q - q^{'})) 1_{Ω} (q^{'}), q, q^{'} \in Z^{d},

(9)

where $sinc$ is the $d$ -dimensional normalized sinc function

sinc (u) := d \prod k = 1 \frac{sin (π u_{k})}{π u_{k}}, u = (u_{1}, \dots, u_{d}) \in Z^{d} .

(10)

The Slepian tapers $m_{1}, \dots, m_{K}$ are thus the top $K$ eigenvectors of $T^{Ω, W}$ , satisfying

T^{Ω, W} m_{k} = λ_{k} m_{k}, k \in {1, \dots, K},

(11)

where $λ_{1} \geq λ_{2} \geq \dots \geq λ_{K}$ .

Figure 2: Eigenvalues of $T^{Ω, W}$ for $d = 1$ , $K = 7$ , and $Ω = {1, \dots, 32}$ .

Unfortunately, serious numerical difficulties arise when attempting solve (11) numerically. Indeed, the first $\sim K$ eigenvalues form a plateau profile, all clustering around $1$ (see Figure 2), resulting in small spectral gaps between successive eigenvalues. For increasing $K$ , small spectral gaps result in an ill-posed eigenvector problem [41, 14], making direct calculation of $m_{1}, \dots, m_{K}$ in finite precision very challenging.

Under certain circumstances, $T^{Ω, W}$ may be replaced with a commuting differential operator whose spectrum does not exhibit the sample plateau. This is the case in one dimension when $Ω = {0, \dots, N - 1}$ , yielding a well-posed eigenvector problem for the Slepian sequences [20, Chapter 2]. For higher dimensions, we may similarly take $Ω$ as a subgrid of $Z^{d}$ , that is, $Ω = {0, \dots, N - 1}^{d}$ . In this case, Slepian tapers are tensor products of one-dimensional Slepian sequences and the resulting estimators achieve good performance [17].

The existence of a commuting differential operator was referred to as a “lucky accident” by Slepian [39]. Similar devices in two dimensions are only available for special cases involving radial symmetries [15, 20, 37] or for polar caps in spherical geometry [36]. For non-symmetric domains, there are no adequate commuting differential operators [8, 16, 29]. Other stable numerical strategies exist only for particular domains [24, 25, 21, 23], or for modified tapers that have analytic expressions [32]. For more general domains, such as the complement of the disk within a rectangular grid (2), no useful symmetries seem to be available. Consequently, the calculation of the Slepian tapers for multitaper estimators on irregular domains is affected by numerical instability. However, as a consequence of the analysis in this article, such instabilities do not preclude the a stable implementation of the multitaper estimator. Indeed, even if the calculation of the Slepian tapers is severely ill-posed, effective numerical proxies can be calculated efficiently (see Section 5).

Multitaper estimators on irregular domains have been studied previously in the literature. Bronez referred to the problem of “irregularly sampling of multidimensional processes” (in reference to the geometry of the set of available samples) and named the corresponding tapers “generalized prolate spheroidal sequences” [10]. More recently, multitaper estimators associated with irregular domains, including the more challenging setting of spherical geometries, have been instrumental in geosciences and climate analysis [20, 18]. In the absence of a commuting operator, practitioners often calculate the tapers by a direct eigenvalue decomposition of $T^{Ω, W}$ , an operation that is admittedly unstable [31, 36, 37, 18], although effective in practice (see Remark 2). Irregular spectra are also relevant in the one-dimensional setting, for example in the field of cognitive radio [19], where the opportunistic occupation of transmission frequencies leads to complex geometries that can be leveraged through carefully designed irregular sampling patterns [12].

4 Analysis of multitaper estimator

Let us now consider the performance of the multitaper estimator ${ˆ S}^{m t}$ . Previous analysis for the one-dimensional case [3] relies on the following aggregated measure for the spectral resolution of the tapers, known as the accumulated spectral window:

ρ (ξ) := \frac{1}{K} K \sum k = 1 {∣ ∣ ∣ ∣ \sum q \in Ω m_{k} [q] e^{2 π i q ξ} ∣ ∣ ∣ ∣}_{k}^{2} ξ \in [- 1 / 2, 1 / 2]^{d} .

(12)

This window determines the smoothness imposed on ${ˆ S}^{m t}$ and therefore controls the bias, variance, and consequently mean squared error (MSE) of the estimator. Indeed, for a multitaper estimator based on Slepian tapers, we have the following estimate, which is a minor extension of a result from Lii and Rosenblatt [26] (the proof is provided in Appendix A).

Proposition 1

Suppose that the spectral density $S$ of $X$ is a $1$ -periodic bounded function. If $ρ$ is defined from the Slepian tapers by (12), the multitaper estimator ${ˆ S}^{m t} (ξ)$ satisfies:

		$Var {{ˆ S}^{m t} (ξ)} := E {{∣ ∣ {ˆ S}^{m t} (ξ) - E {{ˆ S}^{m t} (ξ)} ∣ ∣}^{2}} ≲ \frac{∥ S ∥_{\infty}^{2}}{K},$
		$B i a s {{ˆ S}^{m t} (ξ)} := E {∣ ∣ {ˆ S}^{m t} (ξ) - S (ξ) ∣ ∣} = \| S (ξ) - S * ρ (ξ) \| .$

$_{□}$

Note that the above result can also be proved for an arbitrary set of tapers. Abreu & Romero [3] showed that the spectral window $ρ$ associated with Thomson’s classical (that is, one-dimensional) multitaper estimator resembles a bump function localized in the interval $[- W / 2, W / 2]$ and provided concrete error estimates [3]. Such description of the spectral window, combined with Proposition 1, leads to concrete MSE bounds for the classical multitaper as a function of $K$ , thus elaborating on the more qualitative analysis by Lii and Rosenblatt [26]. These types of estimates are also instrumental in the analysis of multitapering for slowly evolving spectral densities [44].

As a first contribution, we extend the description of the spectral window $ρ$ to arbitrary dimension and general acquisition domains. We let $n_{\partial Ω}$ denote the digital perimeter of $Ω$ :

n_{\partial Ω} = \sum q \in Z^{d} d \sum j = 1 ∣ ∣ 1_{Ω} (q + e_{j}) - 1_{Ω} (q) ∣ ∣,

where ${e_{j} : j = 1, \dots, d}$ is the canonical basis of $Z^{d}$ . In the following, we will assume that $W^{d - 1} n_{\partial Ω} \geq 1$ to avoid degenerate cases.

Theorem 1

Consider the spectral window $ρ$ defined by the Slepian tapers $m_{1}, \dots, m_{K}$ obtained from (6). Assume that $W^{d - 1} n_{\partial Ω} \geq 1$ , where $W$ and $K$ are related by (7). Then

∥ ∥ ρ - \frac{1}{W^{d}} 1_{[- W / 2, W / 2]^{d}} {∥ ∥}_{L^{1} ([- 1 / 2, 1 / 2]^{d})} ≲ \frac{n_{\partial Ω} W^{d - 1}}{K} [1 + log (\frac{n_{Ω}}{n_{\partial Ω}})] .

(13)

$_{□}$

Related results in the context of the short-time Fourier transform can be found in

[1, 2]. The proof of Theorem 1 is postponed to Appendix A.

Combining Proposition 1 and Theorem 1, we obtain MSE bounds for the multitaper estimator on general domains $Ω$ . We also present a simplified expression for the bound, valid in the so-called fine-scale regime. In this regime, one considers a class of acquisition domains $Ω$ for which there is a constant $C > 0$ such that $n_{\partial Ω} ≲ C {n_{Ω}}^{\frac{d - 1}{d}}$ . An instance of this regime occurs, for example, if $Ω$ arises from increasingly fine-scale discretizations of a certain continuous subset of $R^{d}$ , where $C$ depends on the smoothness of the subset.

Theorem 2

Suppose that the spectral density $S$ of $X$ is a $1$ -periodic $C^{2}$ function. Assume that $K n_{\partial Ω}^{d / (d - 1)} \geq n_{Ω}$ . Then the multitaper estimator ${ˆ S}^{m t}$ with $K$ tapers satisfies the following mean squared error bound:

	$MSE {{ˆ S}^{m t} (ξ)}$	$= E {{∣ ∣ {ˆ S}^{m t} (ξ) - S (ξ) ∣ ∣}^{2}}$		(14)
		$≲ ∥ S ∥_{C^{2}}^{2} (\frac{K^{4 / d}}{{n_{Ω}}^{4 / d}} + \frac{n_{\partial Ω}^{2}}{{n_{Ω}}^{2 - 2 / d} K^{2 / d}} {[1 + log (\frac{n_{Ω}}{n_{\partial Ω}})]}^{2} + \frac{1}{K}) .$		(15)

In particular, when $d = 2$ , and, in the fine-scale regime (where $n_{\partial Ω} ≲ \sqrt{n_{Ω}}$ ), the choice $K = ⌈ {n_{Ω}}^{2 / 3} ⌉$ (or, equivalently, $W = {n_{Ω}}^{- 1 / 6}$ ) gives

MSE {{ˆ S}^{m t} (ξ)} ≲ {n_{Ω}}^{- 2 / 3} \cdot {log}^{2} (n_{Ω}) \cdot ∥ S ∥_{C^{2}}^{2} .

(16)

$_{□}$

The proof of Theorem 2 is postponed to Appendix A.

In dimension $d = 1$ , the choice $K = ⌈ {n_{Ω}}^{4 / 5} ⌉$ in (14) (or, equivalently, $W = {n_{Ω}}^{- 1 / 5}$ ) gives

MSE {{ˆ S}^{m t} (ξ)} ≲ {n_{Ω}}^{- 4 / 5} ∥ S ∥_{C^{2}}^{2},

(17)

recovering the result in [3]. In this setting, the minimax risk corresponding to the stronger error measure given by the expected operator norm of the covariance matrix satisfies

inf ˆ S sup S E {sup ξ {∣ ∣ ˆ S (ξ) - S (ξ) ∣ ∣}^{2}} ≍ n^{- 4 / 5} log (n)^{4 / 5} .

(18)

Here, the infimum is taken among all estimators $ˆ S$ based on $n$ consecutive samples, and the first $sup$ is over all spectral densities $S$ with Hölder exponent $2$ satisfying a certain smoothness bound, which determines the implied constant [6, 11]; see also [22]. We are unaware of benchmarks for the spectral estimation problem related to two-dimensional acquisition domains.

The estimates leading to Theorem 1 are also relevant in numerical analysis. For example Proposition 4 below improves on [27].

The above results may be generalized to cover arbitrary frequency profiles for the tapers instead of squares $[- W / 2, W / 2]^{d}$ . For example, in applications where the spectral density does not display any particular anisotropy related to the axes, a disk might be a better choice. For simplicity we do not pursue such generalizations in the present work.

[Shannon number] Theorem 2 and the technical lemmas in Appendix A

provide non-asymptotic bounds for certain heuristic calculations concerning the so-called

Shannon number [37]. These involve the sum of the most significant eigenvalues of the spectral concentration problem for irregular domains, and are typically formulated in the large-scale asymptotic regime (see also [36, Section 7]).

5 Randomized implementation of the multitaper estimator

As discussed in Section 3, calculating the Slepian tapers for arbitrary domains is a difficult task due to the inherent instability of the underlying eigenproblem. To avoid this difficulty, we replace the Slepian tapers with another set of tapers which are easier to calculate but have the same linear span. While these tapers can be very different, it turns out that they yield the same multitaper estimator.

Proposition 2

Let ${g_{1}, \dots, g_{L}}$ and ${{˜ g}_{1}, \dots, {˜ g}_{L}}$ be two orthonormal sets that span the same linear space $V \subseteq ℓ^{2} (Z^{d})$ . Then the corresponding multitaper estimators coincide, that is,

\frac{1}{L} L \sum k = 1 {ˆ S}_{g_{k}} (ξ) = \frac{1}{L} L \sum k = 1 {ˆ S}_{{˜ g}_{k}} (ξ) ξ \in [- 1 / 2, 1 / 2]^{d} .

(19)

$_{□}$

The proof is found in Appendix A.

The above result significantly simplifies our task, since while computing the top $K$ eigenvectors of $T^{Ω, W}$ (i.e., the Slepian tapers) is unstable, this is not the case when computing their linear span. Indeed, there is a large spectral gap between the first $K$ eigenvalues the rest of the spectrum (which Theorem 1 and the estimates in Appendix A validate). As a result, computing the associated subspace is a well-posed problem [41, 14].

We propose to compute these randomized tapers through a block power method applied to $T^{Ω, W}$ . These are then used to compute the multitaper spectral estimator for a given realization of a stationary process. The resulting algorithm is presented as Procedure (P1).

Procedure P1: The randomized multitaper estimator. Data: $∙$ Samples $X [q]$ , for $q \in Ω$ $∙$ A bias-variance trade-off parameter $K$ , $∙$ Number of power method iterations $M$ , $∙$ A target frequency $ξ \in [- 1 / 2, 1 / 2]^{d}$ . Result: ${ˆ S}^{r a n} (ξ)$ , an estimator for $S (ξ)$ . Procedure: [itemindent=0pt] Let $W$ be defined by $K = ⌈ n_{Ω} \times W^{d} ⌉$ .

Draw a random matrix

L \in R^{Ω \times K}

, where

L [q, k] \sim N (0, 1)

, for

q \in Ω

and

k \in {1, \dots, K}

, are drawn independently. Repeat

M

times:

Compute a QR decomposition:

Q \cdot R = T^{Ω, W} \cdot L

. Set

L = Q

. Let

{˜ m}_{1}, \dots, {˜ m}_{K}

be the columns of

L

and

{ˆ S}^{r a n} (ξ) := \frac{1}{K} K \sum k = 1 {ˆ S}_{{˜ m}_{k}} (ξ) .

In exact arithmetic, the column space of $L$ converges to the span of the Slepian tapers ${m_{1}, \dots, m_{K}}$ as we increase $M$ . Consequently, by Proposition 2, ${ˆ S}^{r a n} (ξ)$ converges to ${ˆ S}^{m t} (ξ)$ . Since the gap between the $K$ th and the $(K + 1)$ th eigenvalues of $T^{Ω, W}$ is typically non-negligible [46, 47], convergence is relatively fast and only a moderate number of iterations $M$ is necessary. In fact, a single iteration $M = 1$ is often sufficient for many applications. Furthermore, the only step in Procedure (P1) that depends on the data is the last one. As a result, we may precompute the randomized tapers ${˜ m}_{1}, \dots, {˜ m}_{K}$ ahead of time and then apply them to the data when needed.

The matrix $T^{Ω, W}$ is supported on $Ω \times Ω$ ; in Procedure (P1) it is treated as an element of $R^{Ω \times Ω}$ .

Applying $T^{Ω, W}$ to a vector involves an insertion followed by a convolution and a truncation. Insertion and truncation are diagonal operators, while the convolution is a Toeplitz operator that can be applied using a fast Fourier transform. Multiplication of a vector by $T^{Ω, W}$ is thus achieved in $O (n_{Ω} log n_{Ω})$ time.

[Relation to common practice] The numerical calculation of the Slepian tapers on general domains as eigenvectors of $T^{Ω, W}$ is a severely ill-posed problem. Indeed, the condition number for the calculation of a single eigenvector is inversely proportional to its distance to the rest of the spectrum–see, for example, [33, Equation 3.45] and [35, Section 2.5]–and this spectral gap is small because of the plateau spectral profile of $T^{Ω, W}$ when $K$ is large. However, practitioners resort to such direct methods with remarkable results [31, 36, 37, 18]. Our analysis partly explains this success: unless $K$ is small enough, a standard finite-precision eigenvalue routine will fail to compute the true Slepian tapers, but will still compute an orthonormal basis for their linear span. In other words, the resulting tapers will be linear combinations of the true Slepian tapers, and therefore yield an equally useful multitaper estimator.

6 Numerical results

To empirically evaluate the performance of the proposed randomized multitaper estimator, we perform a few numerical experiments on synthetic data. The randomized multitaper is shown to perform comparably to the standard Slepian multitaper, but with applicability to irregular domains. We also evaluate its accuracy on single-particle cryo-EM data, for both synthetic and experimental images.

6.1 Comparison with tensor Slepian tapers

As discussed in Section 3, tensor products of Slepian functions may be used when $Ω$ is a subgrid of $Z^{d}$ . To illustrate the performance of the randomized multitaper estimator, we therefore first compare it with that of the tensor Slepian multitaper estimator defined on the subgrid shown in Figure 2(a).

(a)

(b)

(c)

Figure 3: (a) A subgrid domain of size

85

-by-

85

within a larger domain of size

128

-by-

128

. (b) The target frequency profile corresponding to

W = 1 / 16

. (c) An irregular domain: the complement of a disk with radius

43

inside a

128

-by-

128

square.

The top row Figure 4 shows a few tensor Slepian functions defined on the subgrid domain of Figure 2(a) with target frequency profile given by 2(b), which corresponds to $W = 1 / 16$ . Below are randomized tapers defined on the same subgrid with the same frequency profile and calculated using Procedure (P1) for $M = 2$ . Both sets of tapers are properly supported on given domain, but their appearance is quite different. Nonetheless, their accumulated spectral windows, shown in Figures 4(a) and 4(b), both closely agree with the target profile in Figure 2(b).

To evaluate the performance of these tapers in a multitaper estimation setting, we generate a Gaussian process with spectral density given by Figure 5(a) on a $128$ -by- $128$ -pixel square. The estimated density of the tensor Slepian multitaper estimator is given in Figure 5(b) while that of the randomized tapers for $M = 2$ iterations is in Figure 5(c). Both agree quite well with the true density. Indeed, the normalized root mean squared error of the tensor Slepian estimator $∥ {ˆ S}^{m t} - S ∥ / ∥ S ∥$ is approximately $2.99 \cdot 10^{- 1}$ while for the randomized tapers, we have $∥ {ˆ S}^{r a n} - S ∥ / ∥ S ∥$ of about $3.00 \cdot 10^{- 1}$ . The deviation between the two estimators is $∥ {ˆ S}^{m t} - {ˆ S}^{r a n} ∥ / ∥ {ˆ S}^{m t} ∥ \approx 1.82 \cdot 10^{- 2}$ . If we increase the number of iterations to $M = 72$ , we obtain equality up to machine precision with $∥ {ˆ S}^{m t} - {ˆ S}^{r a n} ∥ / ∥ {ˆ S}^{m t} ∥ \approx 5.73 \cdot 10^{- 16}$ .

Figure 4: (top row) Three tensor Slepian functions defined on the subgrid domain of Figure 2(a). Three randomized tapers defined on (middle row) the domain of Figure 2(a) and (bottom row) the disk complement domain of Figure 2(c).

(a)

(b)

(c)

Figure 5: The accumulated spectral windows of (a) the tensor Slepian functions, (b) the randomized tapers for the subgrid domain of Figure 2(a), and (c) the randomized tapers for the disk complement domain of Figure 2(c).

(a)

(b)

(c)

(d)

Figure 6: (a) The spectral density of a two-dimensional stochastic process. The density from a

128

-by-

128

realization of the process using (b) tensor Slepian tapers defined on the subgrid of Figure 2(a), (c) the randomized tapers defined on the same subgrid, and (d) the randomized tapers defined on the disk complement of Figure 2(c).

6.2 Illustration for irregular domain

We now replace the subgrid domain with the disk complement domain shown in Figure 2(c). Following the discussion of Section 3, it is computationally intractable to solve the eigenvalue (11) for the desired tapers. We can, however, compute randomized tapers over this domain using Procedure (P1) with $M = 2$ . A few sample tapers are shown in the bottom row Figure 4. Again, their appearance is quite different from the Slepian tapers in the top row, but their accumulated spectral window shown in Figure 4(c) agrees well with the target of Figure 2(b).

Applying these tapers to estimate the spectral density of Figure 5(a) from one realization gives the density depicted in Figure 5(d). The normalized root mean squared error is approximately $2.16 \cdot 10^{- 1}$ .

6.3 Cryo-EM: Synthetic data

(a)

(b)

(c)

(d)

(e)

Figure 7: Simulation of cryo-EM images. (a,b) The clean projection images obtained from the density map of a 70S ribosome. (c) The spectral density

X

of the noise. (d,e) The projection images combined with noise generated using the spectral density.

To evaluate our approach in a real-world application, we consider the estimation of noise power spectra in cryo-EM. In cryo-EM imaging, a solution containing macromolecules of interest are frozen in a thin layer of vitreous ice which is then exposed to an electron beam. A sensor records the transmitted electrons, resulting in a set of tomographic projections depicting the molecules from various viewing angles [13]. To reduce specimen damage, the electron dose is kept low, resulting in exceptionally noisy images, with noise power often exceeding that of the signal by a factor of ten or more. Three-dimensional reconstruction of the molecules, the goal of cryo-EM, therefore requires accurate characterization of the noise model. Since the noise characteristics vary with microscope configuration, ice thickness, and other experimental factors, we cannot rely on ensemble averages for accurate estimation of the noise power spectrum.

To evaluate our approach for the cryo-EM application, we generate a number of projection images using a 70S ribosome density map on a grid with $N = 128$ . These are shown in Figures 6(a) and 6(b). We then generate Gaussian noise with a spectral density $X$ given in Figure 6(c). Adding the noise to the simulated projections, we obtain the images shown in Figure 6(d) and 6(e).

(a)

(b)

Figure 8: Estimation results for cryo-EM. (a) The mask

Ω

for

R = 60

. (b) The inscribed rectangular

Ω^{g r i d}

for the same

R

. (c) The squared bias of the masked periodogram (MPER), tensor multitaper (MT), and randomized multitaper (RT) estimators. (d) The variance of the estimators. (e) The mean squared error (MSE) of the estimators.

Since the central disk of the images contains both the projected molecular density and the noise, we would like to estimate $S$ outside of this disk. Specifically, we would like to restrict our estimator to a set $Ω$ of samples like the one shown in Figure 7(a). A common approach is to define the mask taper $m [q] = {n_{Ω}}^{- 1 / 2} 1_{Ω} (q)$ and use (4) to calculate the tapered periodogram ${ˆ S}_{m} (ξ)$ [45, 7]. As only $K = 1$ taper is used, the result has high variance, so the estimates are averaged over a set of images to obtain an adequate estimate. However, this fails to account for any variability in the noise models of the individual images.

A more accurate estimator is obtained by replacing $Ω$ with a union of rectangular subgrids $Ω^{g r i d} \subset Ω$

Ω^{g r i d} := ⋃ a \in {0, 1}^{2} {q \in Q_{N} : {| q - N a |}_{1} < \frac{N}{2} - \frac{R}{\sqrt{2}}},

(20)

where ${| \cdot |}_{1}$ is the $ℓ^{1}$ -norm. We then apply a standard tensor multitaper estimator to each subgrid and average the results. These tapers are tensor products of one-dimensional Slepian sequences, as described in Section 3. The $Ω^{g r i d}$ corresponding to the $Ω$ of Figure 7(a) is shown in Figure 7(b). Depending on the geometry of $Ω$ , $Ω^{g r i d}$ may discard many points in $Ω$ , increasing variance of the estimate. For small $K$ , this also increases bias, as fewer points in this regime leads to a wider accumulated spectral window $ρ$ .

Performance. We now compare the performance of these baseline estimators, the tapered periodogram ${ˆ S}_{m}$ and the tensor multitaper estimator ${ˆ S}^{m t}$ , to that of our proposed estimator ${ˆ S}^{r a n}$ . For a given $R$ , we define $Ω$ by (2) and $Ω^{g r i d}$ by (20). When $R$ is low, we expect a large bias as the samples are contaminated by the projected density maps, while increasing $R$ results in lower bias but higher variance due to the lower number of available samples. The bias, variance, and mean squared error of the three estimators as a function of $R$ are shown in Figures 7(c), 7(d), and 7(e), respectively.

The bias of the randomized multitaper estimator is higher than the tensor multitaper for low $R$ . This is due to the original mask $Ω$ overlapping with the support of the projection images to a greater extent than $Ω^{g r i d}$ at these radii. As $R$ increases above $60$ , however, the bias for both estimators drops down to the same level. Since it has $K = 1$ , no smoothness is imposed on the tapered periodogram which therefore has lower bias compared to the multitaper estimators.

At low $R$ , the influence of the clean projection images yields high variance for all the estimators, since the images vary according to viewing angle. The effect is exacerbated for the randomized multitaper estimators and the tapered periodogram since $Ω$ has greater overlap with the projections compared to $Ω^{g r i d}$ . As $R$ increases, however, the randomized multitaper estimator ${ˆ S}^{r a n}$ enjoys a lower variance, since it draws upon a larger number of samples compared to the tensor multitaper estimator ${ˆ S}^{m t}$ , reducing the variance by a factor of two. The tapered periodogram, meanwhile, has high variance for all $R$ since it only employs a single taper, providing no variance reduction.

The low bias and variance for high $R$ combine to yield a lower MSE for the randomized multitaper estimator ${ˆ S}^{r a n}$ compared to the other estimators. On average, the randomized multitaper estimator gives an error a factor of two lower than the tensor multitaper, and beats the tapered periodogram by a factor of three.

(a)

(b)

Figure 9: Estimation results for experimental cryo-EM. (a,b) Two sample images from the EMPIAR-10028 dataset containing projection images of an 80S ribosome. (c) The mean squared errors (MSEs) of the masked periodogram (MPER), tensor multitaper (MT), and randomized multitaper (RT) estimators.

6.4 Cryo-EM: Experimental data

We now evaluate performance on images consisting of experimental projections of an 80S ribosome complex from the EMPIAR-10028 dataset [43]. The images are defined on a grid with $N = 360$ . For our evaluation, we choose a subset of the first $1024$ images from the dataset. Two sample images from this subset are shown in Figure 8(a) and 8(b).

To obtain a reasonable approximation of the true power spectrum for each projection image, we process the entire dataset through the RELION software package [34]. This yields estimates of the underlying molecular density and allows us to simulate clean projection images for each of the noisy images. By subtracting the estimated clean images from the noisy images, we obtain a set of images consisting mostly of noise. We may then estimate their power spectra by applying a tensor multitaper estimator over the entire grid $Q_{N}$ . Since these estimates incorporate the noise in the center of the image, they should provide an accurate estimate of the spectral density of the noise in the whole image. We shall therefore use these as a standard against which we compare the estimates obtained from the complement of the disk on the original projection images.

As before, we calculate the tapered periodogram ${ˆ S}_{m}$ on $Ω$ , the tensor multitaper estimator ${ˆ S}^{m t}$ on $Ω^{g r i d}$ , and the randomized multitaper estimator ${ˆ S}^{r a n}$ on $Ω$ . Each estimated power spectrum is compared to the “ground truth” power spectrum obtained for the corresponding image as described above. The resulting MSE’s are plotted for each estimator in Figure 8(c).

As before, the tapered periodogram performs badly since it has $K = 1$ , resulting in high variance and therefore high MSE. The tensor multitaper estimator performs well for low $R$ , but increasing $R$ reduces the area of $Ω^{g r i d}$ , yielding higher variance and higher MSE. A similar behavior is observed for the randomized multitaper estimates, but the reduction in area for $Ω$ is not as drastic for increasing $R$ , so the variance remains small compared to the tensor multitaper estimator. At $R \approx 110$ , the error is minimized and the randomized multitaper outperforms the tensor multitaper by a factor of $1.5$ .

7 Conclusion

We have analyzed the multitaper estimator on arbitrary acquisition domains, providing performance bounds on the mean squared error. Furthermore, we use the fact that the multitaper estimation only depends on the linear span of the tapers to introduce a new implementation based on randomized tapers. These span the same linear space as the Slepian tapers, but are easier to calculate since they do not involve calculating eigenvectors associated with closely spaced eigenvalues. The performance of the randomized multitaper estimator is shown to be comparable to that using Slepian tapers when these are available for rectangular domains. For more general domains, we demonstrate the improved performance of the randomized multitaper estimator compared to Slepian multitaper estimators defined on rectangular subgrids. This is shown for both synthetic examples and on experimental data obtained from cryo-EM imaging.

Future directions of research involve adapting the factor analysis framework proposed for periodogram estimators [4] to multitaper estimators. This would allow for greater variance reduction when a linear structure exists in the variability of spectral density between independent realizations of the point process. Such a situation arises, for example, in the cryo-EM noise estimation task, where a set of underlying noise sources combine at arbitrary strengths to yield the noise process in a given image.

Acknowledgment

The authors are very grateful to Luís Daniel Abreu, Tomasz Hrycak, and Amit Singer, who motivated this article and provided valuable input.

Appendix A Proofs

a.1 Convolution estimates

First, we note that

∥ f * g ∥_{\infty} \leq ∥ f ∥_{\infty} ∥ g ∥_{L^{1} ([- 1 / 2, 1 / 2]^{d})} .

(21)

Furthermore, for a parameter $W \in (0, 1 / 2)$

we consider the periodic extension of the normalized characteristic function

W^{- d} 1_{[- W / 2, W / 2]^{d}}

, and, by a slight abuse of notation, we define its convolution with a 1-periodic function

f

f * W^{- d} 1_{[- W / 2, W / 2]^{d}} (x) = \frac{1}{W^{d}} \int_{[- W / 2, W / 2]^{d}} f (x - y) d y .

An estimate on a second-order Taylor expansion shows that

∥ f - \frac{1}{W^{d}} f * 1_{[- W / 2, W / 2]^{d}} ∥_{\infty} ≲ ∥ f ∥_{C^{2}} W^{2} .

(22)

a.2 Trace and norm of the Toeplitz operators

The estimates derived in this section are also relevant in numerical analysis of Fourier extensions. In particular, they improve on the results in [27] by avoiding assumptions on the (digital) topology of the set $Ω$ . See also [28] for related estimates.

Proposition 3

Let $a, b \in ℓ^{1} (Z^{d})$ with $\sum_{q} a [q] = 1$ . Then

∥ a * b - b ∥_{ℓ^{1} (Z^{d})} \leq ∥ \nabla b ∥_{1} \sum q {| q |}_{\infty} | a (q) |,

where $(\nabla_{k} b) [q] := b [q + e_{k}] - b [q]$ for $k = 1, \dots, d$ , $∥ \nabla b ∥_{1} = \sum_{k = 1}^{d} ∥ \nabla_{k} b ∥_{1}$ , and ${| q |}_{\infty} = max {| q_{1} |, \dots, | q_{d} |}$ . $_{□}$

Proof

For $q \in Z^{d}$ and $k \in {1, \dots, d}$ we let $χ_{k} q \in Z^{d}$ be the truncated vector $(χ_{k} q)_{j} = q_{j} 1_{j \leq k}$ , and also $χ_{0} q = 0$ . For $q^{'} \in Z^{d}$ , let us write

	$(a * b) [q^{'}] - b [q^{'}] = \sum q \in Z^{d} (b [q + q^{'}] - b [q^{'}]) a [- q]$
	$= \sum q \in Z^{d} d \sum k = 1 (b [χ_{k} q + q^{'}] - b [χ_{k - 1} q + q^{'}]) a [- q]$
	$= \sum q \in Z^{d} d \sum k = 1 (b [χ_{k - 1} q + q^{'} + q_{k} e_{k}] - b [χ_{k - 1} q + q^{'}]) a [- q]$
	$= \sum q \in Z^{d} d \sum k = 1 \| q_{k} \| - 1 \sum l = 0 (b [χ_{k - 1} q + q^{'} + sgn (q_{k}) (l + 1) e_{k}] - b [χ_{k - 1} q + q^{'} + sgn (q_{k}) l e_{k}]) a [- q]$
	$= \sum q \in Z^{d} d \sum k = 1 \| q_{k} \| - 1 \sum l = 0 \nabla_{k} b [χ_{k - 1} q + q^{'} + sgn (q_{k}) l e_{k}] a [- q] .$

Therefore,

	$∥ a * b - b ∥_{1}$	$\leq \sum q \in Z^{d} d \sum k = 1 \| q_{k} \| - 1 \sum l = 0 \sum q^{'} \in Z^{d} ∣ ∣ \nabla_{k} b [χ_{k - 1} q + q^{'} + sgn (q_{k}) l e_{k}] ∣ ∣ \| a [- q] \|$
		$\leq \sum q \in Z^{d} d \sum k = 1 ∥ \nabla_{k} b ∥_{1} \| q_{k} \| \| a [- q] \| \leq ∥ \nabla b ∥_{1} \sum q \in Z^{d} {\| q \|}_{\infty} \| a [q] \|,$

as desired, since $| q_{k} | \leq {| q |}_{\infty}$ . (See [28, 1] for related estimates.) $_{■}$

Proposition 4

Let $T^{Ω, W}$ be the matrix in (9). Then

trace [T^{Ω, W}] - t r a c e [(T^{Ω, W})^{2}] ≲ n_{\partial Ω} W^{d - 1} [1 + log (\frac{n_{Ω}}{n_{\partial Ω}})] .

(23)

$_{□}$

Proof

Step 1. (Computations). The function

h [q] := d \prod j = 1 \frac{sin (π W q_{j})}{π q_{j}} = W^{d} sinc (W q), q = (q_{1}, \dots, q_{d}) \in Z^{d},

(24)

determines the Fourier series

Note that, in terms of $h$ , (9) reads: $T_{q, q^{'}}^{Ω, W} = 1_{Ω} (q) h [q - q^{'}] 1_{Ω} (q^{'})$ . We first compute

	$trace [T^{Ω, W}]$	$= \sum q \in Ω h [q - q] = n_{Ω} h_{0} = n_{Ω} ∥ 1_{[- W / 2, W / 2]^{d}} ∥_{1}$		(25)
		$= W^{d} n_{Ω} = W^{d} \sum q \in Z^{d} 1_{Ω} [q] .$		(26)

Second,

	$trace [(T^{Ω, W})^{2}]$	$= \sum q, q^{'} \in Z^{d} 1_{Ω} [q] {∣ ∣ h [q - q^{'}] ∣ ∣}^{2} 1_{Ω} [q^{'}]$		(27)
		$= W^{d} \sum q \in Z^{d} (1_{Ω} * a) [q] 1_{Ω} [q],$		(28)

where

a [q] := W^{- d} {| h [q] |}^{2} .

Step 2. (Truncation errors). Let $L > 0$ and consider the function

~ a [q] := a [q] 1_{| q | \leq L} .

We claim that, for $L \geq 2$ ,

		$∥ a - ~ a ∥_{1} = W^{- d} \sum q \in Z^{d}, \| q \| > L {\| h [q] \|}^{2} ≲ (W L)^{- 1},$		(29)
		$\sum q \in Z^{d} \| q \| \| ~ a [q] \| = W^{- d} \sum q \in Z^{d}, \| q \| \leq L \| q \| {\| h [q] \|}^{2} ≲ \frac{log L}{W} .$		(30)

To show these estimates, we write $h [q] = r [q_{1}] \dots r [q_{d}]$ with

r [q_{1}] = \frac{sin (π W q_{1})}{π q_{1}} = W sinc (W q_{1}), q_{1} \in Z .

We first note the following:


	$\sum q_{1} \in Z, \| q \| > L {\| r [q_{1}] \|}_{1}^{2} ≲ \sum q_{1} \in Z, \| q_{1} \| > L \frac{1}{{\| q_{1} \|}_{1}^{2}} ≲ \frac{1}{L},$
	$\sum q_{1} \in Z, \| q_{1} \| \leq L \| q_{1} \| {\| r [q_{1}] \|}_{1}^{2} ≲ \sum q_{1} \in Z \frac{1}{\| q_{1} \|} ≲ log (L), L$

Multitaper estimation on arbitrary domains

Optimal mean squared error bandwidth for spectral variance estimators in MCMC simulations

How many samples are needed to reliably approximate the best linear estimator for a linear inverse problem?

A recursive procedure for density estimation on the binary hypercube

POT-flavored estimator of Pickands dependence function

Minimum Mean-Squared-Error Autocorrelation Processing in Coprime Arrays

Online Balanced Experimental Design

Hyperspectral Pigment Analysis of Cultural Heritage Artifacts Using the Opaque Form of Kubelka-Munk Theory

1 Introduction

2 Spectral estimation on irregular domains

3 Multitaper estimators

4 Analysis of multitaper estimator

Proposition 1

Theorem 1

Theorem 2

5 Randomized implementation of the multitaper estimator

Proposition 2

6 Numerical results