Multi-target Detection with an Arbitrary Spacing Distribution

TL;DR

This paper introduces autocorrelation and EM-based methods for reconstructing signals in noisy multi-target detection scenarios with arbitrary spacing, especially effective at high noise levels where detection fails.

Contribution

It presents novel approaches that enable signal reconstruction without detection, applicable to arbitrary spacing distributions, and analyzes their sample complexity in high noise conditions.

Findings

Reconstruction feasible without detecting individual signals.

Sample complexity scales as 1/SNR^3 at low SNR.

Methods work with arbitrary spacing distributions.

Abstract

Motivated by the structure reconstruction problem in single-particle cryo-electron microscopy, we consider the multi-target detection model, where multiple copies of a target signal occur at unknown locations in a long measurement, further corrupted by additive Gaussian noise. At low noise levels, one can easily detect the signal occurrences and estimate the signal by averaging. However, in the presence of high noise, which is the focus of this paper, detection is impossible. Here, we propose two approaches---autocorrelation analysis and an approximate expectation maximization algorithm---to reconstruct the signal without the need to detect signal occurrences in the measurement. In particular, our methods apply to an arbitrary spacing distribution of signal occurrences. We demonstrate reconstructions with synthetic data and empirically show that the sample complexity of both methods scales as 1/SNR^3 in the low SNR regime.

Figures

Equations

$$y=s\ast x+\epsilon ,\epsilon \sim \mathcal{N}(0,{\sigma }^2{I}_N),$$

$$X\sim \{\begin{array}{l}\mathcal{N}({\theta }_0,{\sigma }^2{I}_L)\mathrm{if}\eta =0\\ \mathcal{N}({\theta }_1,{\sigma }^2{I}_L)\mathrm{if}\eta =1.\end{array}$$

$$\underset{\sigma \to \infty }{\lim }p(\hat{\eta }=\eta )\le \frac{1}{2},$$

$$a_z^q[l_1,l_2,\dots ,l_{q-1}]=\frac{1}{m}\sum _{i=0}^{m-1}z[i]z[i+l_1]\cdots z[i+l_{q-1}],$$

$$a_z^2[l]=a_z^2[-l],a_z^3[l_1,l_2]=a_z^3[l_2,l_1]=a_z^3[-l_1,l_2-l_1].$$

$${\mathbb{E}}_{\epsilon }\{a_y^1\}$$

$${\mathbb{E}}_{\epsilon }\{a_y^2[l]\}$$

$${\mathbb{E}}_{\epsilon }\{a_y^3[l_1,l_2]\}$$

$$+{\rho }_0{a}_x^1{\sigma }^2(\delta [l_1]+\delta [l_2]+\delta [l_1-l_2]),$$

$$a_y^1$$

$$a_y^2[l]$$

$$a_y^3[l_1,l_2]$$

$$f({\hat{\rho }}_0,\hat{x})$$

$$\displaystyle+\frac{2}{L(L+1)}\sum_{l_{2}=0}^{L-1}\sum_{l_{1}=0}^{l_{2}}\bigg{(}a_{y}^{3}[l_{1},l_{2}]-\hat{\rho}_{0}a_{\hat{x}}^{3}[l_{1},l_{2}]$$

$$\displaystyle\hskip 30.00005pt-\hat{\rho}_{0}a_{\hat{x}}^{1}\sigma^{2}(\delta[l_{1}]+\delta[l_{2}]+\delta[l_{1}-l_{2}])\bigg{)}^{2}.$$

$$\mathrm{RMSE}(\hat{x})=\frac{||\hat{x}-x|{|}_2}{||x|{|}_2},$$

$$\xi [l]=\frac{1}{M-1}\sum _{k=1}^{M-1}\delta [s_{k-1}+l-s_k],$$

$$a_y^1$$

$$a_y^2[l]$$

$$={\rho }_0{a}_x^2[l]+{\rho }_0\sum _{j=L}^{L+l-1}\xi [j]a_x^2[j-l]+{\sigma }^2\delta [l]+{\eta }_2^{\prime }[l]$$

$$a_y^3[l_1,$$

$$={\rho }_0{a}_x^3[l_1,l_2]+{\rho }_0\sum _{j=L}^{L+l_2-l_1-1}\xi [j]a_x^3[j-l_2,j+l_1-l_2]$$

$$\displaystyle\hskip 10.5pt+\rho_{0}\sum_{j=L}^{L+l_{1}-1}\xi[j]a_{x}^{3}[l_{2}-l_{1},j-l_{1}]+\rho_{0}a_{x}^{1}\sigma^{2}\big{(}\delta[l_{1}]$$

$$\displaystyle\hskip 10.5pt+\delta[l_{2}]+\delta[l_{1}-l_{2}]\big{)}+\eta_{3}^{\prime}[l_{1},l_{2}],$$

$${\rho }_1[i]={\rho }_0\xi [i+L],i=0,1,\dots ,L-2.$$

$$f(\hat{x},{\hat{\rho }}_0,{\hat{\rho }}_1)=(a_y^1-{\hat{\rho }}_0{a}_{\hat{x}}^1{)}^2$$

$$\displaystyle+\frac{1}{L}\sum_{l=0}^{L-1}\bigg{(}a_{y}^{2}[l]-\hat{\rho}_{0}a_{\hat{x}}^{2}[l]-\sigma^{2}\delta[l]-\sum_{j=L}^{L+l-1}\hat{\rho}_{1}[j-L]a_{\hat{x}}^{2}[j-l]\bigg{)}^{2}$$

$$\displaystyle+\frac{2}{L(L+1)}\sum_{l_{2}=0}^{L-1}\sum_{l_{1}=0}^{l_{2}}\bigg{(}a_{y}^{3}[l_{1},l_{2}]-\hat{\rho}_{0}a_{\hat{x}}^{3}[l_{1},l_{2}]-\hat{\rho}_{0}a_{\hat{x}}^{1}\sigma^{2}\big{(}\delta[l_{1}]$$

$$\displaystyle+\delta[l_{2}]+\delta[l_{1}-l_{2}]\big{)}-\sum_{j=L}^{L+l_{2}-l_{1}-1}\hat{\rho}_{1}[j-L]a_{\hat{x}}^{3}[j-l_{2},j+l_{1}-l_{2}]$$

$$\displaystyle-\sum_{j=L}^{L+l_{1}-1}\hat{\rho}_{1}[j-L]a_{\hat{x}}^{3}[l_{2}-l_{1},j-l_{1}]\bigg{)}^{2}.$$

$$x^{(n_{\max})}[l]=c_{0}+\sum_{n=1}^{n_{\max}}c_{n}\cos\bigg{(}\frac{2\pi nl}{L}\bigg{)}+d_{n}\sin\bigg{(}\frac{2\pi nl}{L}\bigg{)},$$

$$\Delta x=\{\begin{array}{ll}1& \text{if}{n}_{\max }=(L-1)/2,\\ L/2n_{\max }& \text{otherwise.}\end{array}$$

$$b_{x^{(n_{\max })}}^q[l_1,l_2,\dots ,l_{q-1}]$$

$$=a_{x^{(n_{\max })}}^q[\lfloor l_1\Delta x\rceil ,\lfloor l_2\Delta x\rceil ,\dots ,\lfloor l_{q-1}\Delta x\rceil ],$$

$$b_y^1$$

$$b_y^2[l]$$

$$b_y^3[l_1,l_2]$$

$$\displaystyle\hskip 30.00005pt-a_{y}^{1}\sigma^{2}\big{(}\delta[i_{1}]+\delta[i_{2}]+\delta[i_{1}-i_{2}]\big{)}\bigg{)},$$

$$\displaystyle b_{y}^{1}\$$

$$b_y^2[l]$$

$$b_y^3[l_1,l_2]\}$$

$$+\sum _{j=L^{\prime }}^{L^{\prime }+l_2-l_1-1}{\rho }_1^{(n_{\max })}[j-L^{\prime }]b_{x^{(n_{\max })}}^3[j-l_2,j+l_1-l_2]$$

Full text

I Introduction

We consider the multi-target detection (MTD) problem [10] to estimate a signal $x\in\mathbb{R}^{L}$ from a long, noisy measurement

[TABLE]

where $y\in\mathbb{R}^{N}$ is the linear convolution of an unknown binary sequence $s\in\{0,1\}^{N-L+1}$ with the signal, further corrupted by additive white Gaussian noise with zero mean and variance $\sigma^{2}$, and we assume $L\ll N$. Both $x$ and $s$ are treated as deterministic variables. The signal length $L$ and the noise variance $\sigma^{2}$ are assumed to be known. The non-zero entries of $s$ indicate the starting positions of the signal occurrences in $y$. We require the signal occurrences not to overlap, so consecutive non-zero entries of $s$ are separated by at least $L$ positions. Figure 1 gives an example of the measurement $y$ that contains three signal occurrences at different noise levels.

MTD belongs to the wider class of blind deconvolution problems [25], which have applications in astronomy [24, 33], microscopy [32, 28], system identification [1], and motion deblurring [26, 17], among others. The main difference is that we are only interested in estimating the signal $x$: we treat $s$ as a nuisance variable, while most of the literature on blind deconvolution aims to estimate both $x$ and $s$. This distinction allows us to estimate the signal at higher noise levels, which are usually not addressed in the literature on blind deconvolution, specifically because $s$ cannot be accurately estimated in such regimes. We give a theoretical argument for this claim below, and corroborate it with numerical experiments in Section IV.

This high-noise MTD model has been studied in [9, 10] under the assumption that the signal occurrences either are well separated or follow a Poisson distribution. In applications, however, the signal occurrences in the measurement might be arbitrarily close to each other. In this study, we extend the framework of [9, 10] to allow an arbitrary spacing distribution of signal occurrences by simultaneously estimating the signal and the distance distribution between consecutive signal occurrences.

The solution of the MTD problem is straightforward at high signal-to-noise ratio (SNR), such as the case shown in Figure 1(b). The signal can be estimated by first detecting the signal occurrences and then averaging. In the low SNR regime however, this simple method becomes problematic due to the difficulty of detection, as illustrated by Figure 1(c). More than difficulty, reliable detection becomes impossible beyond a critical noise level. This can be understood from [9, Proposition 3.1], which we reproduce here (the proof is based on the Neyman–Pearson lemma).

Proposition 1

Consider two known vectors $\theta_{0}=x$ and $\theta_{1}=0$ in $\mathbb{R}^{L}$ and a random variable $\eta$ taking value 0 or 1 with equal probability. We observe the random vector $X\in\mathbb{R}^{L}$ with the following distribution:

[TABLE]

For any deterministic estimator $\hat{\eta}$ of $\eta$,

[TABLE]

that is, the probability of success for the best deterministic estimator of $\eta$ is no better than a random guess in the limit of $\sigma\rightarrow\infty$.

In our context, this simple fact can be interpreted as follows: even if we know the true signal $x$ and an oracle provides windows of length $L$ in the observation $y$ which, with equal probability, contain either the noisy signal or pure noise, deciding whether a given window contains one or the other cannot be done significantly better than a random chance decision at high noise level. Results similar in spirit (but distinct in mathematical tools) are derived in [3].

In this work, we suggest to circumvent the estimation of $s$ by estimating the signal $x$ directly. Specifically, we propose two different approaches—autocorrelation analysis and an approximate version of the expectation maximization (EM) algorithm [18]. Since these methods do not require estimating $s$, Proposition 1 does not limit their performance in the low SNR regime.

Autocorrelation analysis relates the autocorrelations calculated from the noisy measurement to the signal. The signal is then estimated by fitting the autocorrelations through least squares. This approach is efficient as it requires only one pass over the data to calculate the autocorrelations; this is of particular importance as the data size grows.

In the second approach, the approximate EM algorithm, the signal is reconstructed by iteratively maximizing the data likelihood, which marginalizes over $s$; importantly, it does not estimate it explicitly. In contrast to autocorrelation analysis, the approximate EM algorithm scans through the whole dataset in each iteration, and hence requires much longer computational time.

In this study, we demonstrate the reconstruction of the underlying signal from the noisy measurement using the two proposed approaches. Our numerical experiments show that the approximate EM algorithm provides slightly more accurate estimates of the signal in the low SNR regime, whereas autocorrelation analysis is considerably faster, especially at low SNR. It is empirically shown that the sample complexity of both approaches scales as SNR*-3* at low SNR, with details discussed later in the text. In the high SNR regime, the sample complexity of both methods scales as SNR*-1*, the same as the sample complexity of the simple method that estimates the signal by first detecting the signal occurrences and then averaging.

The main contributions of this work are as follows.

We formalize MTD for an arbitrary signal spacing distribution, making it a more realistic model for applications. 2. 2.

We propose two algorithms to solve the MTD problem. In particular, our algorithm based on autocorrelations illustrates why, to recover $x$, we need not estimate all of $s$, but rather only a concise summary of it; and why it is possible, in principle at least, to solve this problem for arbitrary noise levels (given a sufficiently long observation). For our second algorithm, we note that the popular EM method is intractable, but we show how to implement an approximation of it, which performs well in practice. 3. 3.

For both algorithms, we design a coarse-to-fine multi-resolution scheme to alleviate issues pertaining to non-convexity. This is related to the ideas of frequency marching which are often used in cryo-electron microscopy (cryo-EM) [6, 34].

II Autocorrelation analysis

In what follows, we discuss autocorrelations of both the signal $x$ (of length $L$) and the measurement $y$ (of length $N$). To keep notation general, we here consider a sequence $z$ of length $m$, and define its autocorrelations of order $q=1,2,\dots$ for any integer shifts $l_{1},l_{2},\dots,l_{q-1}$ as

[TABLE]

where $z$ is zero-padded for indices out of the range $[0,m-1]$. We have $m=L$ when $z$ represents the signal $x$ and $m=N$ when $z$ represents the measurement $y$. Since the autocorrelations only depend on the differences of the integer shifts and are invariant under any permutation of the shifts, for the second- and third-order autocorrelations we have the symmetries

[TABLE]

For applications of higher-order autocorrelations and their symmetries, see for example [22, 35, 31, 4].

In this section, we describe how the autocorrelations of the noisy measurement $y$ are related to the underlying signal $x$. These relations are later used to estimate $x$ without the need to identify the locations of signal occurrences, which are nuisance variables and difficult, if not impossible, to determine reliably at high noise levels. For completeness, we include a brief discussion of the special case where the signals are well separated and refer the reader to [9, 10] for details. The generalization to an arbitrary signal spacing distribution follows.

II-A Well-separated signals

The signals are said to be well-separated when the consecutive non-zero entries of $s$ are separated by at least $2L-1$ positions. Under this condition, the autocorrelations of $y$ with integer shifts $l_{1},l_{2},\dots,l_{q-1}$ within the range $[-(L-1),L-1]$ are unaffected by the relative positions of signal occurrences. As a result, these autocorrelations of $y$ provide fairly direct information about those of $x$, and therefore about $x$ itself.

Due to the presence of Gaussian noise, the entries of $y$ are stochastic. Taking the expectations of the first three autocorrelations of $y$ with respect to the distribution of Gaussian noise, we obtain these relations (see Appendix A):

[TABLE]

where $0\leq l<L$ and $0\leq l_{1}\leq l_{2}<L$. Here, $\rho_{0}=ML/N$ denotes the signal density, where $M$ is the number of signal occurrences in $y$. We assume $M$ grows with $N$ at the constant rate $\rho_{0}/L$. The delta functions, defined by $\delta[0]=1$ and $\delta[l\neq 0]=0$, are due to the autocorrelations of the white Gaussian noise. As indicated by the studies in phase retrieval [15, 7, 12], in general a 1D signal $x$ cannot be uniquely determined by its first two autocorrelations. It is thus necessary (and generically sufficient [10]) to include the third-order autocorrelations to uniquely determine $x$.

The expectations of the autocorrelations of $y$ are estimated by averaging over the given noisy measurement. This data reduction requires only one pass over the data, which is a great computational advantage as the data size grows. As shown in Appendix A, for a given signal $x$, the average over the noisy entries of $y$ gives the relations:

[TABLE]

where $\eta^{1},\eta^{2}[l]$ and $\eta^{3}[l_{1},l_{2}]$ are random variables with zero mean and variances \mathcal{O}\big{(}\sigma^{2}/N\big{)}, \mathcal{O}\big{(}(\sigma^{2}+\sigma^{4})/N\big{)} and \mathcal{O}\big{(}(\sigma^{2}+\sigma^{6})/N\big{)} respectively. For autocorrelations of order $q$, the standard deviations scale as $\sigma^{q}/\sqrt{N}$ at high noise levels. Therefore, we need $\sqrt{N}/\sigma^{q}\gg 1$ in order for the $a_{y}^{q}$ calculated from the noisy measurement to be a good estimator for $\mathop{\mathbb{E}_{\varepsilon}}\{a_{y}^{q}\}$, and thus to establish a reliable relation with the signal $x$ such as (3)-(II-A). Since the SNR is proportional to $\sigma^{-2}$, the sample complexity therefore scales as SNR*-q*. We also expect the error of the reconstructed signal to depend on the errors of the highest-order autocorrelations used in the analysis at high noise levels.

We estimate the signal density $\rho_{0}$ and signal $x$ by fitting the first three autocorrelations of $y$ via non-linear least squares:

[TABLE]

The weights are chosen as the inverse of the number of terms in the respective sums. Since the autocorrelations have symmetries, as indicated in (2), the summations above are restricted to the non-redundant shifts. Due to the errors in estimating $\mathop{\mathbb{E}_{\varepsilon}}\{a_{y}^{3}\}$ with $a_{y}^{3}$, we expect the root-mean-square error of the reconstructed signal $\hat{x}$,

[TABLE]

to scale as $\sigma^{3}$ in the low SNR regime.

We mention that there exist alternative methods to fit the observed autocorrelations. One possibility is to reformulate the problem as a tensor sensing problem: the autocorrelations are then linear functions of an $(L+1)\times(L+1)\times(L+1)$ tensor. Such an approach has been proven to be effective for the related problems of multi-reference alignment [5, 29] and phase retrieval [16]. Similar formulations are also found in tensor optimization [2, 30]. Nevertheless, this method requires lifting the original problem to a higher-dimensional space, which increases the computational burden and scales poorly with the dimension of the problem. In contrast, our least-squares formulation operates in the ambient dimension of the problem and thus might be applicable to high-dimensional setups.

II-B Arbitrary spacing distribution

The condition of well-separated signals can be further relaxed to allow arbitrary spacing distribution by assuming that the signal occurrences in any subset of $y$ follow the same spatial distribution. To this end, we define the pair separation function as follows.

Definition 1

For a given binary sequence $s$ identifying $M$ starting positions of signal occurrences (that is, $\sum_{i}s[i]=M$), the pair separation function $\xi[l]$ is defined as

[TABLE]

where $s_{k}$ indicates the index of the $k^{\mathrm{th}}$ non-zero entry of the sequence $s$. It is not hard to see that $\sum_{l}\xi[l]=1$.

In particular, we force $\xi[0]=\cdots=\xi[L-1]=0$ since we exclude overlapping occurrences; then $\xi[L]$ is the fraction of pairs of consecutive signal occurrences that occur right next to each other (no spacing at all), $\xi[L+1]$ is the fraction of pairs of consecutive signal occurrences that occur with one signal-free entry in between them, etc.

In contrast to the well-separated model, autocorrelations of $y$ may now involve correlating distinct occurrences of $x$, which may be in various relative positions. The crucial observation is that these autocorrelations depend only indirectly on $s$, namely, through $L-1$ entries of the unknown $\xi$, which have much smaller dimensions than $s$. As shown in Appendix A, for a given signal $x$, the autocorrelations of the given measurement $y$ gives the relations:

[TABLE]

where $0\leq l<L$ and $0\leq l_{1}\leq l_{2}<L$. The random variables $\eta_{1}^{\prime}$, $\eta_{2}^{\prime}[l]$ and $\eta_{3}^{\prime}[l_{1},l_{2}]$ have zero mean and variances \mathcal{O}\big{(}\sigma^{2}/N\big{)}, \mathcal{O}\big{(}(\sigma^{2}+\sigma^{4})/N\big{)} and \mathcal{O}\big{(}(\sigma^{2}+\sigma^{6})/N\big{)} respectively. Note that (8)–(II-B) reduce to (3)–(II-A) when $\xi[L]=\xi[L+1]=\cdots=\xi[2L-2]=0$, as required by the condition of well-separated signals. Expressions (8)–(II-B) further simplify upon defining

[TABLE]

After calculating the first three autocorrelations of $y$ from the noisy measurement, we estimate the signal $x$ and the parameters $\rho_{0}$ and $\rho_{1}$ by fitting the autocorrelations of $y$ through the non-linear least squares

[TABLE]

Numerically, we find that this problem can often be solved, meaning that, even though $s$ cannot be estimated, we can still estimate $x$ and the summarizing statistics $\rho_{0}$ and $\rho_{1}$. As discussed in Section II-A, the RMSE of the reconstructed signal is expected to scale as $\sigma^{3}$ in the low SNR regime owing to the errors in estimating $\mathop{\mathbb{E}_{\varepsilon}}\{a_{y}^{3}\}$ with $a_{y}^{3}$.

II-C Frequency marching

To minimize the least squares in (II-A) and (II-B), we use the trust-regions method in Manopt [13] over the product of the Euclidean manifold with the constraints that $\hat{\rho}_{0}$ and $\hat{\rho}_{1}$ are positive. However, as the least squares problems are inherently non-convex, we observe that the iterates of the trust-regions method used for minimization are liable to stagnate in local minima. To alleviate this issue, we adopt the frequency marching scheme [6] in our optimization. The idea behind frequency marching is based on the following heuristics:

The coarse-grained (low-resolution) version of the original problem has a smoother optimization landscape so that it is empirically more likely for the iterates of the optimization algorithm to reach the global optimum of the coarse-grained problem. 2. 2.

Intuitively, the global optimum of the original problem can be reached more easily by following the path along the global optima of a series of coarse-grained problems with incremental resolution.

Our goal is to guide the iterates of the optimization algorithm to reach the global optimum of the original problem by successively solving the coarse-grained problems, which are warm-started with the solution from the previous stage.

The coarse-grained problems are characterized by the order of the Fourier series, $n_{\max}=1,2,\dots,\lfloor L/2\rfloor$, used to express the low-resolution approximate $x^{(n_{\max})}$ by

[TABLE]

where $l=0,1,\dots,L-1$. Instead of the entries of $x$, the least squares are minimized with respect to the Fourier coefficients in our frequency marching scheme. The order $n_{\max}$ is related to the spatial resolution by Nyquist rate:

[TABLE]

The spatial resolution subdivides the signal $x$ into $L^{\prime}=L/\Delta x$ units and represents the “step size” of the shifts for the coarse-grained autocorrelations.

We define the coarse-grained autocorrelations of $x^{(n_{\max})}$ of order $q$ for integer shifts $l_{1},l_{2},\dots,l_{q-1}$ as

[TABLE]

where $\lfloor\cdot\rceil$ rounds the argument to the nearest integer. The coarse-grained autocorrelations of $y$ are given by sub-sampling the original autocorrelations calculated from the full measurement. With $b[l]$ denoting the bin centered at $l\Delta x$, where $l=0,1,\dots,L^{\prime}-1$, we estimate the coarse-grained autocorrelations of $y$ by

[TABLE]

where $0\leq l<L^{\prime}$, $0\leq l_{1}\leq l_{2}<L^{\prime}$, and $B_{2}$ and $B_{3}$ represent the number of terms in the respective sums.

Following the discussion in Section II-B, we relate the first three autocorrelations of $y$ to those of $x^{{(n_{\max})}}$, as defined in (13), by

[TABLE]

The autocorrelations are related by approximation instead of equality to reflect the errors due to the low-resolution approximation. Above, we define $\rho_{1}^{(n_{\max})}$ as the product of the signal density $\rho_{0}$ and the coarse-grained pair separation function:

[TABLE]

where $\xi^{(n_{\max})}$ is defined as

[TABLE]

In each stage of frequency marching, we estimate the Fourier coefficients, which are related to $x^{(n_{\max})}$ through (13), and the parameters $\rho_{0},\rho_{1}^{(n_{\max})}$ by fitting the coarse-grained autocorrelations of $y$ through the non-linear least squares

[TABLE]

Our frequency marching scheme increments the order $n_{\max}$ from $1$ to $\lfloor L/2\rfloor$, and the computed solution of each stage is used to initialize optimization in the next stage.

III Expectation maximization

In this section, as an alternative to the autocorrelations approach, we describe an approximate EM algorithm to address both the cases of well-separated signals and arbitrary spacing distribution. A frequency marching scheme is also designed to help the iterates of the EM algorithm converge to the global maximum of the data likelihood.

III-A Well-separated signals

Given the measurement $y$ that follows the MTD model (1), the maximum marginal likelihood estimator (MMLE) for the signal $x$ is the maximizer of the likelihood function $p(y|x)$. Within the EM framework [18], the nuisance variable $s$ is treated as a random variable drawn from some distribution under the condition of non-overlapping signal occurrences. The EM algorithm estimates the MMLE by iteratively applying the expectation (E) and maximization (M) steps. Specifically, given the current signal estimate $x_{k}$, the E-step constructs the expected log-likelihood function

[TABLE]

where the summation runs over all admissible configurations of the binary sequence $s$. The signal estimate is then updated in the M-step by maximizing $Q(x|x_{k})$ with respect to $x$. The major drawback of this approach is that the number of admissible configurations for $s$ grows exponentially with the problem size. Therefore, the direct application of the EM algorithm is computationally intractable, even for very short measurements.

In our framework of the approximate EM algorithm, we first partition the measurement $y$ into $N_{d}=N/L$ non-overlapping segments, each of length $L$, and denote the $m^{\mathrm{th}}$ segment by $y_{m}$. Overall, the signal can occur in $2L-1$ different ways when it is present in a segment. The signal is estimated by the maximizer of the approximate likelihood function

[TABLE]

where we ignore the dependencies between segments. Our approximate EM algorithm works by applying the EM algorithm to estimate the MMLE of (16), without any prior on the signal. As we will see in Section IV, the validity of the approximation is corroborated by the results of our numerical experiments.

Depending on the position of signal occurrences, the segment $y_{m}$ can be modeled by

[TABLE]

Here, $Z$ first zero-pads $L$ entries to the left of $x$, and $R_{l_{m}}$ circularly shifts the zero-padded sequence by $l_{m}$ positions, that is,

[TABLE]

where $l_{m}=0,1,\dots,2L-1$ and is treated as a random variable. The operator $C$ then crops the first $L$ entries of the circularly shifted sequence, which are further corrupted by additive white Gaussian noise. In this generative model, $l_{m}=0$ represents no signal occurrence in $y_{m}$, and $l_{m}=1,\dots,2L-1$ enumerate the $2L-1$ different ways a signal can appear in a segment.

In the E-step, our algorithm constructs the expected log-likelihood function

[TABLE]

given the current signal estimate $x_{k}$, where

[TABLE]

with the normalization $\sum_{l=0}^{2L-1}p(y_{m}|l,x)=1$. From Bayes’ rule, we have

[TABLE]

which is the normalized likelihood function $p(y_{m}|l,x_{k})$, weighted by the prior distribution $p(l|x_{k})$. In general, the prior distribution $p(l|x_{k})$ is independent of the model $x_{k}$ and can be estimated simultaneously with the signal.

Denoting the prior distribution $p(l)$ by $\alpha[l]$, we rewrite (17) as (up to an irrelevant constant)

[TABLE]

We note that the dependence of $Q(x,\alpha|x_{k},\alpha_{k})$ on the current prior estimates $\alpha_{k}$ lies in $p(l|y_{m},x_{k})$ through (19). The M-step updates the signal estimate and the priors by maximizing $Q(x,\alpha|x_{k},\alpha_{k})$ under the constraint that the priors lie on the simplex $\Delta_{2L}$:

[TABLE]

As shown in Appendix B, we obtain the update rules

[TABLE]

where $0\leq j<L$, and

[TABLE]

where $0\leq l<2L$. We repeat the iterations of the EM algorithm until the estimates stop improving, as judged within some tolerance.

III-B Arbitrary spacing distribution

We extend the approximate EM approach to address arbitrary spacing distribution of signal occurrences by reformulating the probability model: each segment $y_{m}$ can now contain up to two signal occurrences. In this case, the two signals can appear in $L(L-1)/2$ different combinations in a segment, which is explicitly modeled by

[TABLE]

where $L<l^{m}_{1}<2L$ and $0<l^{m}_{2}\leq l^{m}_{1}-L$. Given the signal estimate $x_{k}$ and the shifts $l^{m}_{1}=l_{1},l^{m}_{2}=l_{2}$, the likelihood function $p(y_{m}|l_{1},l_{2},x_{k})$ can be written as

[TABLE]

and we have the normalization condition

[TABLE]

Incorporating the terms with two signal occurrences, the E-step constructs the expected log-likelihood function as

[TABLE]

where the prior $p(l_{1},l_{2})$ is denoted by $\alpha[l_{1},l_{2}]$.

Under the assumption that the signal occurrences in any subset of $y$ follow the same spatial distribution, we can parametrize the priors with the pair separation function $\xi$ (see Definition 1). Recall that $(M-1)\xi[l]$ is the number of pairs of consecutive signal occurrences whose starting positions are separated by exactly $l$ positions. The priors $\alpha[l_{1},l_{2}]$ can be related to the probability that two signal occurrences appear in the combination specified by $(l_{1},l_{2})$ in a segment of length $L$ selected from the measurement $y$, which is estimated by

[TABLE]

where $L<l_{1}<2L$ and $0<l_{2}\leq l_{1}-L$. Here, $(M-1)\xi[l_{1}-l_{2}]$ is the number of segments that realize the configuration of signal occurrences specified by $(l_{1},l_{2})$, and $N-L+1$ indicates the total number of segment choices. The priors $\alpha[l]$ can similarly be related to the probability that a signal occurs in the way specified by $l$ in a segment of length $L$:

[TABLE]

where $0<l\leq L$. An interesting observation is that the number of signal occurrences $M$ can be estimated by

[TABLE]

since the signal occurrences are required not to overlap. The value of the prior $\alpha[0]$ is determined by the normalization

[TABLE]

From (25), (III-B) and (27), we see that the priors are uniquely specified by the positive parameters $\alpha[0]$, $\alpha[1]$, $\rho_{1}[0]$, $\rho_{1}[1]$, $\dots,\rho_{1}[L-2]$, with $\rho_{1}$ defined in (11). Therefore, the normalization (27) can be rewritten as

[TABLE]

In the special case of well-separated signals, where $\xi[i]=0$ for $L\leq i\leq 2L-2$, we have $\alpha[1]=\alpha[2]=\cdots=\alpha[2L-1]$ and the normalization $\alpha[0]+(2L-1)\alpha[1]=1$.

In the M-step, we update the signal estimate and the parameters $\alpha[0],\alpha[1],\rho_{1}$ by maximizing $Q(x,\alpha|x_{k},\alpha_{k})$ under the constraint that (28) is satisfied:

[TABLE]

As shown in (III-B), $Q(x,\alpha|x_{k},\alpha_{k})$ is additively separable for $x$ and $\alpha$, so the constrained maximization (III-B) can be reached by maximizing $Q(x,\alpha|x_{k},\alpha_{k})$ with respect to $x$ and the parameters $\alpha[0],\alpha[1],\rho_{1}$ separately. Maximizing $Q(x,\alpha|x_{k},\alpha_{k})$ with respect to $x$ yields the update rule (see Appendix C for derivation)

[TABLE]

where $0\leq j<L$. To update the parameters $\alpha[0],\alpha[1],\rho_{1}$, we note that the function $Q(x,\alpha|x_{k},\alpha_{k})$ is concave with respect to these parameters and the constraint (28) forms a compact convex set. Therefore, the constrained maximization is achieved using the Frank-Wolfe algorithm [20].

III-C Frequency marching

Because the iterates of the EM algorithm are not guaranteed to converge to the maximum likelihood, we develop a frequency marching scheme to help the iterates converge to the global maximum. Recall that frequency marching converts the original optimization problem into a series of coarse-grained problems with gradually increasing resolution. The coarse-grained version of the EM algorithm is characterized by the spatial resolution

[TABLE]

where $n_{\max}=1,2,\dots,\lfloor(L+1)/2\rfloor$. This spatial resolution models the signal shifts in the coarse-grained problem, with the corresponding expected likelihood function given by

[TABLE]

where $L^{\prime}=L/\Delta_{x}$. This likelihood function is equal to the original one shown in (III-B) restricted to the rounded shifts, so the signal can be updated with (III-B) by ignoring the terms with irrelevant shifts.

Mimicking (25) and (III-B), we construct the expressions for the priors in the coarse-grained problem as

[TABLE]

where $L^{\prime}<l_{1}<2L^{\prime}$, $0<l_{2}\leq l_{1}-L^{\prime}$, $0<l\leq L^{\prime}$, and $\xi^{(n_{\max})}$ is defined in (15). Since the priors are uniquely specified by the positive parameters $\alpha^{(n_{\max})}[0],$ $\alpha^{(n_{\max})}[1],$ $\rho^{(n_{\max})}_{1}[0],$ $\rho^{(n_{\max})}_{1}[1],$ $\dots,\rho^{(n_{\max})}_{1}[L^{\prime}-2],$ with $\rho_{1}^{(n_{\max})}$ defined in (14), we can update the priors by maximizing $Q^{(n_{\max})}(x,\alpha^{(n_{\max})}|x_{k},\alpha^{(n_{\max})}_{k})$ with respect to these parameters under the constraint that the parameters lie on the simplex defined by

[TABLE]

Incrementing from $n_{\max}=1$ to $n_{\max}=\lfloor(L+1)/2\rfloor$, the estimated signal and priors in each stage of frequency marching are used to initialize the optimization problem in the next stage.

IV Numerical Experiments

This section describes the construction of our synthetic data and the performance of the two proposed methods.111The code for all experiments is publicly available at https://github.com/tl578/multi-target-detection. A measurement is generated by first sampling $M$ integers from $[0,N-L]$, with the constraint that any two samples differ by at least $L+W$ entries, for some integer $W$. We use $W=L-1$ for the case of well-separated signals, while $W=0$ is chosen to test our methods for arbitrary spacing distribution of signal occurrences. The sampling is done by generating the random integers uniformly one by one, rejecting any integer that would violate the constraint with respect to previously accepted samples. The $M$ integers indicate the starting positions of the signal occurrences, which are recorded as the non-zero entries of the binary sequence $s$. The noisy measurement $y$ is given by the linear convolution of $s$ with the signal $x$, with each entry of $y$ further corrupted by additive white Gaussian noise of zero mean and variance $\sigma^{2}$. Our synthetic data all have length $N=10^{6}$, and are constructed with the signal shown in Figure 2(a), which has length $L=10$. Also shown in Figures 2(b) and 2(c) are examples of the corrupted signals by different levels of noise. The number of signal occurrences $M$ is adjusted so that the signal density $\rho_{0}=ML/N$ is equal to 0.3 and 0.5 for the cases of well-separated signals and arbitrary spacing distribution respectively. Figure 3 shows the resulting pair separation function for our test case of arbitrary spacing distribution.

For both cases of signal spacing distribution, we reconstruct the signal from synthetic data using the two proposed methods with the frequency marching scheme over a wide range of $\sigma$. When the frequency marching scheme is not implemented, we observe that the algorithms usually suffer from large errors due to the non-convexity of the problems. A total of 20 instances of synthetic data are generated for each value of $\sigma$, and each instance is solved by the two methods using the frequency marching schemes. For our autocorrelation analysis methods, we run the optimization with 10 different random initializations and retain the one that minimizes the cost function, as defined in (II-A) or (II-B). For the EM approach, we also use 10 random initializations and retain the solution with the maximum data likelihood. Given the signal estimate $x_{k}$, we approximate the data likelihood for the cases of well-separated signals and arbitrary spacing distribution by

[TABLE]

and

[TABLE]

respectively. This is the same approximation we use to formulate the EM algorithm.

The RMSEs of the reconstructions for the case of well-separated signals are shown in Figure 4. As predicted in Section II-A, the RMSE for autocorrelation analysis scales in proportion to $\sigma$ at high SNR, and to $\sigma^{3}$ at low SNR. Interestingly, the EM algorithm exhibits the same behavior, which empirically shows that the two approaches share the same scaling of sample complexity in the two extremes of noise levels. Figure 5(a) and 5(b) show an instance of the signals reconstructed from data generated with $\sigma=10^{0.3}$ and $\sigma=10^{0.5}$ respectively. The same scaling is not observed for the reconstructed signal density $\rho_{0}$ (not shown), although the relative errors are generally well below a few percents even in the noisiest cases.

Figure 6 shows the RMSEs of the reconstructed signal for an arbitrary spacing distribution displayed in Figure 3. As a comparison, we also run the estimation algorithms under the (incorrect) assumption of well-separated signals on the same datasets. We first see that autocorrelation analysis assuming well-separated signals obtains poor reconstruction at all noise levels. Although the EM approach that assumes well-separated signals produces more accurate estimates, the nearly constant RMSEs at high SNR indicate the systematic error due to model misspecification. By contrast, the algorithms that assume an arbitrary spacing distribution achieve much better reconstruction, and the resulting RMSEs have the same scaling behaviors at the two extremes of SNR as shown in Figure 4. An instance of the signals reconstructed from data generated with $\sigma=10^{0.3}$ and $\sigma=10^{0.5}$ are shown in Figure 5(c) and 5(d) respectively.

One of the premises of this paper is that, at high noise levels, detection-based methods are destined to fail because detection cannot be done reliably. To support this theoretical argument further, in Figures 4 and 6 we display the performance of the following oracle, named ‘deconv’. This method is given a strictly simpler task to solve: it must estimate $x$ given the observation $y$, the noise level $\sigma$ and the number of signal occurrences $M$, as well as the $\ell_{2}$-distance between the ground truth $x$ and every single window of length $L$ in the observation $y$. Precisely, the oracle has access to the vector $z\in\mathbb{R}^{N-L+1}$ defined for $i=0,\ldots,N-L$ by

[TABLE]

The vector $z$ provides this estimator the advantage of finding the likely locations of signal occurrences. If the norm of $x$ is known, this oracle advantage is equivalent to knowledge of the cross-correlation between $x$ and $y$. Then, the oracle-based estimator ‘deconv’ proceeds as follows: it selects the index $i$ corresponding to the lowest value in $z$, excludes other indices too close to $i$ as prescribed by the signal separation model under consideration, and repeats this procedure $M$ times to produce an estimator of $s$: the $M$ starting positions of signal occurrences in $y$. The signal $x$ is then estimated by averaging the $M$ selected segments of $y$, effectively deconvolving $y$ by the estimator of $s$. To our point: despite the unfair oracle advantage, this estimator is unable to estimate $x$ to a competitive accuracy when the noise level is large.

Also plotted in Figures 4 and 6 are the RMSEs of the estimated signal when the binary sequence $s$ is known. In this case, the RMSEs scale as $\sigma$ and serve as the lower bounds of any algorithm. We can see that, at the low noise levels, our approximate EM algorithm performs nearly as well as the case when $s$ is known.

Figure 7 shows the RMSEs of the reconstructed values of $\rho_{1}$ for the experiments “AA-asd” and “EM-asd” shown in Figure 6. Although our methods are not able to reconstruct $\rho_{1}$ to high precision, the results shown in Figure 6 indicate the necessity to include the pair separation function in the model to achieve good signal reconstruction.

Figure 8 shows the average computation time for the main methods in Figures 4 and 6. EM is slower because it cross-correlates all the observed segments with the signal estimate in each iteration. This issue becomes especially prominent in the case of arbitrary spacing distribution. On the other hand, autocorrelation analysis requires shorter computation time by summarizing the data as autocorrelation statistics with one pass over the data. We note that the distinct computational speeds for autocorrelation analysis and the EM algorithm is also observed in the related problem of multi-reference alignment [8, 14]. The difference in run time and the similar reconstruction quality of the two methods at low SNR make autocorrelation analysis the preferred approach for large datasets.

V Discussion

In this study, we have presented two approaches—autocorrelation analysis and an approximate EM algorithm—to tackle the MTD problem without the need to detect the positions of the underlying signal occurrences. By incorporating the notion of pair separation function, we generalize the solution of the MTD problem from the special case of well-separated signals to allow an arbitrary spacing distribution of signal occurrence. It is empirically shown that the two methods have the same scaling of sample complexity in the two extremes of noise levels, in particular, SNR*-3* at high noise level. Since the optimizations in both methods are non-convex, computing schemes based on frequency marching are designed to help the iterates converge to the global optimum.

Our study of the MTD model is primarily motivated by the goal to reconstruct the structures of small biomolecules using cryo-EM [21, 11]. In a cryo-EM experiment, individual copies of the target biomolecule are dispersed at unknown 2D locations and 3D orientations in a thin layer of vitreous ice, from which 2D tomographic projection images are produced by an electron microscope. To minimize the irreversible structural damage, it is necessary to keep the electron dose low. As a result, the projection images are considerably noisy, and high-resolution structure estimation requires averaging over a large number of noisy projections. In particular, as the size of the molecule gets smaller, the SNR of the data drops correspondingly.

Currently, the analysis workflow of cryo-EM data is roughly divided into two steps: The first step, known as particle picking, detects the locations of the biomolecules in the noisy projection images. The 3D structure of the biomolecule is reconstructed in the second step from the unoriented particle projections. When the sizes of biomolecules get smaller, however, reliable detection of their positions becomes challenging, which in turn hampers successful particle picking [3, 9]. It is estimated from first principles that, in order to obtain a 3 Å resolution reconstruction, the particle size should be at least 45 Å so that the particle occurrences can be accurately detected [23]. The results of [9, 10] and this work suggest that it is possible to bypass the need to detect particle positions but still reconstruct the structure at high resolution.

A key feature of cryo-EM that the MTD model fails to capture is that the “signals” are actually 2D projections of the underlying biomolecules at random 3D orientations, and the orientation distribution is usually non-uniform and unknown. This distinction makes the direct application of our approximate EM algorithm prohibitive, because each observed window now needs to cross-correlate with model projections with all possible 2D in-plane translations and 3D orientations, let alone the case where multiple copies of biomolecules may appear in a window. Further approximations such as representing the data and model projections in low-dimensional subspaces [19, 27] seem necessary to push this approach forward. We note that the use of frequency marching in the MTD problem is also related to this idea: the signal is mapped to low-dimensional subspaces with growing dimensions specified by $n_{\max}$.

On the other hand, the autocorrelations calculated from the noisy measurements in cryo-EM involve averages over both the 3D orientations and 2D in-plane translations. It was pointed out that the 3D structure reconstruction problem by fitting the first three autocorrelations of a biomolecule might be ill-conditioned [9], so information from higher-order autocorrelations seems required. We also expect the cross-terms between 2D projections of different copies to further complicate the reconstruction problem. The resolution of these technical challenges may help extend the use of cryo-EM to smaller biomolecules that are currently believed insurmountable, and is the subject of our ongoing studies.

Acknowledgment

We would like to thank Ayelet Heimowitz, Joe Kileel, Eitan Levin and Amit Moscovich for stimulating discussions and the anonymous reviewers for their valuable comments.

Appendix A Derivations of Relations (3)–(II-A) and (8)–(II-B)

With the noisy measurement $y$ defined in (1), we express its first order autocorrelation as

[TABLE]

Taking the expectation of $a_{y}^{1}$ with respect to the distribution of Gaussian noise, we obtain

[TABLE]

with variance $\mathcal{O}(\sigma^{2}/N)$. We note that this expression applies to both the cases of well-separated signals and arbitrary spacing distribution.

The second order autocorrelation of $y$ with shifts $l=0,1,\dots,L-1$ can be written as

[TABLE]

where $a^{2}_{s\ast x}[l]$ denotes the second order autocorrelation of the linear convolution $s\ast x$ with shifts $l$, and $s_{k}$ indicates the index of the $k^{\mathrm{th}}$ non-zero entry of the sequence $s$. Taking the expectation of $a_{y}^{2}[l]$ with respect to the distribution of Gaussian noise, we have

[TABLE]

The term that scales as \mathcal{O}\big{(}\varepsilon^{2}\big{)} in (A) contributes to the variance by \mathcal{O}\big{(}\sigma^{4}/N\big{)}, while the terms that scale as \mathcal{O}\big{(}\varepsilon\big{)} contribute to the variance by \mathcal{O}\big{(}M\sigma^{2}/N^{2}\big{)}. Since we assume that $M$ grows with $N$ at a constant rate, the resulting variance scales as \mathcal{O}\big{(}(\sigma^{2}+\sigma^{4})/N\big{)}.

For well-separated signals,

[TABLE]

and we obtain the relation in (4). For arbitrary spacing distribution, we need to include the correlation terms between consecutive signal occurrences:

[TABLE]

which gives the relation in (II-B). The quantity $\xi[j]$ denotes the number of consecutive 1’s that are separated by exactly $j$ entries in the sequence $s$, divided by $M-1$ so that $\sum_{j}\xi[j]=1$.

As for the third order autocorrelation of $y$ with shifts $0\leq l_{1}\leq l_{2}<L$, we have

[TABLE]

where $a^{3}_{s\ast x}[l_{1},l_{2}]$ denotes the third order autocorrelation of the linear convolution $s\ast x$ with shifts $(l_{1},l_{2})$, and $x$ is zero-padded for indices out of the range $[0,L-1]$.

Taking the expectation of $a_{y}^{3}[l_{1},l_{2}]$ with respect to the distribution of Gaussian noise, we obtain

[TABLE]

As for the variance, the term that scales as \mathcal{O}\big{(}\varepsilon^{3}\big{)} in (A) contributes to the variance by \mathcal{O}\big{(}\sigma^{6}/N\big{)}, while the terms that scale as \mathcal{O}\big{(}\varepsilon\big{)} contribute to the variance by \mathcal{O}\big{(}M\sigma^{2}/N^{2}\big{)}. Since we assume that $M$ grows with $N$ at a constant rate, the resulting variance scales as \mathcal{O}\big{(}(\sigma^{2}+\sigma^{6})/N\big{)}.

For well-separated signals, we have

[TABLE]

and we obtain the relation in (II-A). By including the correlations between consecutive signal occurrences, we expand $a^{3}_{s\ast x}[l_{1},l_{2}]$ for the case of arbitrary spacing distribution as

[TABLE]

which gives the relation in (II-B).

Appendix B Derivation of (21) and (22)

The constrained maximization in (20) can be achieved with the unconstrained maximization of the Lagrangian

[TABLE]

where $\lambda$ denotes the Lagrange multiplier. We note that the constraints in (20) involve the inequalities that the priors are non-negative. Such constrained maximization in general cannot be achieved by maximizing the Lagrangian, for the inequalities might be violated. As we will see later, however, these inequalities are automatically satisfied at the computed maximum (or local maximum) of the Lagrangian, which justifies this approach.

Since $Q(x,\alpha|x_{k},\alpha_{k})$ is additively separable for $x$ and $\alpha$, we maximize $\mathcal{L}(x,\alpha,\lambda)$ with respect to $x$ and $\alpha$ separately. At the maximum of $\mathcal{L}(x,\alpha,\lambda)$, we have

[TABLE]

where $0\leq j<L$. With $p(y_{m}|l,x)$ give by (18), we can write

[TABLE]

Substituting this expression into (33), we obtain

[TABLE]

Rearranging the terms gives the update rule shown in (21).

In order to update the priors, we maximize $\mathcal{L}(x,\alpha,\lambda)$ with respect to $\alpha$:

[TABLE]

where $0\leq l<2L$. We thus obtain the update rule for $\alpha$ as

[TABLE]

and we can immediately solve $\lambda=N_{d}$ from the normalization $\sum_{l=0}^{2L-1}\alpha[l]=1$.

Appendix C Derivation of (III-B)

As shown in (III-B), $Q(x,\alpha|x_{k},\alpha_{k})$ is additively separable for $x$ and $\alpha$. Moreover, the constraint (28) only involves the values of the parameters $\alpha[0],\alpha[1],\rho_{1}$, so the maximization of $Q(x,\alpha|x_{k},\alpha_{k})$ is unconstrained on $x$. At the maximum of $Q(x,\alpha|x_{k},\alpha_{k})$, we have

[TABLE]

where $0\leq j<L$. With $p(y_{m}|l_{1},l_{2},x_{k})$ given by (III-B), we obtain

[TABLE]

Substituting (B) and (C) into (C), we obtain

[TABLE]

Rearranging the terms, we obtain the update rule in (III-B).