Factor Analysis for Spectral Estimation

TL;DR

This paper introduces a subspace projection method for spectral estimation that improves accuracy for short signals by leveraging a statistical model of fixed sources, with demonstrated success in cryo-electron microscopy data.

Contribution

It proposes a novel subspace-based spectral estimation technique with theoretical guarantees, addressing limitations of existing methods for short signals.

Findings

Enhanced spectral estimation accuracy for short signals.

Theoretical guarantees for the proposed method.

Successful application to cryo-electron microscopy data.

Abstract

Power spectrum estimation is an important tool in many applications, such as the whitening of noise. The popular multitaper method enjoys significant success, but fails for short signals with few samples. We propose a statistical model where a signal is given by a random linear combination of fixed, yet unknown, stochastic sources. Given multiple such signals, we estimate the subspace spanned by the power spectra of these fixed sources. Projecting individual power spectrum estimates onto this subspace increases estimation accuracy. We provide accuracy guarantees for this method and demonstrate it on simulated and experimental data from cryo-electron microscopy.

Equations

$$\mathbf{X}[\vec{ı}]:=\sum _{l=1}^r{a}_l{\mathbf{Z}}_l[\vec{ı}]\vec{ı}\in {\mathbb{Z}}^d\text{mbox},$$

$$R_{\mathbf{Y}}[\vec{ȷ}]:=\mathbb{E}\left[\mathbf{Y}[\vec{ı}]\mathbf{Y}[\vec{ı}+\vec{ȷ}]\right]\vec{ȷ}\in {\mathbb{Z}}^d.$$

$$P_{\mathbf{Y}}(\vec{\xi }):=\sum _{\vec{ı}\in {\mathbb{Z}}^d}{R}_{\mathbf{Y}}[\vec{ı}]{\mathrm{e}}^{-2\pi \mathrm{i}⟨\vec{\xi },\vec{ı}⟩}\vec{\xi }\in [-1/2,1/2{]}^d.$$

$${\hat{P}}_{{\mathbf{Y}}^{(N)}}[\vec{k}]=\frac{1}{N^d}{\left|\sum _{\vec{ı}\in M_N^d}{\mathbf{Y}}^{(N)}[\vec{ı}]{\mathrm{e}}^{-2\pi \mathrm{i}⟨\vec{k},\vec{ı}⟩/N}\right|}^2\vec{k}\in M_N^d\text{mbox}.$$

$$\mathbf{Y}[\vec{ı}]=\sum _{\vec{ȷ}\in {\mathbb{Z}}^d}W[\vec{ı}]\psi [\vec{ı}-\vec{ȷ}]\text{mbox},$$

$$\mathbb{E}{\hat{P}}_{{\mathbf{Y}}^{(N)}}[\vec{k}]=P_{\mathbf{Y}}(\vec{k}/N)+O(N^{-1/2})\text{mbox}.$$

$$\delta_{N}[\vec{k}_{1},\vec{k}_{2}]=\left\{\begin{array}[]{ll}2&\vec{k}_{1}=\vec{k}_{2}\mbox{~{}and~{}}\vec{k}_{1}\in\frac{N}{2}\mathbb{Z}^{d}\\ 1&\vec{k}_{1}=\vec{k}_{2}\mbox{~{}and~{}}\vec{k}_{1}\notin\frac{N}{2}\mathbb{Z}^{d}\\ 0&\vec{k}_{1}\neq\vec{k}_{2}\end{array}\right.\mbox{.}$$

$$\mathrm{Cov}\left[{\hat{P}}_{{\mathbf{Y}}^{(N)}}[{\vec{k}}_1],{\hat{P}}_{{\mathbf{Y}}^{(N)}}[{\vec{k}}_2]\right]$$

$$={\delta }_N[{\vec{k}}_1,{\vec{k}}_2]P_{\mathbf{Y}}({\vec{k}}_1/N{)}^2+O(N^{-1/2})\text{mbox}.$$

$$v_{N,W}^{(\vec{l})}[\vec{ı}]=\prod _{p=1}^d{v}_{N,W}^{(l_p)}[i_p]\vec{ı}=(i_1,\dots ,i_d)\in M_N^d\text{mbox},$$

$${\hat{P}}_{{\mathbf{Y}}^{(N)}}^{(\mathrm{MT})}[\vec{k}]=\frac{1}{K^d}\sum _{\vec{l}\in \{1,\dots ,K{\}}^d}{\hat{P}}_{v_{N,W}^{(\vec{l})}\cdot {\mathbf{Y}}^{(N)}}[\vec{k}]\text{mbox},$$

$${\mathbf{X}}_s[\vec{ı}]:=\sum _{l=1}^r{a}_{s,l}{\mathbf{Z}}_{s,l}[\vec{ı}]\text{mbox},$$

$$P_{{\mathbf{X}}_s|{\vec{a}}_s}(\vec{\xi })=\sum _{l=1}^r{a}_{s,l}^2{P}_{{\mathbf{Z}}_l}(\vec{\xi })\text{mbox}.$$

$${\mu }_n[\vec{k}]=\frac{1}{n}\sum _{s=1}^n{\hat{P}}_{{\mathbf{X}}_s^{(N)}}[\vec{k}]$$

$$C_n[{\vec{k}}_1,{\vec{k}}_2]=\frac{1}{n}\sum _{s=1}^n{\hat{P}}_{{\mathbf{X}}_s^{(N)}}[{\vec{k}}_1]{\hat{P}}_{{\mathbf{X}}_s^{(N)}}[{\vec{k}}_2]\text{mbox},$$

$${\mu }_n[\vec{k}]=P_{\mathbf{X}}(\vec{k}/N)+{ϵ}_1[\vec{k}]$$

$$\frac{1}{1+\delta [{\vec{k}}_1,{\vec{k}}_2]}{C}_n[{\vec{k}}_1,{\vec{k}}_2]-{\mu }_n[{\vec{k}}_1]{\mu }_n[{\vec{k}}_2]$$

$$=\mathrm{Cov}\left[P_{\mathbf{X}|\vec{a}}({\vec{k}}_1/N),P_{\mathbf{X}|\vec{a}}({\vec{k}}_2/N)\right]+{ϵ}_2[{\vec{k}}_1,{\vec{k}}_2]\text{mbox},$$

$${\mathbb{E}}_{{\mathbf{Z}}_1,\dots ,{\mathbf{Z}}_r|\vec{a}}{\hat{P}}_{{\mathbf{X}}^{(N)}}[\vec{k}]=P_{\mathbf{X}|\vec{a}}(\vec{k}/N)+O(N^{-1/2})\text{mbox},$$

$$\mathbb{E}{\hat{P}}_{{\mathbf{X}}^{(N)}}[\vec{k}]=P_{\mathbf{X}}(\vec{k}/N)+O(N^{-1/2})\text{mbox}.$$

$${\mu }_n[\vec{k}]=\frac{1}{n}\sum _{s=1}^n{\hat{P}}_{{\mathbf{X}}_s^{(N)}}[\vec{k}]=\mathbb{E}{\hat{P}}_{{\mathbf{X}}^{(N)}}[\vec{k}]+{ϵ}_1[\vec{k}]\text{mbox},$$

$${\mathrm{Cov}}_{{\mathbf{Z}}_1,\dots ,{\mathbf{Z}}_r|\vec{a}}\left[{\hat{P}}_{{\mathbf{X}}^{(N)}}[{\vec{k}}_1],{\hat{P}}_{{\mathbf{X}}^{(N)}}[{\vec{k}}_2]\right]$$

$$=\delta [{\vec{k}}_1,{\vec{k}}_2]P_{\mathbf{X}|\vec{a}}({\vec{k}}_1/N{)}^2+O(N^{-1/2})\text{mbox}.$$

$${\mathbb{E}}_{{\mathbf{Z}}_1,\dots ,{\mathbf{Z}}_r|\vec{a}}\left[{\hat{P}}_{{\mathbf{X}}_s^{(N)}}[{\vec{k}}_1]{\hat{P}}_{{\mathbf{X}}_s^{(N)}}[{\vec{k}}_2]\right]$$

$$=(1+\delta [{\vec{k}}_1,{\vec{k}}_2])P_{\mathbf{X}|\vec{a}}({\vec{k}}_1/N)P_{\mathbf{X}|\vec{a}}({\vec{k}}_2/N)+O(N^{-1/2})\text{mbox}.$$

$$\mathbb{E}\left[{\hat{P}}_{{\mathbf{X}}_s^{(N)}}[{\vec{k}}_1]{\hat{P}}_{{\mathbf{X}}_s^{(N)}}[{\vec{k}}_2]\right]$$

$$=(1+\delta [{\vec{k}}_1,{\vec{k}}_2])\mathbb{E}\left[P_{\mathbf{X}|\vec{a}}({\vec{k}}_1/N)P_{\mathbf{X}|\vec{a}}({\vec{k}}_2/N)\right]+O(N^{-1/2})\text{mbox}.$$

$$\frac{1}{1+\delta [{\vec{k}}_1,{\vec{k}}_2]}\mathbb{E}\left[{\hat{P}}_{{\mathbf{X}}^{(N)}}[{\vec{k}}_1]{\hat{P}}_{{\mathbf{X}}^{(N)}}[{\vec{k}}_2]\right]-P_{\mathbf{X}}({\vec{k}}_1/N)P_{\mathbf{X}}({\vec{k}}_2/N)$$

$$=\mathrm{Cov}\left[P_{\mathbf{X}|\vec{a}}({\vec{k}}_1/N),P_{\mathbf{X}|\vec{a}}({\vec{k}}_2/N)\right]+O(N^{-1/2})\text{mbox}.$$

$${\Sigma }_n[{\vec{k}}_1,{\vec{k}}_2]=\frac{1}{1+\delta [{\vec{k}}_1,{\vec{k}}_2]}{C}_n[{\vec{k}}_1,{\vec{k}}_2]-{\mu }_n[{\vec{k}}_1]{\mu }_n[{\vec{k}}_2]\text{mbox},$$

$${\hat{P}}_{{\mathbf{X}}^{(N)}}^{(\mathrm{MT})}=P_{\mathbf{X}|\vec{a}}+R^{(\mathrm{MT})}\text{mbox},$$

$${\hat{P}}_{{\mathbf{X}}^{(N)}}^{(\mathrm{MT},\mathrm{lin})}[\vec{k}]=\sum _{l=1}^r{v}_l[\vec{k}]⟨v_l,{\hat{P}}_{{\mathbf{X}}^{(N)}}^{(\mathrm{MT})}⟩\text{mbox}.$$

Full text

I Introduction

Estimating the power spectrum of a stationary field arises in many applications [1, 2, 3, 4]. For example, noise statistics play an important role in maximum likelihood estimation of signal parameters. Hence, its power spectrum must be estimated.

The most basic estimator of the power spectrum is the periodogram, which is asymptotically unbiased but has very high variance. To remedy this, spectral smoothness constraints or parametric models are often imposed [5, 6]. The multitaper estimator [7, 6, 8] is particularly popular, which multiplies the data with fixed tapers, computes their periodograms, and averages. This reduces variance at the expense of smoothing the estimated power spectrum, potentially increasing bias.

In this work, we consider the problem of estimating the power spectra of multiple independent signals. Given identically distributed signals, an ensemble average of their periodograms or multitaper estimates provides a low-variance power spectrum estimate. However, in certain applications, signals are not identically distributed, but combine a small number of identically distributed stochastic sources. For example, a non-stationary signal may be divided up into approximately stationary parts, each of which is a linear combination of some fixed random signals. Alternatively, noise signals may arise from measurements subject to differing experimental conditions which are described by combining various noise sources. This is the case for cryo-electron microscopy (cryo-EM) images, where large molecules are frozen in a thin layer of ice and then imaged by exposing them to an electron beam and recording the transmitted electrons [9, 10]. Variations in experimental conditions result in different noise characteristics. To whiten the projection images, accurate estimates of their noise power spectrum are therefore necessary.

While traditional methods such as multitaper estimates are applied in this context, an estimator which takes into account the low-dimensional variability of the noise distributions could yield improved accuracy. We propose a factor analysis model in which the $d$-dimensional field $\mathbf{X}:\mathbb{Z}^{d}\rightarrow\mathbb{R}$ is a linear combination of $r$ fixed fields

[TABLE]

where $a_{1},\ldots,a_{r}$ are independent random variables and $\mathbf{Z}_{1},\ldots,\mathbf{Z}_{r}$ are independent linear random fields. Typical values for $d$ include $d=1$ for stationary processes and $d=2$ for random images.

It follows that the power spectrum of $\mathbf{X}$ conditioned on $a_{1},\ldots,a_{r}$ is a linear combination of the power spectra of $\mathbf{Z}_{1},\ldots,\mathbf{Z}_{r}$. That is, the conditional power spectra reside in a low-dimensional subspace and are therefore determined by a small number of factors. We introduce a method for estimating this subspace and these factors given multiple realizations of $\mathbf{X}$. The method allows for estimation of the number of factors, which is small for most applications. Linear projection of individual power spectrum estimates onto the factor subspace then yields improved accuracy.

Section II discusses the power spectrum estimation problem for a single signal and describes the basic periodogram and multitaper methods. The factor analysis model and associated power spectrum estimation method are introduced and analyzed in Section III. Finally, Section IV presents numerical results on simulations and experimental data taken from cryo-EM images. Software implementing our proposed method and reproducing the figures of this paper is available at http://github.com/janden/fase/.

II Single Power Spectrum Estimation

Given a zero-mean stationary random field $\mathbf{Y}:\mathbb{Z}^{d}\rightarrow\mathbb{R}$ with finite second moments, we define its autocovariance as

[TABLE]

Its Fourier transform is the power spectrum of $\mathbf{Y}$:

[TABLE]

The power spectrum $P_{\mathbf{Y}}$ is symmetric and non-negative.

While $\mathbf{Y}$ is defined over $\mathbb{Z}^{d}$, we are only given samples over some finite domain $M_{N}^{d}$, where $M_{N}=\{-\lceil N/2\rceil+1,\ldots,\lfloor N/2\rfloor\}$. Let us denote this restriction of $\mathbf{Y}$ to $M_{N}^{d}$ by $\mathbf{Y}^{(N)}$. Our task is then to estimate $P_{\mathbf{Y}}$ given $\mathbf{Y}^{(N)}$.

II-A Periodogram

The most basic estimator for the power spectrum of a stationary field is the periodogram, given by

[TABLE]

Since $\mathbf{Y}^{(N)}$ is real, $\widehat{P}_{\mathbf{Y}^{(N)}}$ is symmetric and only needs to be computed on half of $M_{N}^{d}$. Let us therefore choose a subdomain $M_{N,+}^{d}$ of $M_{N}^{d}$ such that $M_{N,+}^{d}\cup-M_{N,+}^{d}=M_{N}~{}(\mathrm{mod}~{}N)$ and $M_{N,+}^{d}\cap-M_{N,+}^{d}=\{\vec{0}\}$.

To study the properties of the periodogram, we first define a linear stationary field.

Definition 1**.**

A stationary field $\mathbf{Y}:\mathbb{Z}^{d}\rightarrow\mathbb{R}$ is linear if

[TABLE]

where $W:\mathbb{Z}^{d}\rightarrow\mathbb{R}$ is a stationary, zero-mean, white noise field such that $\mathbb{E}|W[\vec{\imath}]|^{2}=1,\mathbb{E}|W[\vec{\imath}]|^{4}<\infty$ and $\psi:\mathbb{Z}^{d}\rightarrow\mathbb{R}$ is a kernel satisfying $\sum_{\vec{\imath}\in\mathbb{Z}^{d}}|\psi[\vec{\imath}]|\|\vec{\imath}\|_{1}^{1/2}<\infty$.

For such a field, the periodogram is an asymptotically unbiased estimator for the power spectrum:

Lemma 1**.**

Let $\mathbf{Y}:\mathbb{Z}^{d}\rightarrow\mathbb{R}$ be a linear field and let $\mathbf{Y}^{(N)}$ be its restriction $M_{N}^{d}$. Then we have, for $\vec{k}\in M_{N,+}^{d}$,

[TABLE]

The proof is obtained as a generalization of the case $d=1$ found in Proposition 10.3.1 of [5]. Recasting the periodogram as the Fourier transform of sums of pairs $\mathbf{Y}^{(N)}[\vec{\imath}]\mathbf{Y}^{(N)}[\vec{\imath}+\vec{\jmath}]$, the bias is due to not weighting these sums by the number of available pairs. Reweighting eliminates the bias, but at a significant increase in variance, increasing the mean squared error of the estimator. For this reason, the biased periodogram estimator is often preferred over its unbiased variant [6].

Despite its low bias for large $N$, the periodogram is not suitable for power spectrum estimation due its high variance:

Lemma 2**.**

Let $\mathbf{Y}$ and $\mathbf{Y}^{(N)}$ be defined as in Lemma 1 and

[TABLE]

We then have, for all $\vec{k}_{1},\vec{k}_{2}\in M_{N,+}^{d}$,

[TABLE]

The proof again follows the one-dimensional case, Proposition 10.3.2 in [5]. It relies on the fact that as $N\rightarrow\infty$, $\widehat{P}_{\mathbf{Y}^{(N)}}[\vec{k}]$ converges to a $\chi^{2}$-distributed random variable. Asymptotically, the standard deviation of $\widehat{P}_{\mathbf{Y}^{(N)}}[\vec{k}]$ is therefore approximately proportional to its expectation $P_{\mathbf{Y}}(\vec{k}/N)$.

II-B The Multitaper Method

The multitaper method for spectral estimation was introduced by D. J. Thomson as an alternative to the periodogram with reduced variance at the cost of increased bias [7]. It relies on multiplying the sample $\mathbf{Y}^{(N)}$ by a number of data tapers, computing their periodograms, and averaging.

Given a target spectral resolution $\frac{1}{2N}<W<\frac{1}{2}$, the data tapers for $d=1$ are given by $K=\lfloor 2NW\rfloor$ discrete prolate spheroidal sequences $v_{N,W}^{(l)}[i]$ for $i\in M_{N}$, where $l\in 1,\ldots,K$. Abusing notation slightly, we denote their tensor products by

[TABLE]

where $\vec{l}=(l_{1},\ldots,l_{d})\in\{1,\ldots,K\}^{d}$ [11]. The associated multitaper estimator is then given by

[TABLE]

where $(v_{N,W}^{(\vec{l})}\cdot\mathbf{Y}^{(N)})[\vec{\imath}]=v_{N,W}^{(\vec{l})}[\vec{\imath}]\mathbf{Y}^{(N)}[\vec{\imath}]$ for all $\vec{\imath}\in M_{N}^{d}$.

The resolution parameter $W$ specifies the bias-variance trade-off of the estimator. As $W$ is increased, variance is reduced, but the estimate is smoothed, potentially increasing the bias of the estimator.

III Power Spectrum Factor Analysis

In certain applications, we do not have just one signal whose power spectrum we would like to estimate, but a set of signals. To estimate their power spectra, we can apply the methods described in the previous section to each signal individually. If in addition the signals are identically distributed, we can reduce the variance in the estimate by an ensemble average.

On the other hand, if the signals are not identically distributed, an ensemble average will destroy the variability in the power spectra. A different approach is applicable if the signals can be written as linear combinations of a small set of fixed random fields. In this case, factor analysis allows us to identify the subspace spanned by the power spectra of these fixed fields, allowing for increased accuracy by projecting multitaper estimates onto this space.

III-A Motivation: Cryo-Electron Microscopy

Single-particle cryo-EM images are typically very noisy, with signal-to-noise ratios typically below $1/10$ [9, 10]. Many reconstruction algorithms assume that this noise is close to white, but this is rarely the case. As a result, a whitening step is typically performed prior to reconstruction, in which the noise power spectrum is estimated and used to compute a whitening filter that is then applied. Power spectrum estimation is typically done by extracting image patches containing mostly noise and averaging their periodograms [12]. Other methods estimate the noise power spectrum as part of the reconstruction algorithm [13].

However, the noise distribution typically varies between projection images. Indeed, different projection images are obtained under different experimental conditions (ice thickness, electron dose, microscope defocus, electron scattering properties of the molecules, and so on) which affect the characteristics of the noise [15, 16, 17]. The variation in the noise signal can be modeled as a random linear combination of fixed noise sources due to background electrons, electrons scattered by the ice, and those inelastically scattered by the molecule. Experimental images with different noise power spectra are shown in Figure 1. As we will see in the following, the power spectra of these noise signals can be accurately estimated using the proposed factor analysis.

III-B Problem Setup

Let us denote the zero-mean stationary random field of interest by $\mathbf{X}:\mathbb{Z}^{d}\rightarrow\mathbb{R}$. We model it in (1) as a linear combination of $r$ fixed independent linear fields $\mathbf{Z}_{1},\ldots,\mathbf{Z}_{r}$, with independent random coefficients $a_{1},\ldots,a_{r}$.

Instead of a single field $\mathbf{X}$, we are concerned with a set of $n$ independent copies $\mathbf{X}_{1},\ldots,\mathbf{X}_{n}$ of $\mathbf{X}$ defined by

[TABLE]

for $s=1,\ldots,n$. Here $\mathbf{Z}_{1,l},\ldots,\mathbf{Z}_{n,l}$ and $a_{1,l},\ldots,a_{n,l}$ are independent and identically distributed copies of $\mathbf{Z}_{l}$ and $a_{l}$, respectively, for $l=1,\ldots,r$. The given data consists of restrictions $\mathbf{X}^{(N)}_{1},\ldots,\mathbf{X}^{(N)}_{n}$ of $\mathbf{X}_{1},\ldots,\mathbf{X}_{n}$ to $M_{N}^{d}$.

The coefficient vectors $\vec{a}_{1},\ldots,\vec{a}_{n}$ are unknown, but fixed for each signal. As a result, the power spectrum of interest for $\mathbf{X}_{s}$ is the conditional power spectrum $P_{\mathbf{X}_{s}|\vec{a}_{s}}$, where we have replaced the expectation $\mathbb{E}$ in (2) by the conditional expectation $\mathbb{E}_{\mathbf{Z}_{s,1},\ldots,\mathbf{Z}_{s,r}|\vec{a}_{s}}$. Our problem is therefore estimating $P_{\mathbf{X}_{1}|\vec{a}_{1}},\ldots,P_{\mathbf{X}_{n}|\vec{a}_{n}}$ given $\mathbf{X}^{(N)}_{1},\ldots,\mathbf{X}^{(N)}_{n}$.

One approach is to use the multitaper method described in Section II-B applied to each signal $\mathbf{X}^{(N)}_{1},\ldots,\mathbf{X}^{(N)}_{n}$. However, we can obtain improved estimates by exploiting the fact that, for $s=1,\ldots,n$,

[TABLE]

That is, the conditional power spectra are contained in a subspace of dimension at most $r$ spanned by $P_{\mathbf{Z}_{1}},\ldots,P_{\mathbf{Z}_{r}}$. These power spectra can thus be described using a small set (at most $r$) of common factors. We can exploit this to increase the accuracy of a power spectrum estimate by projecting it onto this low-dimensional subspace. In the following, we propose a method for estimating this subspace.

III-C Power Spectrum Covariance Estimation

The non-trivial eigenvectors of the covariance matrix $\mathrm{Cov}\left[P_{\mathbf{X}_{s}|\vec{a}_{s}}\right]$ span the subspace containing the conditional power spectra. To simplify our expressions, we define $P_{\mathbf{X}|\vec{a}}$ in the same way as $P_{\mathbf{X}_{s}|\vec{a}_{s}}$ for $s=1,\ldots,n$ and note that these are all identically distributed. We consider estimating the covariance matrix $\mathrm{Cov}\left[P_{\mathbf{X}|\vec{a}}\right]$ using the following theorem:

Theorem 1**.**

Let $\mathbf{X}$ and $\mathbf{X}_{1},\ldots,\mathbf{X}_{n}$ be as in (1) and (11), respectively. Furthermore, let $\mathbf{X}^{(N)}_{1},\ldots,\mathbf{X}^{(N)}_{n}$ be the restrictions of $\mathbf{X}_{1},\ldots,\mathbf{X}_{n}$ to $M_{N}^{d}$. In this case, the quantities

[TABLE]

and

[TABLE]

satisfy

[TABLE]

and

[TABLE]

where $\vec{k},\vec{k}_{1},\vec{k}_{2}\in M_{N,+}^{d}$. The expected magnitudes $\mathbb{E}|\epsilon_{1}[\vec{k}]|$ and $\mathbb{E}|\epsilon_{2}[\vec{k}_{1},\vec{k}_{2}]|$ are each bounded by $O(n^{-1/2}+N^{-1/2})$.

Proof.

First, we condition Lemma 1 on $\vec{a}$ to obtain

[TABLE]

for $\vec{k}\in M_{N,+}^{d}$. Taking expectation with respect to $\vec{a}$ then gives

[TABLE]

Averaging $\widehat{P}_{\mathbf{X}^{(N)}_{1}}[\vec{k}],\ldots,\widehat{P}_{\mathbf{X}^{(N)}_{n}}[\vec{k}]$ then gives us,

[TABLE]

according to the central limit theorem, where $\mathbb{E}|\epsilon_{1}[\vec{k}]|=O(n^{-1/2})$. Plugging (18) into (19) then proves (15).

Conditioning Lemma 2 with respect to $\vec{a}$ gives

[TABLE]

Consequently

[TABLE]

Taking the expectation with respect to $\vec{a}$ now gives

[TABLE]

We now have

[TABLE]

Replacing the left-hand side expectations with sums using the central limit theorem yields (16). ∎

The theorem allows us to estimate the covariance matrix of the conditional power spectra. Traditionally, this would be done by forming $C_{n}[\vec{k}_{1},\vec{k}_{2}]-\mu_{n}[\vec{k}_{1}]\mu_{n}[\vec{k}_{2}]$. However, due to the $\chi^{2}$ distribution of the periodogram coordinates, the variances of $C_{n}$ are larger than they should be, so a correction is needed. That correction takes the form of $(1+\delta[\vec{k}_{1},\vec{k}_{2}])^{-1}$.

Let us denote this covariance matrix estimate by

[TABLE]

for $\vec{k}_{1},\vec{k}_{2}\in M_{N,+}^{d}$. Calculating its top eigenvectors allows us to estimate the low-dimensional subspace containing the conditional power spectrum $P_{\mathbf{X}|\vec{a}}$. Specifically, these eigenvectors define a linear space that approximates the subspace spanned by $P_{\mathbf{Z}_{1}},\ldots,P_{\mathbf{Z}_{r}}$.

III-D Linear Projection Estimators

We now use the information gained in the previous section to obtain a better estimate of the conditional power spectrum $P_{\mathbf{X}|\vec{a}}$ from the multitaper estimate $\widehat{P}_{\mathbf{X}^{(N)}}^{\mathrm{(MT)}}$. It satisfies

[TABLE]

where the size of the residual $R^{\mathrm{(MT)}}$ depends on $N$ and the regularity of $P_{\mathbf{X}|\vec{a}}$ with respect to the target multitaper resolution $W$ [6]. We know that $P_{\mathbf{X}|\vec{a}}$ is in the subspace spanned by the non-trivial eigenvectors of $\mathrm{Cov}\left[P_{\mathbf{X}|\vec{a}}\right]$. We can use this to reduce the variance of the residual term.

While the exact covariance is not known to us, we have an estimator $\Sigma_{n}$. Let $v_{1},\ldots,v_{r}$ be the leading eigenvectors of $\Sigma_{n}$. We note that $r$ does not need to be known beforehand, but can be estimated from the eigenspectrum of $\Sigma_{n}$. Here, we set $r$ to be the number of dominant eigenvalues of $\Sigma_{n}$ by locating the “knee” of the spectrum, where there is a significant drop in amplitude from the $r$th eigenvalue to the $(r+1)$th. In addition, since we expect $\mu_{n}$ to be contained in the span of $v_{1},\ldots,v_{r}$, we can also check that the projection of $\mu_{n}$ to that span preserves most of its energy.

Projecting $\widehat{P}_{\mathbf{X}^{(N)}}^{\mathrm{(MT)}}$ onto the span of $v_{1},\ldots,v_{r}$ yields

[TABLE]

This significantly reduces the estimation error $R$ since it will typically be almost orthogonal to $v_{l}$ and thus $\langle v_{l},R\rangle\approx 0$ for $l=1,\ldots,r$. The orthogonality follows from the fact that $R$ is approximately isotropically distributed (its coordinates are close to independent random variables), so its inner product with a small set of fixed vectors is small with high probability. Note that $\widehat{P}^{\mathrm{(MT,lin)}}$ is not guaranteed to be non-negative. To remedy this, negative components can be set to zero. However, this does not greatly affect the error of the estimate.

IV Numerical Results

To evaluate our method, we test it on noisy images simulated using (1) for $r=2$. The first noise source $\mathbf{Z}_{1}$ has its energy concentrated in the low frequencies, with a power spectrum of $P_{\mathbf{Z}_{1}}(\vec{\xi})=2\mathrm{rect}(4\|\vec{\xi}\|)$, where $\mathrm{rect}$ is the rectangular function with support $[-1/2,1/2]$. The second source $\mathbf{Z}_{2}$ has a power spectrum of $P_{\mathbf{Z}_{2}}(\vec{\xi})=1/(1+4\|\vec{\xi}\|)$. These are combined with normally distributed coefficients $a_{1},a_{2}\sim N(0,1)$ in (1) to yield $\mathbf{X}$. Figure 2(a) shows two sample images of size $N=32$.

We calculate noise images of different sizes $N$ to study the effect of this parameter on the accuracy of the estimation. Figure 2(b) shows the top eigenvalues of $\Sigma_{n}$ for a dataset of size $n=1024$ with image size $N=32$. There is a good separation of the top two eigenvalues from the bulk of the spectrum, with a spectral gap of $\lambda_{2}/\lambda_{3}$ of about $3.4$.

Figure 2(c,d) shows the mean absolute error (MAE) of $\widehat{P}^{\mathrm{(MT)}}_{\mathbf{X}^{(N)}_{1}},\ldots,\widehat{P}^{\mathrm{(MT)}}_{\mathbf{X}^{(N)}_{n}}$ and $\widehat{P}^{\mathrm{(MT,lin)}}_{\mathbf{X}^{(N)}_{1}},\ldots,\widehat{P}^{\mathrm{(MT,lin)}}_{\mathbf{X}^{(N)}_{n}}$ with respect to the conditional power spectra $P_{\mathbf{X}_{1}|\vec{a}_{1}},\ldots,P_{\mathbf{X}_{n}|\vec{a}_{n}}$ as functions of $n$ for image sizes $N=32$ and $N=64$, respectively. The unprojected multitaper estimator was computed with $W=1/16$ while the others used $W=1/64$. For comparison, we plot the baseline estimator of assigning $\mu_{n}$ to each power spectrum and the oracle projection of the multitaper estimate, where the top $r$ eigenvectors of the population covariance $\mathrm{Cov}[P_{\mathbf{X}|\vec{a}}]$ are used in the definition of $\widehat{P}^{\mathrm{(MT,lin)}}$. For small values of $n$, the covariance estimation fails, so the projected estimator $\widehat{P}^{\mathrm{(MT,lin)}}$ performs worse compared to the unprojected variant $\widehat{P}^{\mathrm{(MT)}}$. As $n$ increases, however, the linear projection improves results, eventually converging to the error obtained using the subspace obtained from the population covariance.

We also evaluate our algorithm using two experimental datasets from cryo-EM, one containing $n=10000$ images of size $N=130$ depicting a 70S ribosome [14] and another containing $n=105247$ images with $N=360$ depicting an 80S ribosome [18]. The top eigenvalues of $\Sigma_{n}$ are shown for the two datasets in Figure 3(a,b). As in the simulation, we see that a few eigenvalues dominate, with the estimation making up a separate bulk distribution. For the first dataset, a value of $r=2$ seems appropriate, while for the second dataset, we take $r=3$. This number of components is in line with most noise models for cryo-EM, which typically include two or three noise components [15, 16, 17].

Two sample conditional power spectra estimates for each dataset are shown in Figure 3(c,d). The projection images corresponding to the two power spectra in 3(c) are the ones shown previously in Figure 1. We see that the estimated power spectra are quite reasonable, with the lower-frequency noise image corresponding to the power spectrum concentrated in the low frequencies while the image with “whiter” noise has a flatter estimated power spectrum.

V Conclusion

We have introduced a factor analysis model for noise fields. These arise when the nature of the noise varies between measurements such as in cryo-EM, where the noise is affected by various experimental factors. To estimate the individual power spectra of each noise image we introduce a new estimation method which relies on approximating the covariance of these spectra. The method is shown to provide accurate results in both simulated and experimental datasets.

Acknowledgment

The authors thank Fred Sigworth for discussions about cryo-EM noise models, which inspired this work. They also thank Frederik Simons for his helpful comments.