Multitaper estimation on arbitrary domains

TL;DR

This paper introduces a new approach for multitaper spectral density estimation on arbitrary domains, providing theoretical performance bounds and demonstrating improved accuracy in cryo-electron microscopy applications.

Contribution

It shows that multitaper estimators depend only on the span of tapers, allowing the use of proxy tapers computed with standard algorithms, and establishes performance bounds for arbitrary domains.

Findings

Reduces mean squared error by over 50% in cryo-EM data

Validates use of proxy tapers for arbitrary domains

Provides theoretical bounds for estimator performance

Abstract

Multitaper estimators have enjoyed significant success in estimating spectral densities from finite samples using as tapers Slepian functions defined on the acquisition domain. Unfortunately, the numerical calculation of these Slepian tapers is only tractable for certain symmetric domains, such as rectangles or disks. In addition, no performance bounds are currently available for the mean squared error of the spectral density estimate. This situation is inadequate for applications such as cryo-electron microscopy, where noise models must be estimated from irregular domains with small sample sizes. We show that the multitaper estimator only depends on the linear space spanned by the tapers. As a result, Slepian tapers may be replaced by proxy tapers spanning the same subspace (validating the common practice of using partially converged solutions to the Slepian eigenproblem as tapers). These proxies may consequently be calculated using standard numerical algorithms for block diagonalization. We also prove a set of performance bounds for multitaper estimators on arbitrary domains. The method is demonstrated on synthetic and experimental datasets from cryo-electron microscopy, where it reduces mean squared error by a factor of two or more compared to traditional methods.

Figures

Equations

$$\Vert f{\Vert }_{C^2}:=\max \left\{\underset{x\in {\mathbb{R}}^d}{\sup }\left|f(x)\right|,\underset{x\in {\mathbb{R}}^d}{\sup }\left|{\partial }_{x_j}f(x)\right|,\underset{x\in {\mathbb{R}}^d}{\sup }\left|{\partial }_{x_j}{\partial }_{x_{j^{\prime }}}f(x)\right|:j,j^{\prime }=1,\dots ,d\right\},$$

$$1_{\Omega}(x):=\left\{\begin{array}[]{ll}1,&\mbox{if}~{}x\in{\Omega}\\ 0,&\mbox{otherwise.}\end{array}\right.$$

$$\mathrm{Cov}[q,q^{\prime }]=\mathbb{E}\left\{\mathcal{X}[q]\mathcal{X}[q^{\prime }]\right\}q,q^{\prime }\in {\mathbb{Z}}^d.$$

$$\mathrm{Cov}[q,q^{\prime }]=r[q-q^{\prime }]q,q^{\prime }\in {\mathbb{Z}}^d.$$

$$S(\xi )=\sum _{q\in {\mathbb{Z}}^d}r[q]{\mathrm{e}}^{2\pi \mathrm{i}⟨q,\xi ⟩},\xi \in [-1/2,1/2{]}^d.$$

$${\mathcal{Q}}_N=\{0,\dots ,N-1{\}}^2$$

$$\Omega :=\left\{q\in {\mathcal{Q}}_N:\left|q-(N/2,N/2)\right|>R\right\},$$

$$m[q]=0,\text{mbox}ifq\notin \Omega .$$

$${\hat{S}}_m(\xi ):={\left|\sum _{q\in \Omega }m[q]\mathcal{X}[q]{\mathrm{e}}^{2\pi \mathrm{i}⟨q,\xi ⟩}\right|}^2,\xi \in [-1/2,1/2{]}^d.$$

$${\hat{S}}^{\mathrm{mt}}(\xi ):=\frac{1}{K}\sum _{k=0}^{K-1}{\hat{S}}_{m_k}(\xi )\xi \in [-1/2,1/2{]}^d.$$

$$\text{mbox}maximize{\int }_{[-W/2,W/2{]}^d}{\left|\sum _{q\in \Omega }m[q]{\mathrm{e}}^{2\pi \mathrm{i}⟨q,\xi ⟩}\right|}^{2}d\xi ,$$

$$\text{mbox}subjectto:\sum _{q\in \Omega }{\left|m[q]\right|}^2=1,\text{mbox}and\mathrm{supp}(m)\subseteq \Omega ,$$

$$W=(K/n_{\Omega }{)}^{1/d},$$

$$K=\lceil n_{\Omega }\times W^d\rceil ,$$

$$\text{mbox}maximize{\int }_{[-W/2,W/2{]}^2}{\left|\sum _{(q_1,q_2)\in \Omega }m[q_1,q_2]{\mathrm{e}}^{2\pi \mathrm{i}(q_1{\xi }_1+q_2{\xi }_2)}\right|}^{2}d{\xi }_{1}d{\xi }_2,$$

$$\text{mbox}subjectto:\sum _{(q_1,q_2)\in \Omega }{\left|m[q_1,q_2]\right|}^2=1,\text{mbox}and\mathrm{supp}(m)\subseteq \Omega ,$$

$$\displaystyle T^{\Omega,W}[q,q^{\prime}]=\left\{\begin{array}[]{ll}W^{d}\,\operatorname{sinc}\big{(}W(q-q^{\prime})\big{)},&\mbox{if}~{}q,q^{\prime}\in{\Omega},\\ 0,&\mbox{otherwise},\end{array}\right.$$

$${\mathrm{sinc}}_d(u):=\prod _{k=1}^d\frac{\sin (\pi u_k)}{\pi u_k},u=(u_1,\dots ,u_d)\in {\mathbb{Z}}^d.$$

$$T^{\Omega ,W}{m}_k={\lambda }_k{m}_k,k\in \{0,\dots ,K-1\},$$

$$\sum _{(q_1^{\prime },q_2^{\prime })\in \Omega }\frac{\sin (\pi W(q_1-q_1^{\prime }))\sin (\pi W(q_2-q_2^{\prime }))}{{\pi }^2(q_1-q_1^{\prime })(q_2-q_2^{\prime })}{m}_k[q_1^{\prime },q_2^{\prime }]={\lambda }_k{m}_k[q_1,q_2],(q_1,q_2)\in \Omega ,$$

$$\rho (\xi ):=\frac{1}{K}\sum _{k=0}^{K-1}{\left|\sum _{q\in \Omega }{m}_k[q]{\mathrm{e}}^{2\pi \mathrm{i}⟨q,\xi ⟩}\right|}^2\xi \in [-1/2,1/2{]}^d.$$

$$\mathrm{Var}\left\{{\hat{S}}^{\mathrm{mt}}(\xi )\right\}:=\mathbb{E}\left\{{\left|{\hat{S}}^{\mathrm{mt}}(\xi )-\mathbb{E}\left\{{\hat{S}}^{\mathrm{mt}}(\xi )\right\}\right|}^2\right\}\le \frac{C}{K}\underset{\xi \in [-1/2,1/2{]}^d}{\sup }S(\xi {)}^2$$

$$\mathrm{Bias}\left\{{\hat{S}}^{\mathrm{mt}}(\xi )\right\}:=\mathbb{E}\left\{\left|{\hat{S}}^{\mathrm{mt}}(\xi )-S(\xi )\right|\right\}=\left|S(\xi )-S\ast \rho (\xi )\right|.$$

$$n_{\partial \Omega }=\sum _{q\in {\mathbb{Z}}^d}\sum _{j=1}^d\left|1_{\Omega }(q+e_j)-1_{\Omega }(q)\right|,$$

$${\int }_{[-1/2,1/2{]}^d}\left|\rho (\xi )-\frac{1}{W^d}{1}_{[-W/2,W/2{]}^d}(\xi )\right|d\xi \le C\frac{n_{\partial \Omega }{W}^{d-1}}{K}\left[1+\log \left(\frac{n_{\Omega }}{n_{\partial \Omega }}\right)\right],$$

$$\mathrm{MSE}\left\{{\hat{S}}^{\mathrm{mt}}(\xi )\right\}$$

$$\le C\Vert S{\Vert }_{C^2}^2\left(\frac{K^{4/d}}{{n_{\Omega }}^{4/d}}+\frac{n_{\partial \Omega }^2}{{n_{\Omega }}^{2-2/d}{K}^{2/d}}{\left[1+\log \left(\frac{n_{\Omega }}{n_{\partial \Omega }}\right)\right]}^2+\frac{1}{K}\right).$$

$$\mathrm{MSE}\left\{{\hat{S}}^{\mathrm{mt}}(\xi )\right\}\le C{n_{\Omega }}^{-2/3}\cdot {\log }^2\left(n_{\Omega }\right)\cdot \Vert S{\Vert }_{C^2}^2.$$

$$\mathrm{MSE}\left\{{\hat{S}}^{\mathrm{mt}}(\xi )\right\}\le C{n_{\Omega }}^{-4/5}\Vert S{\Vert }_{C^2}^2,$$

$$A{n}^{-4/5}\log (n{)}^{4/5}\le \underset{\hat{S}}{\inf }\underset{S}{\sup }\mathbb{E}\left\{\underset{\xi }{\sup }{\left|\hat{S}(\xi )-S(\xi )\right|}^2\right\}\le B{n}^{-4/5}\log (n{)}^{4/5},$$

$$\frac{1}{L}\sum _{k=0}^{L-1}{\hat{S}}_{g_k}(\xi )=\frac{1}{L}\sum _{k=0}^{L-1}{\hat{S}}_{{\tilde{g}}_k}(\xi )\xi \in [-1/2,1/2{]}^d.$$

$${\hat{S}}^{\mathrm{pmt}}(\xi ):=\frac{1}{K}\sum _{k=0}^{K-1}{\hat{S}}_{{\tilde{m}}_k}(\xi ).$$

Full text

{AMS}

62G05, 62M15, 62M40, 33E10, 15B05, 15A18, 42B10, 94A12, 65F15, 92C55

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 $\mathcal{X}$ over $\mathbb{Z}^{d}$ is characterized by its spectral density $S$, defined as the Fourier series of its autocovariance function [36]. Spectral estimation is the task of recovering $S$ given one or more realizations of $\mathcal{X}$ over a finite subset ${\Omega}$ of $\mathbb{Z}^{d}$. The challenge in the estimation problem is two-fold: (i) stochastic fluctuations in $\mathcal{X}$ introduce variance into estimates of $S$, and (ii) spatial constraints restrict us to using only samples from the domain ${\Omega}$. When many realizations are available, the first point may be addressed through ensemble averaging. For several applications, however, such as geosciences, cosmology, and electron microscopy 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 ${\Omega}$ into disjoint subsets, compute spectral density estimates for each, and average the result. In one dimension, with ${\Omega}=\{0,\ldots,N-1\}$ for some $N>0$, this is typically done by partitioning ${\Omega}$ 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) [7]. 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 [53], refines this method by instead computing periodograms over all of ${\Omega}$, but with the samples first multiplied by a set of $K$ sequences called tapers. A favorable trade-off between bias and variance is obtained for specific sequences known as discrete prolate spheroidal sequences, or Slepian tapers [50].

Originally defined for $d=1$, multitaper estimators naturally generalize to higher dimensions [48]. The Slepian tapers, however, are defined as solutions to an ill-posed eigenvalue problem. Special geometries, including rectangular grids [24, Chapter 2] and domains involving other symmetries [19, 24, 46, 44], may be analyzed by means of commuting differential operators, yielding well-posed problems that allow for explicit calculation of the tapers. There are other stable algorithms for calculating Slepian tapers, but these are defined only for specific domains [30, 31, 25, 29] or for certain modified tapers [39, 26, 18, 45].

At times, however, more general acquisition domains are needed, such as the complement of a disk (see Figure 1). 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 ${\Omega}$ of the earth [47, 44, 43, 46, 38]. For these general domains ${\Omega}$, all available methods inherit the instability of the underlying eigenvalue problem [46]. That being said, remarkable results have been obtained by direct application of standard eigenvalue solvers and using the resulting eigenvectors as tapers, despite these not being true Slepian tapers [38, 44, 46, 22]. The behavior and performance of these pseudo-multitaper estimators, however, is not well understood.

An alternative solution is to partition this irregular domain into smaller rectangular regions and apply the tensor Slepian multitaper estimator to those. Since each subdomain is a rectangle, Slepian functions are readily calculated on these domains. This strategy, however, lacks the simplicity of 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. This explains the success of standard eigendecomposition algorithms applied to the ill-posed eigenvalue problem since a typical failure mode results in vectors with the same span as the true eigenvectors. To ensure this fortuitous behavior, we propose replacing the standard eigendecomposition with a block eigendecomposition, extracting a basis for the subspace spanned by the leading $K$ eigenvectors. 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 performance bounds for the one-dimensional case [3]. We validate these bounds numerically, showing how they correctly predict the error decay as a function of the acquisition geometry.

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 proxy tapers. Finally, Section 6 illustrates the performance of the proposed estimator using numerical experiments.111Python code for reproducing the results of this work may be found at https://github.com/janden/pmte

Notation

In the following, $:=$ specifies the definition of a quantity. For a vector $x\in\mathbb{R}^{d}$, we let \left|x\right|:=\big{(}\sum_{k=1}^{d}\left|x_{k}\right|^{2}\big{)}^{1/2} be its Euclidean norm. The $C^{2}$-norm of a twice continuously differentiable function $f:\mathbb{R}^{d}\to\mathbb{C}$ is

[TABLE]

whenever finite. We say that $f$ is 1-periodic if $f(x+k)=f(x)$ for all $k\in\mathbb{Z}^{d}$. The indicator function $1_{\Omega}:S\to\mathbb{R}$ of a subset ${\Omega}\subset S$ is defined as

[TABLE]

The inner product between two vectors $q$ and $q^{\prime}$ in $\mathbb{R}^{d}$ is given by $\left<u,v^{\prime}\right>=\sum_{j=1}^{d}u_{j}v^{\prime}_{j}$. Finally, $\mathbb{E}(X)$ denotes the expected value of the random variable $X$ and $\lceil x\rceil$ denotes the smallest integer $k$ such that $k\geq x$.

2 Spectral estimation on irregular domains

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

[TABLE]

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

[TABLE]

The covariance information is equivalently encoded in the spectral density

[TABLE]

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

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

[TABLE]

of pixels, and $\mathcal{X}[q]$ represents additive noise at the pixel $q\in{\mathcal{Q}_{N}}$. A good 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 [16]. To reduce specimen damage, electron dose is kept low, resulting in exceptionally noisy images, as shown in Figure 2. 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

[TABLE]

as illustrated in Figure 1.

Another instance is found in geosciences, where $\mathcal{X}$ is a physical quantity on an approximately flat portion of the earth, and measurements are only available on a subregion ${\Omega}$, corresponding, for example, to a continent [46]. A more sophisticated model replaces the Cartesian grid for a grid on the sphere [44, 38].

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 a single-shot estimator of S. Instead of averaging estimates from several realizations of $\mathcal{X}$, a single-shot estimator only takes a single realization as input. The total number of samples is therefore ${n_{\Omega}}$. This low sample size often results in high variance, so we need a method which regularizes this estimate. The aim of the present work is to study the performance of one such method, the multitaper estimator, and provide an explicit implementation.

3 Multitaper estimators

The inherent variability of $\mathcal{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 [11, 36].

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

[TABLE]

We now define the tapered periodogram of $\mathcal{X}$ with taper $m$ as

[TABLE]

Given an orthonormal family $m_{0},\ldots,m_{K-1}\in\ell^{2}(\mathbb{Z}^{d})$ of tapers supported on ${\Omega}$, we define the corresponding multitaper estimator as the average

[TABLE]

The error of $\widehat{S}^{\mathrm{mt}}$ in estimating $S$ depends on the choice of tapers $m_{0},\ldots,m_{K-1}$.

The tapers of the traditional multitaper estimator are defined using the following spectral concentration problem:

[TABLE]

where

[TABLE]

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

[TABLE]

and controls the amount of smoothness imposed on the estimate $\widehat{S}^{\mathrm{mt}}$. More concretely, the spectral concentration problem (6) for $d=2$ reduces to

[TABLE]

with $W=\sqrt{K/{n_{\Omega}}}$.

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

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

[TABLE]

where $q,q^{\prime}\in\mathbb{Z}^{d}$ and $\operatorname{sinc}$ is the $d$-dimensional normalized sinc function

[TABLE]

The matrix $T^{\Omega,W}$ is indexed by vectors $q,q^{\prime}\in\mathbb{Z}^{d}$, which means that it maps the space of $d$-dimensional arrays to itself. This is done by identifying a $d$-dimensional array with a one-dimensional array through lexicographical ordering of indices. The Slepian tapers $m_{0},\ldots,m_{K-1}$ are thus the top $K$ eigenvectors of $T^{\Omega,W}$, satisfying

[TABLE]

where $\lambda_{0}\geq\lambda_{2}\geq\ldots\geq\lambda_{K-1}$. For $d=2$, we may substitute (12) and (13) into (14), obtaining

[TABLE]

and $m_{k}[q_{1},q_{2}]=0$ for $(q_{1},q_{2})\notin{\Omega}$.

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

Under certain circumstances, $T^{\Omega,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 ${\Omega}=\{0,\ldots,N-1\}$, yielding a well-posed eigenvector problem for the Slepian sequences [24, Chapter 2]. For $d>1$, we may similarly take ${\Omega}$ as a subgrid of ${{\mathbb{Z}}^{d}}$, that is, ${\Omega}=\{0,\ldots,N-1\}^{d}$. In this case, Slepian tapers are tensor products of one-dimensional Slepian sequences [21].

The existence of a commuting differential operator was referred to as a “lucky accident” by Slepian [49]. Similar devices in two dimensions are only available for special cases involving radial symmetries [19, 24, 46] or for polar caps in spherical geometry [44]. For non-symmetric domains, there are no adequate commuting differential operators [10, 20, 35]. Other stable numerical strategies exist only for particular domains [30, 31, 25, 29] or for modified tapers that have analytic expressions [39]. For more general domains 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. As a consequence of the analysis in this work, however, such instabilities do not preclude a stable implementation of the multitaper estimator. Indeed, even if the calculation of the Slepian tapers is severely ill-posed, effective numerical proxies are available (see Section 5).

Multitaper estimators on irregular domains have been studied previously, notably by Bronez, who referred to the problem of “irregularly sampling of multidimensional processes” (referring to the geometry of the set of available samples) and named the corresponding tapers “generalized prolate spheroidal sequences” [12]. More recently, multitaper estimators associated with irregular domains, including the more challenging setting of spherical geometries, have been instrumental in geosciences and climate analysis [24, 22]. Irregular spectra are also relevant in the one-dimensional setting, such as in the field of cognitive radio [23], where the opportunistic occupation of transmission frequencies leads to complex geometries that can be leveraged through carefully designed irregular sampling patterns [14]. In the absence of a commuting operator, practitioners often calculate the tapers by a direct eigenvalue decomposition of $T^{\Omega,W}$, an operation that is admittedly unstable but effective in practice [38, 44, 46, 22]. A potential reason is that, even when the eigendecomposition fails, it may still yield vectors with the same span as the eigenvectors, which, as shown below in Proposition 5.1, perform the same for multitaper estimation (see Section 5).

4 Analysis of multitaper estimator

Let us now consider the performance of the multitaper estimator $\widehat{S}^{\mathrm{mt}}$. 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:

[TABLE]

This window determines the smoothness imposed on $\widehat{S}^{\mathrm{mt}}$ 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 [32] (the proof is provided in Appendix A).

Proposition 4.1**.**

Suppose that the spectral density $S$ of $\mathcal{X}$ is a $1$-periodic bounded function. If $\rho$ is defined from the Slepian tapers by (16), the multitaper estimator $\widehat{S}^{\mathrm{mt}}(\xi)$ satisfies

[TABLE]

for some constant $C\geq 0$, and

[TABLE]

Note that the above result can also be proved for an arbitrary set of tapers. In [3] it was shown that the spectral window $\rho$ 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 4.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 [32]. These types of estimates are also instrumental in the analysis of multitapering for slowly evolving spectral densities [56].

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

[TABLE]

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

Theorem 4.2**.**

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

[TABLE]

*for some constant $C\geq 0$. *

Related results in the context of the short-time Fourier transform can be found in [1, 2]. The proof of Theorem 4.2 is postponed to Appendix A.

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

Theorem 4.3**.**

Suppose that the spectral density $S$ of $\mathcal{X}$ is a $1$-periodic $C^{2}$ function. Assume that $Kn_{\partial\Omega}^{d/(d-1)}\geq{n_{\Omega}}$. Then there exists a constant $C\geq 0$, such that the multitaper estimator $\widehat{S}^{\mathrm{mt}}$ with $K$ tapers satisfies the mean squared error bound

[TABLE]

In particular, when $d=2$, and, in the fine-scale regime (where $n_{\partial\Omega}$ is of the order $\sqrt{{n_{\Omega}}}$), the choice $K=\lceil{n_{\Omega}}^{2/3}\rceil$ (or, equivalently, $W={n_{\Omega}}^{-1/6}$) gives

[TABLE]

The new constant in (20) depends on the constant comparing $n_{\partial\Omega}$ and $\sqrt{{n_{\Omega}}}$. The proof of Theorem 4.3 is postponed to Appendix A.

Remark 1**.**

For $d=1$, the choice $K=\lceil{n_{\Omega}}^{4/5}\rceil$ in (18) (or, equivalently, $W={n_{\Omega}}^{-1/5}$) gives

[TABLE]

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

[TABLE]

*for some $A,B>0$. Here, the infimum is taken among all estimators $\widehat{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 constants $A,B$ [8, 13]; see also [27]. We are unaware of benchmarks for the spectral estimation problem related to two-dimensional acquisition domains. (Expressions for the MSE depending on the signal and noise power are however available, see, e.g. [43, 15, 37].) *

Remark 2**.**

*Theorem 4.3 follows by combining Proposition 4.1 and Theorem 4.2 and thus provides individual estimates for the bias and variance of the multitaper estimator. A small variation of the proofs yields similar results for eigenvalue weighted estimators – see [3, Theorem 2.2]. *

Remark 3**.**

*The estimates leading to Theorem 4.2 are also relevant in numerical analysis. For example Proposition A.3 below improves on [33]. *

Remark 4**.**

*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. *

Remark 5** (Shannon number).**

*Theorem 4.3 and the technical lemmas in Appendix A provide non-asymptotic bounds for certain heuristic calculations concerning the so-called Shannon number [46]. 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 [44, Section 7]). *

5 Proxy Slepian tapers

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. 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, [40, Equation 3.45] and [42, Section 2.5]—and this spectral gap is small because of the plateau spectral profile of $T^{\Omega,W}$ for large $K$. However, practitioners resort to such direct methods with remarkable results [38, 44, 46, 22]. One possible explanation is that, unless $K$ is small enough, a standard finite-precision eigenvalue routine may fail to compute the true Slepian tapers, but will still compute an orthonormal basis for their linear span. (See also [37] for examples where standard eigenvalue routines succeed or fail, and a discussion of possible workarounds.) As we shall see below, however, these tapers still yield the desired multitaper estimator, since the latter only depends on the tapers through their span.

Proposition 5.1**.**

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

[TABLE]

The proof is found in Appendix A.

While this explains the partial success of standard eigenvalue routines when applied to the ill-posed Slepian eigenproblem, the above result also indicates a more straightforward approach. Instead of trying to solve the standard diagonalization of $T^{\Omega,W}$—which may fail but still give usable tapers—we perform a block diagonalization. In other words, we calculate a basis for the span of the top $K$ eigenvectors of $T^{\Omega,W}$. Indeed, there is a large spectral gap between the first $K$ eigenvalues the rest of the spectrum (which Theorem 4.2 and the estimates in Appendix A validate to some extent). As a result, computing the associated subspace is a well-posed problem [52, 17].

We propose to compute these proxy Slepian tapers by applying a block power method to $T^{\Omega,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 Algorithm 1.

In exact arithmetic, the column space of $G$ converges to the span of the Slepian tapers $\{m_{0},\ldots,m_{K-1}\}$ as we increase $T$. Consequently, by Proposition 5.1, $\widehat{S}^{\mathrm{pmt}}(\xi)$ converges to $\widehat{S}^{\mathrm{mt}}(\xi)$. Since the gap between the $K$th and the $(K+1)$th eigenvalues of $T^{\Omega,W}$ is typically non-negligible [58, 59], subspace convergence is relatively fast, when errors are measured by the maximal singular value of the difference between the orthogonal projections onto the true and calculated subspaces (operator norm) [17, Theorem 8.2.2]. Although convergence of the column space of $G$ in operator norm is sufficient for accurate reproduction of the Slepian multitaper estimate, it is not necessary. Indeed, the estimates in Section A indicate that the multitaper estimator is robust under small errors in the calculation of the span of the Slepian tapers, when measured by the average of the singular values of the difference between the orthogonal projections onto the true and calculated subspaces (normalized trace norm). A moderate number of iterations $T$ may be sufficient to control such average and thus accurately calculate the multitaper estimator. Indeed, a single iteration $T=1$ is often sufficient for many applications. For special domains, such as rectangles, precise non-asymptotic estimates on the spectrum of $T^{\Omega,W}$ allow one to quantify the previous remarks, as in [58, 59, 28]. We expect that a similar analysis of the complexity of Algorithm 1 for general domains is also possible.

We also note that the only step in Algorithm 1 that depends on the data is the last one. As a result, we may precompute the proxy tapers $\widetilde{m}_{0},\ldots,\widetilde{m}_{K-1}$ once and apply them to the data as many times as necessary.

Remark 6**.**

*The matrix $T^{\Omega,W}$ is supported on ${\Omega}\times{\Omega}$; in Algorithm 1 it is treated as an element of $\mathbb{R}^{{\Omega}\times{\Omega}}$. *

Remark 7**.**

*Applying $T^{\Omega,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^{\Omega,W}$ is thus achieved in $\mathcal{O}({n_{\Omega}}\log{n_{\Omega}})$ time. *

Remark 8**.**

*Several authors have proposed numerical recipes to produce tapers adapted to irregular acquisition domains, resorting for example to QR decompositions and Karhunen–Loève expansions; see, e.g., [26, 18, 51]. We are unaware of corresponding performance results. *

6 Numerical results

To empirically evaluate the performance of the proposed proxy multitaper estimator, we perform a few numerical experiments on synthetic data. First, we numerically validate the theoretical results of Theorems 4.2 and 4.3. We then compare the proxy tapers to standard Slepian tapers on rectangular domains, where they are shown to perform similarly. Since the proxy tapers are easily computed on arbitrary domains, we also demonstrate its behavior on the complement of a disk. We also evaluate their mean squared error on single-particle cryo-EM data, for both synthetic and experimental images.

6.1 Empirical error analysis

To evaluate Theorem 4.2, we consider random fields on a grid ${\mathcal{Q}_{N}}$ with $N=128$ and $d=2$. The mask is a disk of radius $R$ given by ${\Omega}=\{q\in{\mathcal{Q}_{N}}:\left|q-(N/2,N/2)\right|<R\}$ and the bandwidth is fixed at $W=1/4$. For each radius, we then compute a set of proxy tapers $\widetilde{m}_{0},\ldots,\widetilde{m}_{K-1}$ and calculate their spectral window $\rho$, which we compare to the “ideal” window $W^{-d}1_{[-W/2,W/2]^{d}}$. The result is shown in Figure 4(a). Computing the average slope in the logarithmic plot shows that the error decays approximately as $R^{-1.0}$. This is close to the decay predicted by Theorem 4.2, whose leading term is $n_{\partial\Omega}W^{d-1}K^{-1}\approx n_{\partial\Omega}{n_{\Omega}}^{-1}\,W^{-1}$ and is proportional to $R^{-1}\,W^{-1}$.

We use the same set of masks to evaluate Theorem 4.3. The tapers are generated as before, but with $W={n_{\Omega}}^{-1/6}$, as recommended in Theorem 4.3. We simulate $M=200$ images $\mathcal{X}_{1},\ldots,\mathcal{X}_{M}$ by sampling a random field with a $1$-periodic $C^{2}$ spectral density $S$ given by

[TABLE]

where $J_{\ell}$ is the $\ell$th Bessel function. Each image $\mathcal{X}_{\nu}$ then yields a proxy multitaper spectral density estimate $\widehat{S}^{\mathrm{pmt}}_{\nu}$, all of which are used to estimate the mean squared error as

[TABLE]

for all $\xi$ on an $N\times N$ grid. These are then summarized by taking the maximum over $\xi$. The result is shown, as a function of $R$, in Figure 4(b). Dividing the error by $\log^{2}({n_{\Omega}})$ and computing the average slope in the logarithmic plot, we obtain a decay of about $R^{-1.5}$. Up to the logarithm factor and the norm of $S$, Theorem 4.3 gives an error bound of ${n_{\Omega}}^{-2/3}$, which is proportional to $R^{-4/3}$ and therefore close to the empirical decay.

6.2 Comparison with tensor Slepian tapers

As discussed in Section 3, tensor products of Slepian functions may be used when ${\Omega}$ is a subgrid of $\mathbb{Z}^{d}$. To illustrate the performance of the proxy tapers, we therefore first compare them to standard tensor Slepian tapers defined on the subgrid shown in Figure 5(a).

The top row Figure 6 shows a few tensor Slepian tapers defined on the subgrid domain of Figure 5(a) with target frequency profile given by 5(b), which corresponds to $W=1/8$. Below are proxy tapers defined on the same subgrid with the same frequency profile and calculated using Algorithm 1 for $T=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 7(a) and 7(b), both agree well with the target profile in Figure 5(b).

To evaluate the performance of these tapers in a multitaper estimation setting, we generate a Gaussian process with spectral density given by Figure 8(a) on a $128$-by-$128$-pixel square. The estimated density of the tensor Slepian multitaper estimator is given in Figure 8(b) while that of the proxy tapers for $T=2$ iterations is in Figure 8(c). Both agree quite well with the true density. Indeed, the normalized root mean squared errors $\|\widehat{S}^{\mathrm{mt}}-S\|/\|S\|$ and $\|\widehat{S}^{\mathrm{pmt}}-S\|/\|S\|$ for the tensor Slepian tapers and proxy tapers are both approximately $1.67\cdot 10^{-1}$. The deviation between the two estimators is $\|\widehat{S}^{\mathrm{mt}}-\widehat{S}^{\mathrm{pmt}}\|/\|\widehat{S}^{\mathrm{mt}}\|\approx 3.05\cdot 10^{-3}$. If we increase the number of iterations to $T=72$, we obtain equality up to machine precision with $\|\widehat{S}^{\mathrm{mt}}-\widehat{S}^{\mathrm{pmt}}\|/\|\widehat{S}^{\mathrm{mt}}\|\approx 5.42\cdot 10^{-16}$.

6.3 Illustration for irregular domain

We now replace the subgrid domain with the disk complement domain shown in Figure 5(c). Following the discussion of Section 3, it is computationally challenging to solve the eigenvalue (14) for the desired tapers (although standard eigenvalue solvers may still provide useful tapers, as observed in Section 5). We can, however, compute proxy tapers over this domain using Algorithm 1 with $T=2$. A few sample tapers are shown in the bottom row Figure 6. Again, their appearance is quite different from the Slepian tapers in the top row, but their accumulated spectral window shown in Figure 7(c) agrees well with the target of Figure 5(b).

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

6.4 Cryo-EM: Synthetic data

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 [16]. 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 a good characterization of the noise model. This is especially important for methods which estimate the covariance structure of the underlying images [9, 5]. Since the noise characteristics vary with microscope configuration, ice thickness, and other experimental factors, we cannot rely on ensemble averages for low-variance 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 9(a) and 9(b). We then generate Gaussian noise with a spectral density $\mathcal{X}$ given in Figure 9(c). Although the noise in cryo-EM projections is dominated by Poisson-distributed shot noise [6, 54], the number of registered electrons is large enough that a Gaussian approximation suffices for our purposes. The power spectral density was chosen to approximate those found in experimental datasets [4], with high energy in the low frequencies that decays quickly to zero at high frequencies. Adding the noise to the simulated projections, we obtain the images shown in Figure 9(d) and 9(e).

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 ${\Omega}$ of samples like the one shown in Figure 10(a). A common approach is to define the mask taper

[TABLE]

and use (4) to calculate the tapered periodogram $\widehat{S}_{m}(\xi)$ [57, 9]. 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 better estimator is obtained by replacing ${\Omega}$ with a union of rectangular subgrids ${\Omega^{\mathrm{grid}}}\subset{\Omega}$

[TABLE]

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 ${\Omega^{\mathrm{grid}}}$ corresponding to the ${\Omega}$ of Figure 10(a) is shown in Figure 10(b). Depending on the geometry of ${\Omega}$, ${\Omega^{\mathrm{grid}}}$ may discard many points in ${\Omega}$, 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 $\rho$.

Performance

We now compare the performance of these baseline estimators, the tapered periodogram $\widehat{S}_{m}$ and the tensor multitaper estimator $\widehat{S}^{\mathrm{mt}}$, to that of our proposed estimator $\widehat{S}^{\mathrm{pmt}}$. For a given $R$, we define ${\Omega}$ by (2) and ${\Omega^{\mathrm{grid}}}$ by (26). 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. We estimate the bias, variance, and MSE of the three estimators by computing them on $M=1000$ synthetic images of size $N=128$ generated as described in the beginning of the section. For the tensor multitaper and proxy taper estimators, we set $W=1/8$. We plot the resulting bias, variance, and MSE estimates as a function of $R$ in Figures 10(c), 10(d), and 10(e).

The bias of the proxy multitaper estimator is higher than the tensor multitaper for low $R$. This is due to the original mask ${\Omega}$ overlapping with the support of the projection images to a greater extent than ${\Omega^{\mathrm{grid}}}$ 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 proxy multitaper estimators and the tapered periodogram since ${\Omega}$ has greater overlap with the projections compared to ${\Omega^{\mathrm{grid}}}$. As $R$ increases, however, the proxy multitaper estimator $\widehat{S}^{\mathrm{pmt}}$ enjoys a lower variance, since it draws upon a larger number of samples compared to the tensor multitaper estimator $\widehat{S}^{\mathrm{mt}}$, 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 proxy multitaper estimator $\widehat{S}^{\mathrm{pmt}}$ compared to the other estimators. On average, the proxy multitaper estimator gives an error a factor of $1.7$ lower than the tensor multitaper.

6.5 Cryo-EM: Experimental data

We now evaluate performance on images consisting of experimental projections of an 80S ribosome complex from the EMPIAR-10028 dataset [55]. The images are defined on a grid with $N=360$. For our evaluation, we choose a subset of $M=120$ images from the dataset and set $W=1/16$. Two sample images from this subset are shown in Figure 2.

To obtain a reasonable approximation of the true power spectrum for each projection image, we process the entire dataset through the RELION software package [41]. 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 ${\mathcal{Q}_{N}}$. Since these estimates incorporate the noise in the center of the image, they should provide a good 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 $\widehat{S}_{m}$ on ${\Omega}$, the tensor multitaper estimator $\widehat{S}^{\mathrm{mt}}$ on ${\Omega^{\mathrm{grid}}}$, and the proxy multitaper estimator $\widehat{S}^{\mathrm{pmt}}$ on ${\Omega}$. Each estimated power spectrum is compared to the “ground truth” power spectrum obtained for the corresponding image as described above, which gives an estimated MSE for each method. Similarly to the simulation results, the error for the tapered periodogram $\widehat{S}_{m}$ is dominated by the variance, resulting in an MSE of around $10$ largely independent of $R$. The MSE’s for the tensor multitaper and proxy multitaper estimators are plotted for each estimator in Figure 11.

The tensor multitaper estimator performs well for low $R$, but increasing $R$ reduces the area of ${\Omega^{\mathrm{grid}}}$, yielding higher variance and higher MSE. A similar behavior is observed for the proxy multitaper estimates, but the reduction in area for ${\Omega}$ is not as drastic for increasing $R$, so the variance remains small compared to the tensor multitaper estimator. At $R\approx 96$, the error is minimized and the proxy multitaper outperforms the tensor multitaper by a factor of $1.7$.

7 Conclusion

We have analyzed the multitaper estimator on arbitrary acquisition domains, providing performance bounds on the mean squared error. Furthermore, we show that the multitaper estimate only depends on the tapers through their linear span. This explains the success of applying standard eigenvalue algorithms to the ill-posed Slepian eigenproblem, since a common mode of failure is for these to yield a set of vectors with the same span as the desired eigenvectors. Using the resulting vectors as tapers therefore yields results close to those obtained with the true Slepian tapers. We also present a more straightforward approach of calculating these proxy Slepian tapers using the block power method. The performance of the resulting proxy multitaper estimator is shown to be comparable to that using Slepian tapers when these are available for rectangular domains. We also illustrate the performance of the proxy multitaper estimator for more general domains and compared it to tensor Slepian multitaper estimators on rectangular subgrids. Numerical results are obtained 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 random field. 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.

Acknowledgments

The authors are very grateful to Luís Daniel Abreu, Tomasz Hrycak, Frederik Simons, and Amit Singer, who motivated this article and provided valuable input. They would also like to thank the anonymous reviewers for their helpful remarks. J. L. R. gratefully acknowledges support from the Austrian Science Fund (FWF): P 29462 and Y 1199, and from the WWTF grant INSIGHT (MA16-053). The Flatiron Institute is a division of the Simons Foundation.

Appendix A Proofs

The $p$-norm of a vector $x\in\mathbb{R}^{d}$ is denoted $\left|x\right|_{p}$. For a function $f\in\ell^{1}(\mathbb{Z}^{d})$, we denote its $\ell^{1}$-norm $\sum_{q}|f[q]|$ by $|f|$.

For two non-negative functions $f,g:X\to[0,+\infty)$, we write $f\lesssim g$ if there exists a constant $C\geq 0$ such that $f\leq Cg$. We also write $f\asymp g$, if $f\lesssim g$ and $g\lesssim f$.

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

[TABLE]

The $L^{\infty}$-norm and $L^{1}$-norm of a measurable, 1-periodic function $f:\mathbb{R}^{d}\to\mathbb{C}$ are, respectively,

[TABLE]

For clarity, we sometimes write $\lVert f\rVert_{L^{\infty}([-1/2,1/2]^{d})}$ and $\lVert f\rVert_{L^{1}([-1/2,1/2]^{d})}$.

A.1 Convolution estimates

First, we note that

[TABLE]

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$ as

[TABLE]

An estimate on a second-order Taylor expansion shows that

[TABLE]

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 [33] by avoiding assumptions on the (digital) topology of the set ${\Omega}$. See also [34] for related estimates.

Proposition A.1**.**

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

[TABLE]

*where $\left(\nabla_{k}b\right)[q]:=b[q+e_{k}]-b[q]$ for $k=1,\ldots,d$, $\lVert\nabla b\rVert_{1}=\sum_{k=1}^{d}\lVert\nabla_{k}b\rVert_{1}$, and $\left|q\right|_{\infty}=\max\{\left|q_{1}\right|,\ldots,\left|q_{d}\right|\}$. *

Proof A.2**.**

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

[TABLE]

Therefore,

[TABLE]

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

Proposition A.3**.**

Let $T^{\Omega,W}$ be the matrix in (12). Then

[TABLE]

Proof A.4**.**

Step 1. (Computations)**. The function

[TABLE]

determines the Fourier series

[TABLE]

Note that, in terms of $\mathrm{h}$, (12) reads: $T^{\Omega,W}_{q,q^{\prime}}=1_{\Omega}(q)\mathrm{h}[q-q^{\prime}]1_{\Omega}(q^{\prime})$. We first compute

[TABLE]

Second,

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

We first note the following:

[TABLE]

To show (36), we exploit the fact that the $d$-ball of radius $L$ contains a $d$-cube of side $2\cdot d^{-1/2}L$ and estimate

[TABLE]

Similarly, for $L\geq 2$

[TABLE]

which gives (37) and establishes the remaining claim.

Step 3. (Final estimates)**. Let $b:=\lVert\tilde{a}\rVert_{1}^{-1}\tilde{a}$. Noting that $\sum_{q}b[q]=1$, we may combine (33) with (35) and use Proposition A.1 to form the bound

[TABLE]

where we have also used the fact that $\lVert 1_{\Omega}*(a-b)\rVert_{1}\leq\lVert 1_{\Omega}\rVert_{1}\lVert a-b\rVert_{1}$.

Let $\varepsilon:=\lVert a-\tilde{a}\rVert_{1}$. Since $\lVert a\rVert_{1}=1$ and $\tilde{a}[q]\leq a[q]$ for all $a\in{\mathbb{Z}}$, we have $\varepsilon=1-\lVert\tilde{a}\rVert_{1}$. Then, by (36), $\varepsilon\lesssim(WL)^{-1}$. Hence, there exist a constant $C_{0}>2$ such that $\varepsilon<1/2$, if $WL\geq C_{0}$. Let us assume that for the moment that $WL\geq C_{0}$, so that $\varepsilon<1/2$ and $\lVert\tilde{a}\rVert_{1}=1-\varepsilon>1/2$, and estimate

[TABLE]

Similarly, by (37),

[TABLE]

Substituting these estimates into (38) gives

[TABLE]

The right-hand side is minimized at $L=\tfrac{{n_{\Omega}}}{n_{\partial\Omega}}$, which yields (30), provided that $WL=W\tfrac{{n_{\Omega}}}{n_{\partial\Omega}}\geq C_{0}$. On the other hand, if $W\tfrac{{n_{\Omega}}}{n_{\partial\Omega}}\leq C_{0}$, (30) is trivially true because

[TABLE]

A.3 Expectation and variance of the tapered estimators

Proof A.5** (Proof of Proposition 4.1).**

The fact that $\rho$ depends only on the linear span the tapers follows from the proof of Proposition 5.1 below. The variance bound can be proved as in [32, Theorem 2]. Only the orthogonality of the tapers is important here (see also [24, Chapter 3]). For the bias of a taper $m\in\ell^{1}({{\mathbb{Z}}^{d}})$, a direct calculation yields:

[TABLE]

*By averaging this expression over all tapers, we then obtain the bias of the multitaper estimator. *

A.4 Description of the spectral window

Proof A.6** (Proof of Theorem 4.2).**

We use the notation $M_{k}(\xi)=\sum_{q\in{{\mathbb{Z}}^{d}}}m_{k}[q]\mathrm{e}^{2\pi\mathrm{i}\left<q,\xi\right>}$.

Step 1*. We use that $\{m_{0},\ldots,m_{{n_{\Omega}}-1}\}$ is an orthonormal basis of $\ell^{2}(\Omega)$ and compute*

[TABLE]

since an orthonormal change of basis preserves the norm of the function. We then note that

[TABLE]

This lets us form the estimate

[TABLE]

Similarly,

[TABLE]

Hence,

[TABLE]

Step 2*. Using (39), we estimate,*

[TABLE]

Step 3*. We now proceed as in [3]:*

[TABLE]

By Proposition A.3,

[TABLE]

*where constant vanishes since $n_{\partial\Omega}W^{d-1}\geq 1$. Finally, we combine this with (40) to obtain (17). *

A.5 Proof of Theorem 4.3

Proof A.7**.**

Using that $W\asymp\left(\tfrac{K}{{n_{\Omega}}}\right)^{1/d}$, we invoke Theorem 4.2, (29) and (28) to obtain

[TABLE]

Combining this with Proposition 4.1, we obtain

[TABLE]

Finally for $d=2$, if $K\approx{n_{\Omega}}^{2/3}$, and ${n_{\Omega}}\leq n_{\partial\Omega}^{2}\leq C{n_{\Omega}}$, for some constant $C>0$, we obtain

[TABLE]

A.6 Proof of Proposition 5.1

Let $U\in\mathbb{R}^{L\times L}$ be an orthogonal matrix such that

[TABLE]

Using that $U^{\mathrm{T}}U=I$, we compute

[TABLE]

as claimed.