Bispectrum Inversion with Application to Multireference Alignment

TL;DR

This paper introduces a novel approach for multireference alignment using invariant features like the bispectrum, enabling direct signal estimation from noisy data with provable recovery guarantees.

Contribution

It proposes new algorithms based on invariant features for signal recovery in MRA, including non-convex optimization, frequency marching, semidefinite relaxation, and integer programming, with theoretical guarantees.

Findings

Algorithms recover signals exactly without noise.

Methods are robust to high noise levels.

Invariant features enable stable estimation with fewer data passes.

Abstract

We consider the problem of estimating a signal from noisy circularly-translated versions of itself, called multireference alignment (MRA). One natural approach to MRA could be to estimate the shifts of the observations first, and infer the signal by aligning and averaging the data. In contrast, we consider a method based on estimating the signal directly, using features of the signal that are invariant under translations. Specifically, we estimate the power spectrum and the bispectrum of the signal from the observations. Under mild assumptions, these invariant features contain enough information to infer the signal. In particular, the bispectrum can be used to estimate the Fourier phases. To this end, we propose and analyze a few algorithms. Our main methods consist of non-convex optimization over the smooth manifold of phases. Empirically, in the absence of noise, these non-convex algorithms appear to converge to the target signal with random initialization. The algorithms are also robust to noise. We then suggest three additional methods. These methods are based on frequency marching, semidefinite relaxation and integer programming. The first two methods provably recover the phases exactly in the absence of noise. In the high noise level regime, the invariant features approach for MRA results in stable estimation if the number of measurements scales like the cube of the noise variance, which is the information-theoretic rate. Additionally, it requires only one pass over the data which is important at low signal-to-noise ratio when the number of observations must be large.

Tables

Table 1. Table I: Comparison of main MRA approaches. F ​ ( N ) 𝐹 𝑁 F(N) denotes the complexity of inverting the bispectrum. For instance, for the FM algorithm, F ​ ( N ) = O ​ ( N 2 ) 𝐹 𝑁 𝑂 superscript 𝑁 2 F(N)=O(N^{2}) . Storage requirements include the possibility of streaming computations where possible.

Method	Computational complexity	Storage requirement	Comments
Template alignment	$O\left(MN\log N\right)$	$O\left(N\right)$	Fails at moderate SNR (see Figure I.1)
Angular synchronization	$O\left(M^{2}N\log N\right)$	$O\left(M^2\right)$	Fails at low SNR
Expectation maximization	$O\left(TMN\log N\right)$	$O\left(MN\right)$	Empirically accurate; $\#$iterations $T$ grows with noise level
Invariant features (this paper)	$O\left(M{N}^2+F\left(N\right)\right)$	$O\left(N^2\right)$	Under mild conditions, accurate estimation if $M$ grows as ${\sigma }^6$

Equations

$${\xi }_j$$

$${\hat{r}}_j=\arg \underset{k}{\max }\Re \left\{\sum _{n=0}^{N-1}{\xi }_1[n]\overline{{\xi }_j[n+k]}\right\},j=2,\cdots ,M,$$

$$c_1$$

$$c_2[n_1]$$

$$c_3[n_1,n_2]$$

$${\hat{\mu }}_x=\frac{1}{M}\sum _{j=1}^M\left(\frac{1}{N}\sum _{n=0}^{N-1}{\xi }_j[n]\right)\sim \mathcal{N}\left({\mu }_x,\frac{{\sigma }^2}{NM}\right).$$

$$P_x[k]$$

$$\mathrm{DFT}\left(R_{s}x\right)[k]$$

$$\mathbb{E}\left\{P_{{\xi }_j}[k]\right\}$$

$${\hat{P}}_x[k]$$

$${\hat{\sigma }}^2$$

$$B_x[k_1,k_2]=y[k_1]\overline{y[k_2]}y[k_2-k_1].$$

$$B_{R_{s}x}[k_1,k_2]$$

$$\cdot \left(y[k_2-k_1]e^{2\pi i(k_1-k_2)s/N}\right)$$

$$=B_x[k_1,k_2].$$

$$B_x=y{y}^{\ast }\circ T(y),$$

$$\mathbb{E}\left\{B_{{\xi }_j}\right\}=B_x+{\sigma }^2{N}^2{\mu }_{x}A,$$

$$A_{\mathbb{R}}=\left[\begin{array}{cccccc}3& 1& 1& 1& \dots & 1\\ 1& 1& 0& 0& \dots & 0\\ 1& 0& 1& 0& \dots & 0\\ 1& 0& 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& 0& 0& 0& \dots & 1\end{array}\right],A_{\mathbb{C}}=\left[\begin{array}{cccccc}2& 1& 1& 1& \dots & 1\\ 0& 1& 0& 0& \dots & 0\\ 0& 0& 1& 0& \dots & 0\\ 0& 0& 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& 0& \dots & 1\end{array}\right].$$

$${\hat{B}}_{x-{\mu }_x}=\frac{1}{M}\sum _{j=1}^M{B}_{{\xi }_j-{\hat{\mu }}_x},$$

$${\tilde{B}}_x[k_1,k_2]=\tilde{y}[k_1]\overline{\tilde{y}[k_2]}\tilde{y}[k_2-k_1].$$

$${\hat{B}}_x=\frac{1}{M}\sum _{j=1}^M{B}_{{\xi }_j}-{\sigma }^2{N}^2{\hat{\mu }}_{x}A$$

$$:=\frac{\Vert B_{x-{\mu }_x}-\frac{1}{M}{\sum }_{j=1}^M{B}_{{\xi }_j-{\hat{\mu }}_x}{\Vert }_{\mathrm{F}}}{\Vert B_{x-{\mu }_x}{\Vert }_{\mathrm{F}}},$$

$$\underset{z\in {\mathbb{C}}^N}{\min }{\Vert W\circ \left(\tilde{B}-z{z}^{\ast }\circ T(z)\right)\Vert }_{\mathrm{F}}^2\text{ subject to }|z[k]|=1,\forall k.$$

$${\Vert W\circ \left(\tilde{B}-z{z}^{\ast }\circ T(z)\right)\Vert }_{\mathrm{F}}^2=\Vert W\circ \tilde{B}{\Vert }_{\mathrm{F}}^2$$

$$+{\Vert W\circ z{z}^{\ast }\circ T(z)\Vert }_{\mathrm{F}}^2-2⟨W\circ \tilde{B},W\circ z{z}^{\ast }\circ T(z)⟩,$$

$$⟨U,V⟩$$

$$⟨W\circ \tilde{B},W\circ z{z}^{\ast }\circ T(z)⟩=⟨z,M(z)z⟩,$$

$$M(z):=W^{(2)}\circ \tilde{B}\circ \overline{T(z)},$$

$$\underset{z\in {\mathbb{C}}^N}{\max }f(z)=⟨z,M(z)z⟩\text{ subject to }$$

$$\begin{array}{rl}\hat{y}=\arg \underset{z\in {\mathbb{C}}^N}{\max }\Re \left\{z^{\ast }M(z)z\right\}\text{mbox}subjectto& |z[k]|=1,\forall k,\\ \text{(if x is real) }& z[k]=\overline{z[-k]},\forall k,\end{array}$$

$$\tilde{B}\circ \overline{T({\hat{y}}_{k-1})}\approx {\hat{y}}_{k-1}{\hat{y}}_{k-1}^{\ast }.$$

$$\arg \underset{z\in {\mathbb{C}}^N}{\max }\Re \left\{z^{\ast }\left(\tilde{B}\circ \overline{T({\hat{y}}_{k-1})}\right)z\right\}\text{ subject to }|z[\ell ]|=1,\forall \ell ,$$

$$\mathrm{phase}\left(B_x[N-1,1]B_x[1,2]\cdot B_x[1,2]\cdots B_x[1,N-1]\right).$$

Full text

I Introduction

We consider the problem of estimating a discrete signal from multiple noisy and translated (i.e., circularly shifted) versions of itself, called multireference aligment (MRA). This problem occurs in a variety of applications in biology [1, 2, 3, 4], radar [5, 6], image registration and super-resolution [7, 8, 9], and has been the subject of recent theoretical analysis [10, 11]. The MRA model reads

[TABLE]

where $\varepsilon_{j}$ are i.i.d. normal random vectors with variance $\sigma^{2}$ and the underlying signal $x$ is in $\mathbb{R}^{N}$ or in $\mathbb{C}^{N}$. Operator $R_{r_{j}}$ rotates the signal ${x}$ circularly by $r_{j}$ locations, namely, $(R_{r_{j}}{x})[n]={x}[n-r_{j}]$, where indexing is zero-based and considered modulo $N$ (throughout the paper). While both $x$ and the translations $\{r_{j}\}$ are unknown, we stress that the goal here is merely to estimate $x$. This estimation is possible only up to an arbitrary translation.

A chief motivation for this work arises from the imaging technique called single particle Cryo-Electron Microscopy (Cryo-EM), which allows to visualize molecules at near-atomic resolution [12, 13]. In Cryo-EM, we aim to estimate a three dimensional (3D) object from its two-dimensional (2D) noisy projections, taken at unknown viewing directions [14, 15]. While typically the recovery process involves alignment of multiple observations in a low signal-to-noise ratio (SNR) regime, the underlying goal is merely to estimate the 3D object. In this manner, with the unknown shifts corresponding to the unknown viewing directions, MRA can be understood as a simplified model for Cryo-EM.

Existing approaches for MRA can be classified into two main categories. The first class of methods aims to estimate the set of translations $\{r_{j}\}$ first. Given this set, estimating $x$ can be achieved easily by aligning all observations $\xi_{j}$ and then averaging to reduce the noise. The second class, which we favor in this paper, consists of methods which aim to estimate the signal directly, without estimating the shifts.

Considering the first class, one intuitive approach to estimating the translations is to fix a template observation, say $\xi_{1}$, and to estimate the relative translations by cross-correlation. This is called template matching. Specifically, $r_{j}$ is estimated as

[TABLE]

where $\Re\{z\}$ and $\overline{z}$ denote the real part and the conjugate of a complex number $z$. This approach requires only one pass over the data: for each observation, the best shift can be computed in ${O}(N\log N)$, and the aligned observations can be averaged online. This results in a total computational cost of ${O}(MN\log N)$, see Table I. While this approach is simple and efficient, it necessarily fails below a critical SNR—see Figure I.1 for a representative example.

The issue with template matching is that we rely on aligning each observation to only one template: this is error prone at low SNR. Instead, to derive a more robust estimator, one can look for the most suitable alignment among all pairs of observations. The $M^{2}$ relative shifts thus computed must then be reconciled into a compatible choice of $M$ shifts for the individual observations. This is a discrete version of the angular synchronization problem, see [16, 17, 18, 19, 20, 21]. The computational complexity of aligning all pairs individually is ${O}(M^{2}N\log N)$, while storing the results uses ${O}(M^{2})$ memory.

Alternative algorithms for estimating the translations are based on different SDP relaxations [22, 23], iterative template alignment [24], zero phase representations [5] and neural networks [6]. The statistical limits of alignment tasks were derived for a variety of setups and noise models, see for instance [25, 26, 27, 28]. For example, for a continuous, 2D version of the MRA model, it was shown that the Cramér–Rao lower bound (CRLB) for translation estimation is proportional to the noise variance $\sigma^{2}$ [25]; crucially, it does not improve with $M$, even if the underlying signal is known. This is motivation to consider the second category of MRA methods, where shifts are not estimated.

Section VI elaborates on expectation maximization (EM) which tries to compute the maximum marginalized likelihood estimator (MMLE) of the signal—marginalization is done over the shifts. This method acknowledges the difficulty of alignment by working not with estimates of the shifts themselves, but rather with estimates of the probability distributions of the shifts. As a result, EM achieves excellent numerical performance in practice. However its computational complexity is high and its performance is not understood in theory.

It has been shown recently that the sample complexity of MRA, under assumption that shifts are distributed uniformly, is proportional to $\sigma^{6}$ in the low SNR regime. In other words, the number of measurements $M$ needs to scale like $\sigma^{6}$ to retain a constant estimation error [11].

In this work, we propose a framework which achieves this sample complexity by estimating the sought signal $x$ directly using features that are invariant under translations. For instance, the mean of $x$ is invariant under translation and can be estimated easily from the mean of all observations. We further use the power spectrum and the bispectrum of the observations—which are Fourier-transform based invariants—to estimate the magnitudes and phases of the signal’s Fourier transform, respectively.

For any fixed noise level (which may be arbitrarily large), these features can be estimated accurately provided sufficiently many measurements are available. Hence, our approach allows to deal with any noise level. Besides achieving the sample complexity, the computational complexity and memory requirements of the methods we describe are relatively low. Indeed, the only operations whose computational cost grows with $M$ are computations of averages over the data. These can be performed on-the-fly and are easily parallelizable. We mention that a recent tensor decomposition algorithm also achieves this estimation rate [29].

Given estimators for the mean and power spectrum of $x$, estimating the DC component and Fourier magnitudes of $x$ is straightforward. In this paper, we thus focus on the task of recovering the Fourier phases of $x$ from an estimator of its bispectrum. We propose two non-convex optimization algorithms on the manifold of phases for this task, which we call bispectrum inversion. We also discuss three additional algorithms which do not require initialization (and hence could be used to initialize others), based on frequency marching, SDP relaxation and integer programming. The first two methods recover the phases exactly in the absence of noise.

Beyond MRA, the bispectrum plays a central role in a variety of signal processing applications. For instance, it is a key tool to separate Gaussian and non-Gaussian processes [30, 31]. It is also used to investigate the cosmic background radiation [32, 33], seismic signal processing [34], image deblurring [35], feature extraction for radar [36], analysis of EEG signals [37], MIMO systems [38] and classification [39] (see also [40, 41, 42, 43, 44] and references therein). In Section III, we review previous works on bispectrum inversion [45, 46, 34]. Reliable algorithms to invert the bispectrum, as studied here, may prove useful in some of these applications.

The paper is organized as follows. Section II discusses the invariant feature approach for MRA. Section III presents the non-convex algorithms on the manifold of phases for bispectrum inversion. Section IV is devoted to additional algorithms that can be used to initialize the non-convex algorithms. Section V analyzes one of the proposed non-convex algorithm. Section VI elaborates on the EM approach for MRA, Section VII shows numerical experiments and Section VIII offers conclusions and perspective.

Throughout the paper we use the following notation. Vectors $x$ in $\mathbb{R}^{N}$ or $\mathbb{C}^{N}$ and $y\in\mathbb{C}^{N}$ denote the underlying signal and its discrete Fourier transform (DFT), respectively. In the sequel, all indices are understood modulo $N$, namely, in the range $0,\ldots,N-1$. The phase of a complex scalar $a$, defined as $a/|a|$ if $a\neq 0$ and zero otherwise, is denoted by $\mathrm{phase}(a)$ or $\tilde{a}$. The conjugate-transpose of a vector $z$ is denoted by $z^{*}$. We use ${}^{\prime}\circ^{\prime}$ to denote the Hadamard (entry-wise) product, $\mathbb{E}$ for expectation, Tr$(Z)$ for the trace and $\|Z\|_{\mathrm{F}}$ for the Frobenius norm of a matrix $Z$. We reserve $T(z)$ for circulant matrices determined by their first row $z$, i.e., $T(z)[k_{1},k_{2}]=z[k_{2}-k_{1}]$, and $\mathcal{H}^{N}$ for the set of Hermitian matrices of size $N\times N$.

II Multireference Alignment via Invariant Features

We propose to solve the MRA problem directly using features that are invariant under translations. Unlike pairwise alignment, this approach fuses information from all $M$ observations together—not just of pairs—and it only aims to recover the signal itself—not the translations. The essence of this idea was discussed as a possible extension in [22, Appendix A]. The invariant features can be understood either as auto-correlation functions or as their Fourier transform. In this work, we make use of the first three invariants defined as

[TABLE]

for $n_{1},n_{2}=0,\ldots,N-1$. It is clear that $c_{1},c_{2},c_{3}$ are invariant under circular shifts of $x$. For higher-order invariants based on auto-correlations, see for instance [47].

The first feature is the mean of the signal which is the auto-corrleation function of order one (i.e., $c_{1}$ in (II)). The distribution of the mean of $\xi_{j}$ is then given by $\mu_{\xi_{j}}\sim\mathcal{N}\left(\mu_{x},\frac{\sigma^{2}}{N}\right)$ and we can estimate $\mu_{x}$ as

[TABLE]

Estimating the signal’s mean supplies only limited information about the signal itself. Thus, we consider also the auto-correlation function of order two (i.e., $c_{2}$ in (II)). Its Fourier transform, the power spectrum, is explicitly defined as

[TABLE]

for all $k$, where $y$ is the DFT of $x$. An alternative way to understand the invariance of the power spectrum under shifts is through the effect of shifts on the DFT of a signal:

[TABLE]

Thus, shifts only affect the phases of the DFT, so that $P_{R_{s}x}=P_{x}$ for any shift $R_{s}$. Furthermore, owing to independence of the noise with respect to the signal itself and to the shift,

[TABLE]

where the second term is the power spectrum of the noise $\varepsilon_{j}$. Therefore, we estimate the power spectrum of $x$ as:

[TABLE]

It can be shown that ${\hat{P}_{x}}$ is unbiased and its variance is dominated by $\frac{\sigma^{4}}{M}$ for large $\sigma$. Hence, $\hat{P}_{x}\rightarrow P_{x}$ as $M\rightarrow\infty$. In particular, accurate estimation of the power spectrum requires $M$ to scale like $\sigma^{4}$. In the sequel, we assume that $\sigma$ is known.111If $\sigma$ is not known, it can be estimated from the data as

$\displaystyle\hat{\sigma}^{2}$ $\displaystyle=\frac{1}{N}.\operatorname{variance}\left(\sum_{n=0}^{N-1}\xi_{j}[n]\right)_{j=1,\ldots,M}.$

Recovering a signal from its power spectrum is commonly referred to as phase retrieval. This problem received considerable attention in recent years, see for instance [48, 49, 50, 51, 52, 53, 54]. It is well known that almost no one-dimensional signal can be determined uniquely from its power spectrum. Therefore, we use the power spectrum merely to estimate the signal’s Fourier magnitudes. As explained next, we use the auto-correlation of third order and its Fourier transform, the bispectrum, to estimate the Fourier phases.

Since phase retrieval is in general ill posed, we use the auto-correlation function of order three (that is, $c_{3}$ in (II)) through its Fourier transform, the bispectrum, to estimate the Fourier phases of the sought signal. The bispectrum is a function of two frequencies $k_{1},k_{2}=0,\dots,N-1$ and is defined as [55]:

[TABLE]

Note that, if $y[0]\neq 0$, the power spectrum is explicitly included in the bispectrum since $P_{x}[k]={B_{x}[k,k]}/{y[0]}$. The fact that the bispectrum is invariant under shifts can also be deduced from (II.3). Indeed, for any shift $R_{s}$,

[TABLE]

In matrix notation, we express this as

[TABLE]

where $T(y)$ is a circulant matrix whose first row is $y$, that is, $T(y)[k_{1},k_{2}]=y[k_{2}-k_{1}]$. Observe that if $x$ is real, then $y[k]=\overline{y[-k]}$ so that $T(y)$ and $B_{x}$ are Hermitian matrices. Simple expectation calculations lead to the conclusion that

[TABLE]

where $A=A_{\mathbb{R}}$ or $A=A_{\mathbb{C}}$ depending on $x\in\mathbb{R}^{N}$ or $x\in\mathbb{C}^{N}$ and

[TABLE]

Since the bias term is proportional to $\mu_{x}$, we propose to estimate ${\hat{B}_{x-\mu_{x}}}$ by averaging over $B_{\xi_{j}-\mu_{x}}$ for all $j$. This estimator is unbiased and its variance is controlled by $\frac{\sigma^{6}}{M}$ for large $\sigma$. Therefore, $M$ is required to scale like $\sigma^{6}$ to ensure accurate estimation. In practice, $\mu_{x}$ is not known exactly. Thus, we estimate the bispectrum by

[TABLE]

which is asymptotically unbiased. For finite $M$ and large $\sigma$, bias induced by the approximation $\hat{\mu}_{x}\approx\mu_{x}$ is significantly smaller than the standard deviation of (II.8).

The bispectrum contains information about the Fourier phases of $x$ because, defining $\tilde{y}[k]=\mathrm{phase}(y[k])$ and $\tilde{B}_{x}[k_{1},k_{2}]=\mathrm{phase}({B}_{x}[k_{1},k_{2}])$ where $\mathrm{phase}$ extracts the phase of a complex number (and returns 0 if that number is 0), we have

[TABLE]

In matrix notation, the normalized bispectrum takes the form $\tilde{B}_{x}=\tilde{y}\tilde{y}^{*}\circ T(\tilde{y})$.

Contrary to the power spectrum, the bispectrum is usually invertible. Indeed, in the absence of noise, the bispectrum determines the sought signal uniquely under moderate conditions:

Proposition II.1**.**

For $N\geq 5$, let $x\in\mathbb{C}^{N}$ be a signal whose DFT $y$ obeys $y[k]\neq 0$ for $k=1,\dots,K$, possibly also for $k=0$, and zero otherwise. Up to integer time shifts, $x$ is determined exactly by its bispectrum provided $K\geq\frac{N+1}{2}$.

For $N\geq 5$, let $x\in\mathbb{R}^{N}$ be a real signal whose DFT $y$ obeys $y[k]\neq 0$ for $k=1,\ldots,K$ and $k=N-1,\ldots,N-K$, possibly also for $k=0$, and zero otherwise. Up to integer time shifts, $x$ is determined exactly by its bispectrum provided $\frac{N}{3}\leq K\leq\frac{N-1}{2}$.

Proof.

This is a direct corollary of Lemmas V.1 and -B.1. ∎

We stress that the bispectrum estimator in (II.8) is not a bispectrum itself, since the set of bispectra is not a linear space: $\hat{B}_{x-\mu_{x}}$ is not invertible as such [42]. Algorithms we propose aim to find a stable inverse, in the sense that the recovered signal will have a bispectrum which is close to the estimated bispectrum in $\mathbb{C}^{N\times N}$. The following propositions combined argue formally that this can be done in the MRA model. The proofs in Appendix -A are constructive.

Proposition II.2** (Stable bispectrum inversion).**

There exists an estimator $\hat{x}$ with the following property. For any signal $x$ in $\mathbb{R}^{N}$ or $\mathbb{C}^{N}$ whose DFT is non-vanishing, there exist a precision $\delta=\delta(x)>0$ and a sensitivity $L=L(x)<\infty$ such that if an estimator $\hat{B}_{x}$ of $B_{x}$ satisfies $\|\hat{B}_{x}-B_{x}\|_{\mathrm{F}}\leq\delta$, then $\hat{x}=\hat{x}(\hat{B}_{x})$ satisfies $\min_{r=0\ldots N-1}\|x-R_{r}\hat{x}\|_{2}\leq L\|\hat{B}_{x}-B_{x}\|_{\mathrm{F}}$.

Proposition II.3** (Bispectrum estimation).**

For any signal $x$ in $\mathbb{R}^{N}$ or $\mathbb{C}^{N}$ whose DFT is non-vanishing, for any required precision $\delta>0$ and for any probability $p<1$, there exists a constant $C=C(x,p,\delta)<\infty$ such that, for any noise level $\sigma>0$, if the number of observations $M$ exceeds $C\cdot(\sigma^{2}+\sigma^{6})$, the estimator

[TABLE]

satisfies $\|\hat{B}_{x}-B_{x}\|_{\mathrm{F}}\leq\delta$ with probability at least $p$.

We mention that uniqueness in the continuous setup was considered in [56]. The more general setting of bispectrum over compact groups was considered in [57, 58, 59, 60].

The MRA model here assumes i.i.d. Gaussian noise. However, the estimation is performed by averaging in the bispectrum domain, where noise affecting individual entries is correlated. Consequently, one may want to use a more robust estimator, such as the median. Yet, computing the median of complex matrices is computationally expensive, while computing the average can be performed efficiently and on-the-fly, that is, without requiring to store all observations. For Gaussian noise, we have noticed numerically that using the mean or the median for bispectrum estimation leads to comparable estimation errors (experiments not shown). In other noise models, e.g., with outliers, it might be useful to consider the median or the median of means method, see for instance [61].

Figure II.1 presents the relative estimation error of the power spectrum and bispectrum as a function of the number of observations $M$. For the bispectrum, the relative error is computed as

[TABLE]

and similarly for the power spectrum. As expected, the slope of all curves is approximately $1/2$ in logarithmic scale, implying that the estimation error decreases as $O(1/\sqrt{M})$. The invariant features approach for MRA is summarized in Algorithm 1.

Consider the case in which the number of samples $M$ may be very large whereas the size of the object is fixed, namely $N\ll M$. This case is of interest in many applications, such as cryo-EM [14, 15]. In this regime, the invariant features approach has two important advantages over methods that rely on estimating the translations. First, in the invariant features approach, we average over the $M$ observations (which is computationally cheap), and then apply a more complex algorithm (say, to recover a signal from its bispectrum) whose input size is a function of $N$ but is independent of $M$. Hence, the overall complexity of this approach can be relatively low. Second, the alignment-based method requires storing all $M$ observations, namely, $MN$ samples, which is unnecessary in the invariant features approach. There, for each observation, we just need to compute its invariants, to be averaged over all observations: this can be done online (in streaming mode) and in parallel.

III Non-convex Algorithms for bispectrum inversion

After estimating the first Fourier coefficient $y[0]$ as $N\hat{\mu}_{x}$, our approach for MRA by invariant features consists of two parts. We use the power spectrum to estimate the signal’s Fourier magnitudes and the bispectrum for the phases. The first part is straightforward: $|y[k]|$ can be estimated as $\sqrt{\hat{P}_{x}[k]}$ if $\hat{P}_{x}[k]\geq 0$, and as 0 otherwise. Hereafter, we focus on estimation of the phases of the DFT, $\tilde{y}$.

In the literature, two main approaches were suggested to invert the discrete bispectrum. The first is based on estimating the frequencies one after the other by exploiting simple algebraic relations [45, 46]. The second approach suggests to estimate the signal by least-squares solution and phase unwrapping [46, 34]. We improve these methods and suggest a few new algorithms. The algorithms are split into two sections. This section is devoted to two new non-convex algorithms based on optimization on the manifold of phases. Both of these algorithms require initialization. While experimentally it appears that random initialization works well, for completeness, in the next section we propose three additional algorithms which do not need initialization and hence could be used to initialize the non-convex algorithms.

III-A Local non-convex algorithm over the manifold of phases

In this section, similarly to (II.9), we let $\tilde{B}$ denote our estimate of the phases of $B_{x}$. Since $\tilde{B}\approx\tilde{y}\tilde{y}^{*}\circ T(\tilde{y})$, one way to model recovery of the Fourier phases $\tilde{y}$ is by means of the non-convex least-squares optimization problem

[TABLE]

The matrix $W\in\mathbb{R}^{N\times N}$ is a weight matrix with nonnegative entries. These weights can be used to indicate our confidence in each entry of $\tilde{B}$. Expanding the squared Frobenius norm yields

[TABLE]

where

[TABLE]

is the real inner product associated to the Frobenius norm. Under the constraints on $z$, the first two terms are constant and the inner product term is equivalent to

[TABLE]

with

[TABLE]

where we use the notation $W^{(2)}:=W\circ W$. One possibility is to choose $W=\sqrt{|\hat{B}_{x-\mu_{x}}|}$, where the absolute value and the square root are taken entry-wise, so that $M(z)=\hat{B}_{x-\mu_{x}}\circ\overline{T(z)}$. Hence, optimization problem (III.1) is equivalent to

[TABLE]

We can also impose $z[0]=\mathrm{phase}(\hat{\mu}_{x})$. If $x$ is real, we have the additional symmetry constraints $z[k]=\overline{z[-k]}$.

Since the cost function $f$ is continuous and the search space is compact, a solution exists. Of course, the solution is not unique, in accordance with the invariance of the bispectrum under integer time-shifts of the underlying discrete signal. This is apparent through the fact that the cost function $f$ is invariant under the corresponding (discrete) transformations of $z$. This is true independently of the data $\tilde{B}$ and $W$. The proof is in Appendix -C.

Lemma III.1**.**

The cost function $f$ is invariant under transformations of $z$ that correspond to integer time-shifts of the underlying signal.

To solve this non-convex program, we use the Riemannian trust-region method (RTR) [62], whose usage is simplified by the toolbox Manopt [63]. RTR enjoys global convergence to second-order critical points, that is, points which satisfy first- and second-order necessary optimality conditions [64] and a quadratic local convergence rate. Empirically, in the noiseless case it appears that the algorithm recovers the target signal with random initialization, all local minima are global (with minor technicality for even $N$ in the real case) and all second-order critical points have an escape direction, that is, saddles are “strict”. Numerical experiments demonstrate reasonable robustness in the face of noise. This algorithm is summarized in Algorithm 2 and studied in detail in Section V.

III-B Iterative phase synchronization algorithm

In this section we present an alternative heuristic to the non-convex algorithm on the manifold of phases. This algorithm is based on iteratively solving the phase synchronization problem. Suppose we get an estimation of $\tilde{y}$, say $\hat{y}_{k-1}$. If $\hat{y}_{k-1}\approx\tilde{y}$ is non-vanishing, then this estimation should approximately satisfy the bispectrum relation:

[TABLE]

The underlying idea is now to push the current estimation towards $\tilde{y}$ by finding a rank-one approximation of $\tilde{B}\circ\overline{T(\hat{y}_{k-1})}$ with unit modulus entries. This problem can be formulated as:

[TABLE]

where we treat the matrix $\tilde{B}\circ\overline{T(\hat{y}_{k-1})}$ as a constant. This problem is called phase synchronization. Many algorithms have been suggested to solve the phase synchronization problem. Among them are the eigenvector method, SDP relaxation, projected power method, Riemannian optimization and approximate message passing [16, 20, 17, 18, 19]. Notice that the solution of (III.5) is only defined up to a global phase, namely, if $z$ is optimal, then so is $ze^{i\phi}$ for any angle $\phi$. To resolve this ambiguity, we require knowledge of the phase of the mean, $\tilde{y}[0]$ (which is easy to estimate from the data) and we pick the global phase of $\hat{y}_{k}$ such that $\hat{y}_{k}[0]=\tilde{y}[0]$.

The $k$th iteration of our algorithm thus (tries to) solve the phase synchronization problem with respect to the matrix $M_{k-1}:=\tilde{B}\circ\overline{T(\hat{y}_{k-1})}$, where $\hat{y}_{k-1}$ is the solution of the previous estimation. Assuming the signal is real, we also impose at each iteration the conjugate-reflection property of $\hat{y}_{k}[\ell]=\overline{\hat{y}_{k}[-\ell]}$ for all $\ell$ so that $M_{k}$ is Hermitian. In the numerical experiments in Section VII, we solve (III.5) by the Riemannian trust-region method described in [17]. Empirically, the performance of this algorithm and Algorithm 2 is indistinguishable. The algorithm is summarized in Algorithm 3.

IV Initialization-free Algorithms

The previous section was devoted to non-convex algorithms to invert the bispectrum. In this section we present three additional algorithms based on frequency marching (FM), SDP relaxation and phase unwrapping. These algorithms do not require initialization and therefore could be used to initialize the non-convex algorithms.

We prove that FM and the SDP recover the Fourier phases exactly in the absence of noise under the assumption that we can fix $\tilde{y}[1]$. If the signal has non-vanishing DFT, $\tilde{y}[1]$ can be estimated from the bispectrum using the fact that $\tilde{y}[1]^{N}$ equals

[TABLE]

(Any $N$th root can be used for $\tilde{y}[1]$, corresponding to the $N$ possible shifts of $x$.) In all cases, we argue that forcing $\tilde{y}[1]=1$ is acceptable if $N$ is large. Indeed, recall that a shift by $\ell$ entries in the space domain is equivalent to modulating the $k$th Fourier coefficient by $e^{-2\pi i\ell k/N}$. In particular, it means that the phase $\tilde{y}[1]$ can be shifted by $e^{-2\pi i\ell/N}$ for an arbitrary $\ell\in\mathbb{Z}$. Thus, for signals of length $N\gg 1$, the phase $\tilde{y}[1]$ can be set arbitrarily with only small error. In the numerical experiments of Section VII, we give the correct value of $\tilde{y}[1]$ to the algorithms in order to assess their best possible behavior.

We begin by discussing the FM algorithm, which is a simple propagation method: it is exact in the absence of noise. Notwithstanding, its estimation for the low-frequency coefficients is sensitive to noise. Because of its recursive nature, error in the low frequencies propagates to the high frequencies, resulting in unreliable estimation. The other two algorithms are more computationally demanding but appear more robust.

IV-A Frequency marching algorithm

The FM algorithm is a simple propagation algorithm in the spirit of [45, 46] that aims to estimate $\tilde{y}$ one frequency at a time. This algorithm has computational complexity $O(N^{2})$ and it recovers $\tilde{y}$ exactly for both real and complex signals in the absence of noise, assuming $\tilde{y}[1]$ is known.

Let us denote $\tilde{B}[k_{1},k_{2}]=e^{i\Psi[k_{1},k_{2}]}$ and $\tilde{y}[k]=e^{i\psi[k]}$. Accordingly, we can reformulate (II.5) as

[TABLE]

where the modulo is taken over the sum of all three terms. Using this relation, we can start to estimate the missing phases. The first unknown phase, $\psi[2]$, can be estimated by:

[TABLE]

where $\hat{\psi}[2]$ refers to the estimator of $\psi[2]$ (defined modulo $2\pi$). We can estimate the next phase in the same manner:

[TABLE]

For higher frequencies, we have more measurements to rely on. For the fourth entry, we now can derive two estimators as follows:

[TABLE]

and

[TABLE]

In the noiseless case, it is clear that ${\psi}[4]=\hat{\psi}^{(1)}[4]=\hat{\psi}^{(2)}[4]$. In a noisy environment, we can reduce the noise by averaging the two estimators, where averaging is done over the set of phases (namely, over the rotation group SO(2)) as explained in Appendix -D. Specifically,

[TABLE]

We can iterate this procedure. To estimate phase $q$, we want to consider all entries of $\tilde{B}[k,\ell]=\tilde{y}[k]\overline{\tilde{y}[\ell]}\tilde{y}[\ell-k]$ such that exactly one of the indices $k,\ell,$ or $\ell-k$ is equal to $q$ and all other indices are in $1,\ldots,q-1$, so that all other phases involved have already been estimated. A simple verification shows that only entries $\tilde{B}[p,q],\thinspace p=1,\dots,q-1,$ have that property. Furthermore, because of symmetry in the bispectrum (V.4), half of these entries are redundant so that only entries $\tilde{B}[p,q],\thinspace p=1,\dots,\left\lfloor\frac{q}{2}\right\rfloor$ remain. As a result, estimation of the $k$th phase relies on averaging over $\left\lfloor\frac{k}{2}\right\rfloor$ equations, as summarized in Algorithm 4, with the following simple guarantee. The above construction yields the following proposition.

Proposition IV.1**.**

Let $\tilde{B}=\tilde{B}_{x}$ be the normalized bispectrum as defined in (II.9) and assume that $\tilde{y}[1]$ is known. If $y[k]\neq 0$ for $k=1,\dots,K$, then Algorithm 4 recovers the Fourier phases $\tilde{y}[k],\thinspace k=1,\dots,K$ exactly.

We note in closing that, if the signal $x$ is real, symmetries in the phases $\tilde{y}$ and $\tilde{B}$ can be exploited easily in FM.

IV-B Semidefinite programming relaxation

In this section we assume that the DFT $y$ is non-vanishing so that the bispectrum relation can be manipulated as

[TABLE]

where $\overline{T(\tilde{y})}$ is its entry-wise conjugate. The developments are easily adapted if the signal has zero mean. Similarly to the FM algorithm, we assume that $\tilde{y}[0]$ and $\tilde{y}[1]$ are available. We aim to estimate $\tilde{y}$ by a convex program. As a first step, we decouple the bispectrum equation and write the problem of estimating $\tilde{y}$ as the following non-convex optimization problem:

[TABLE]

where $\mathcal{H}^{N}$ is the set of Hermitian matrices of size $N$ and $W\in\mathbb{R}^{N\times N}$ is a real weight matrix with positive entries. In particular, in the numerical experiments we set $W=|B|$.

In the absence of noise, the minimizers of (IV.1) satisfy the bispectrum equation. However, in general these cannot be computed in polynomial time. In order to make the problem tractable, we relax the non-convex coupling constraint $Z=zz^{*}$ to the convex constraint $Z\succeq zz^{*}$ (that is, $Z-zz^{*}$ is positive semidefinite). The convex relaxation is then given by

[TABLE]

Upon solving (IV.2), which can be done in polynomial time with interior point methods, the phases $\tilde{y}$ are estimated from $\mathrm{phase}(z)$. In practice, we use CVX to solve this problem [65]. The algorithm is summarized in Algorithm 5. We note that problem (IV.2) is not a standard SDP, in that its cost function is nonlinear.

In the noiseless case, the SDP relaxation (IV.2) recovers the missing phases exactly. Interestingly, the proof is not so much based on optimality conditions as it is on an algebraic property of circulant matrices. The proof of the following property is given in Appendix -E.

Lemma IV.2**.**

Let $\hat{u}$ be the DFT of a vector $u\in\mathbb{C}^{N}$ obeying $u[k]=\overline{u[-k]}$, so that $\hat{u}$ is real. If $u[0]=u[1]=1$ and $\hat{u}$ is non-negative, then $u[k]=1$ for all $k$.

The following theorem is a direct corollary of Lemma IV.2. The main proof idea is as follows. Consider $u=\overline{\tilde{y}}\circ z$ where $(Z,z)$ is optimal for the SDP; then, the constraints ensure $u[0]=u[1]=1$. Furthermore, one can see via the Schur complement that the constraints force $T(u)$ to be positive semidefinite. Since the eigenvalues of $T(u)$ are the DFT of $u$, it follows that $\hat{u}$ is non-negative, so that the lemma above applies and $u\equiv 1$, or, equivalently, $z=\tilde{y}$. Details of the proof are in Appendix -F.

Theorem IV.3**.**

For a real signal with non-vanishing DFT $y$, if all weights in $W$ are positive, $\tilde{y}[0]$ and $\tilde{y}[1]$ are known and the objective value of (IV.2) attains 0 (which is the case in the absence of noise), then the SDP has a unique solution given by $z=\tilde{y}$ and $Z=zz^{*}$.

We close with an important remark about the symmetry breaking purpose of constraint $z[1]=\tilde{y}[1]$ in the SDP. Because the signal $x$ can be recovered only up to integer time shifts, even in the noiseless case, without this constraint there are at least $N$ distinct solutions $(z,Z)$ to the SDP. Because SDP is a convex program, any point in the convex hull of these $N$ points is also a solution. Thus, if the symmetry is not broken, the set of solutions contains many irrelevant points. Furthermore, interior point methods tend to converge to a center of the set of solutions, which in this case is never one of the desired solutions.

IV-C Phase unwrapping by integer programming algorithm

The next algorithm is based on solving an over-determined system of equations involving integers. Let us denote $\tilde{y}[k]=e^{i\psi[k]}$ and $\tilde{B}[k_{1},k_{2}]=e^{i\Psi[k_{1},k_{2}]}$ so the normalized bispectrum model is given by

[TABLE]

By taking the logarithm, we get the algebraic relation

[TABLE]

where, as a result of phase wrapping, $\chi$ takes on integer values. Let $\Psi_{\textrm{vec}}$ and $\chi_{\textrm{vec}}$ be the column-stacked versions of $\Psi$ and $\chi$, respectively. Then, the model reads

[TABLE]

where the sparse matrix $A\in\mathbb{R}^{N^{2}\times N}$ encodes the right hand side of (IV.3). It can be verified that $A$ is of rank $N-1$ (see for instance [66]), with null space corresponding to the time-shift-induced ambiguity on the phases (II.3). Note that both the integer vector $\chi_{\textrm{vec}}$ and the phases $\psi$ are unknown. Given $\chi_{\textrm{vec}}$, the phases $\psi$ can be obtained easily by solving

[TABLE]

for some $\ell_{p}$ norm. Observe that any error in estimating $\chi$ may cause a big estimation error of $\psi$ in (IV.5). These errors can be thought of as outliers. Hence, we choose to use least unsquared deviations (LUD), $p=1$, which is more robust to outliers. The more challenging task is to estimate the integer vector $\chi_{\textrm{vec}}\in\mathbb{Z}^{N^{2}}$. To this end, we first eliminate $\psi$ from (IV.4) as follows. Let $C\in\mathbb{R}^{(N^{2}-(N-1))\times N^{2}}$ be a full rank matrix such that $CA=0$, that is, the columns of $C^{T}$ are in the null space of $A^{T}$. Matrix $C$ can be designed by at least two methods. One, suggested in [67], exploits the special structure of $A$ to design a sparse matrix composed of integer values. Another, which we use here, is to take $C$ to have orthonormal rows which form a basis of the kernel of $A^{T}$. Numerical experiments (not shown) indicate that the latter approach is more stable. Next, we multiply both sides of (IV.4) from the left by $C$ to get

[TABLE]

Therefore, the integer recovery problem can be formulated as

[TABLE]

where we minimize over all integers. Note that $C\Psi_{\textrm{vec}}$ is a known vector. The problem is then equivalent to finding a lattice vector with the basis $C$ which is as close as possible to the vector $-C\Psi_{\textrm{vec}}/(2\pi)$. While the problem is known to be NP-hard, we approximate the solution of (IV.6) with the LLL (Lenstra–Lenstra–Lovasz) algorithm, which can be run in polynomial time [68]. The LLL algorithm computes a lattice basis, called a reduced basis, which is approximately orthogonal. It uses the Gram–Schmidt process to determine the quality of the basis. For more details, see [69, Ch. 17].

We note that (IV.6) is under-determined as the matrix $C$ is of rank $N^{2}-\mbox{rank}(A)=N^{2}-(N-1)$. While the LLL algorithm works with under-determined systems, in our case we can solve it for a determined system since we can fix the first $N-1$ entries of $\chi_{\textrm{vec}}$ to be zero.222We omit the proof of this property here and only mention that it is based on the derivation in [67]. Once we have estimated $\chi_{\textrm{vec}}$, we solve (IV.5) with $p=1$. This approach is summarized in Algorithm 6.

V Analysis of optimization over phases

In this section, we study the non-convex optimization problem (III.4) and give more implementation details to solve it, since numerical experiments identify this as the method of choice for MRA from invariant features among all methods compared. We start by considering the general case of a complex signal $x\in\mathbb{C}^{N}$ and consider the real case in Appendix -B. Recall that we aim to maximize

[TABLE]

where the inner product is defined by (III.2), $W$ is a real weighting matrix and $W^{(2)}:=W\circ W$. The optimization problem lives on a manifold, that is, a smooth nonlinear space. Indeed, the smooth cost function $f(z)$ is to be maximized over the set

[TABLE]

which is a Cartesian product of $N$ unit circles in the complex plane (a torus). Theory and algorithms for optimization on manifolds can be found in the monograph [71]. We follow this formalism here. Details can also be found in [17], which deals with the similar problem of phase synchronization, using similar techniques. For the numerical experiments below, we use the toolbox Manopt which provides implementations of various optimization algorithms on manifolds [63].

Under mild conditions, the global optima of (III.4) correspond exactly to $\tilde{y}$ up to integer time shifts. This fact is proven in Appendix -G.

Lemma V.1**.**

For $N\geq 3$, let $x\in\mathbb{C}^{N}$ be a signal whose DFT $y$ is nonzero for frequencies $k$ in $\{1,\ldots,K\}$, possibly also for $k=0$, and zero otherwise. Up to integer time shifts, $x$ is determined exactly by its bispectrum $B$ provided $K\geq\frac{N+1}{2}$. Furthermore, the global optima of (III.4) correspond exactly to the relevant phases of $y$—up to the effects of integer time shifts—provided $W[k,\ell]$ is positive when $B[k,\ell]\neq 0$.

The problem at hand is

[TABLE]

This is smooth but non-convex, so that in general it is hard to compute the global optimum. We derive first- and second-order necessary optimality conditions. Points which satisfy these conditions are called critical and second-order critical points, respectively. Known algorithms converge to critical points (e.g., Riemannian gradient descent) and even to second-order critical points (e.g., Riemannian trust-regions) regardless of initialization [71, 62, 64]. Empirically, despite non-convexity, the global optimum appears to be computable reliably in favorable noise regimes.

As we proceed to consider optimization algorithms for (III.4), the gradient of $f$ will come into play:

[TABLE]

where $M^{\mathrm{adj}}\colon\mathbb{C}^{N\times N}\to\mathbb{C}^{N}$ is the adjoint of $M$ with respect to the inner product $\left\langle{\cdot},{\cdot}\right\rangle$. Formally, the adjoint is defined such that, for any $z\in\mathbb{C}^{N},X\in\mathbb{C}^{N\times N}$,

[TABLE]

Specifically, in Appendix -H we show that

[TABLE]

where $T_{k}$ is a circulant matrix with ones in its $k$th (circular) diagonal and zero otherwise, namely,

[TABLE]

As it turns out, under the symmetries of the problem at hand, there is no need to evaluate $M^{\mathrm{adj}}$ explicitly. Indeed, $\tilde{B}$ obeys

[TABLE]

This property is preserved when $\tilde{B}$ is obtained by averaging bispectra of multiple observations, as in (II.8). Assuming the same symmetry for the real weights $W$, we find below that $M^{\mathrm{adj}}(zz^{*})=M(z)z$. See Appendix -I.

Lemma V.2**.**

If $\tilde{B}[k_{2}-k_{1},k_{2}]=\tilde{B}[k_{1},k_{2}]$ and $W[k_{2}-k_{1},k_{2}]=W[k_{1},k_{2}]$ for all $k_{1},k_{2}$, then $M^{\mathrm{adj}}(zz^{*})=M(z)z$ for all $z\in\mathbb{C}^{N}$.

Thus, under the symmetries assumed in Lemma V.2, the gradient of $f$ simplifies and we get a simple expression for the Hessian as well:

[TABLE]

For unconstrained optimization, the first-order necessary optimality conditions are $\nabla f(z)=0$. In the presence of the constraint $z\in\mathcal{M}$, the conditions are different. Namely, following [71, eq. (3.37)], since $\mathcal{M}$ is a submanifold of $\mathbb{C}^{N}$, first-order necessary optimality conditions state that the orthogonal projection of the gradient $\nabla f(z)$ to the tangent space to $\mathcal{M}$ at $z$ must vanish. The result of this projection is called the Riemannian gradient. Formally, the tangent space is obtained by linearizing (differentiating) the constraints $|z[k]|^{2}=\left\langle{z[k]},{z[k]}\right\rangle=1$ for all $k$, yielding

[TABLE]

Orthogonal projection of $u\in\mathbb{C}^{N}$ to the tangent space $\mathrm{T}_{z}\mathcal{M}$ can be computed entry-wise by subtracting from each $u[k]$ its component aligned with $z[k]$. Let $\operatorname{Proj}_{z}\colon\mathbb{C}^{N}\to\mathrm{T}_{z}\mathcal{M}$ denote this projection. This operation admits a compact matrix notation as

[TABLE]

where $\operatorname{ddiag}\colon\mathbb{C}^{N\times N}\to\mathbb{C}^{N\times N}$ sets all non-diagonal entries of a matrix to zero. Equipped with this notion and the expression for $\nabla f(z)$ (V.5), it follows that the Riemannian gradient of $f$ at $z$ on $\mathcal{M}$ is

[TABLE]

with

[TABLE]

Lemma V.3**.**

If $z\in\mathcal{M}$ is optimal for (V.1), then $\operatorname{grad}\!f(z)=0$; equivalently, $\mathrm{diag}(\nabla f(z)z^{*})=\nabla f(z)\circ\overline{z}$ is real.

Proof.

See [72, Rem. 4.2 and Cor. 4.2]. For the equivalence, notice that $\operatorname{Proj}_{z}(u)=0$ if and only if $u[k]=\Re\{u[k]\overline{z[k]}\}z[k]$ for all $k$, and multiply by $\overline{z[k]}$ on both sides using $|z[k]|=1$. ∎

A point $z$ which satisfies these conditions is called a critical point. Likewise, we can define a notion of Riemannian Hessian as the linear, self-adjoint operator on $\mathrm{T}_{z}\mathcal{M}$ which captures infinitesimal changes in the Riemannian gradient around $z$. Without getting into technical details, we follow [71, eq. (5.15)] and define (with $\operatorname{D}\!$ the directional derivative operator):

[TABLE]

where $\operatorname{D}\!D(z)[\dot{z}]$ is a real, diagonal matrix. Its contribution to the Hessian is zero, since $(\operatorname{D}\!D(z)[\dot{z}])z$ vanishes under the projection $\operatorname{Proj}_{z}$. Hence,

[TABLE]

The Riemannian Hessian intervenes in the second-order necessary optimality conditions as follows.

Lemma V.4**.**

If $z\in\mathcal{M}$ is optimal for (V.1), then $\operatorname{grad}\!f(z)=0$ and $\operatorname{Hess}\!f(z)\preceq 0$, that is, for all $\dot{z}\in\mathrm{T}_{z}\mathcal{M}$ we have

[TABLE]

Proof.

See [72, Rem. 4.2 and Cor. 4.2]. In the equality, we used the fact that $\operatorname{Proj}_{z}$ is self-adjoint and $\dot{z}\in\mathrm{T}_{z}\mathcal{M}$. ∎

A point $z$ which satisfies these conditions is called a second-order critical point. With unit weights, the following lemma shows that second-order critical points $z$, in the noiseless case, cannot have an arbitrarily bad objective value $f(z)$. This result is weak, however, since empirically it is observed that in the noiseless case local optimization methods consistently converge to global optima whose value are $N^{2}$, suggesting that all second-order critical points are global optima in this simplified scenario. While we do not have a proof for this stronger conjecture, we provide the lemma below because it is analogous to [17, Lemma 14] which, in that reference, is a key step toward proving global optimality of second-order critical points.

Lemma V.5**.**

In the absence of noise and with unit weights, a second-order critical point $z$ of (V.1) satisfies

[TABLE]

In particular, this implies

[TABLE]

Proof.

See Appendix -J for the proof of the inequality. It follows from two key considerations. First, because $z$ is a critical point, Lemma V.3 indicates that $\nabla f(z)\circ\overline{z}$ is real. Second, because $z$ is second-order critical, the Riemannian Hessian at $z$ must be negative semidefinite by Lemma V.4. Applied to all tangent directions at $z$ which perturb only one phase at a time implies the desired inequality. The fact that $3f(z)=\left\langle{z},{\nabla f(z)}\right\rangle$ follows from (V.5). ∎

One final ingredient that is necessary to optimize $f$ over $\mathcal{M}$ is a means of moving away from a current iterate $z\in\mathcal{M}$ to the next by following a tangent vector $\dot{z}$. A simple means of achieving this is through a retraction [71, Def. 4.1.1]. For $\mathcal{M}$, an obvious retraction is the following:

[TABLE]

With the formalism of (V.1) and the above derivations, we can now run a local Riemannian optimization algorithm. As an example, the Riemannian gradient ascent algorithm would iterate the following:

[TABLE]

where $\eta^{(t)}>0$ is an appropriately chosen step size and $z^{(0)}\in\mathcal{M}$ is an initial guess. It is relatively easy to choose the step sizes such that the sequence $z^{(t)}$ converges to critical points regardless of $z^{(0)}$, with a linear local convergence rate [71, §4]. In practice, we prefer to use the Riemannian trust-region method (RTR) [62], whose usage is simplified by the toolbox Manopt [63]. RTR enjoys global convergence to second-order critical points [64] and a quadratic local convergence rate.

In this section, the analysis focused on complex signals. For real signals, we can follow the same methodology while taking the symmetry in the Fourier domain into account. This analysis is given in Appendix -B.

VI Expectation maximization

In this section, we detail the expectation maximization algorithm (EM) [73] applied to MRA. As the numerical experiments in Section VII demonstrate, EM achieves excellent accuracy in estimating the signal. However, compared to the invariant features approach proposed in this paper, it is significantly slower and requires many passes over the data (thus excluding online processing).

Let $X=\left[\xi_{1},\ldots,\xi_{M}\right]$ be the data matrix of size $N\times M$, following the MRA model (I.1). The maximum marginalized likelihood estimator (MMLE) for the signal $x$ given $X$ is the maximizer of the likelihood function $L(x;X)=p(X|x)$ (the probability density of $X$ given $x$). This density could in principle be evaluated by marginalizing the joint distribution $p(X,r|x)$ over the unknown shifts $r\in\{0,\ldots,N-1\}^{M}$. This, however, is intractable as it involves summing over $N^{M}$ terms.

Alternatively, EM tries to estimate the MMLE as follows. Given a current estimate for the signal $x_{k}$, consider the expected value of the log-likelihood function, with respect to the conditional distribution of $r$ given $X$ and $x_{k}$:

[TABLE]

This step is called the E-step. Then, iterate by computing the M-step:

[TABLE]

For the MRA model, this can be done in closed form. Indeed, the log-likelihood function follows from the i.i.d. Gaussian noise model:

[TABLE]

To take the expectation with respect to $r$, we need to compute $w_{k}^{\ell,j}$: for each observation $j$, this is the probability that the shift $r_{j}$ is equal to $\ell$, given $X$ and assuming $x=x_{k}$. This also follows easily from the i.i.d. Gaussian noise model:

[TABLE]

(with appropriate scale so that $\sum_{\ell=0}^{N-1}w_{k}^{\ell,j}=1$). This allows to write $Q$ down explicitly:

[TABLE]

This is a convex quadratic expression in $x$ with maximizer

[TABLE]

In words: given an estimator $x_{k}$, the next estimator is obtained by averaging all shifted versions of all observations, weighted by the empirical probabilities of the shifts. Considering all shifts of all observations would, in principle, induce an iteration complexity of $O(MN^{2})$, but fortunately, for each observation, the matrix of its shifted versions is circulant, which makes it possible to use FFT to reduce the overall computational cost to $O(MN\log N)$. See the available code for details. We note that Matlab naturally parallelizes the computations over $M$.

In practice, we set $x_{0}\sim\mathcal{N}(0,I_{N})$ to be a random guess. Furthermore, for $M\geq 3000$, we first execute 3000 batch iterations, where the EM update is computed based on a random sample of 1000 observations (fresh sample at each iteration). This inexpensively transforms the random initialization into a ballpark estimate of the signal. The algorithm then proceeds with full-data iterations until the relative change between two consecutive estimates drops below $10^{-5}$ (in $\ell_{2}$-norm, up to shifts).

VII Numerical Experiments

This section is devoted to numerical experiments, examining all proposed algorithms. Code for all algorithms and to reproduce the experiments is available online.333https://github.com/NicolasBoumal/MRA The experiments were conducted as follows. The true signal $x$ of length $N=41$ is a fixed window of height 1 and width 21. With this signal, the signal-to-noise ratio is $\frac{\|x\|^{2}}{\|\varepsilon\|^{2}}\approx\frac{1}{2\sigma^{2}}$. We generated a set of $M$ shifted noisy versions of $x$ as

[TABLE]

where each shift was randomly drawn from a uniform distribution over $\{0,\dots,N-1\}$ and $\varepsilon_{j}\sim\mathcal{N}(0,\sigma^{2}I)$ for all $j$. The relative recovery error for a single experiment is defined as

[TABLE]

where $\hat{x}$ is the estimation of the signal. All results are averaged over 20 repetitions. While we present here results for a specific signal, alternative signal models (e.g., random signals) showed similar numerical behavior.

The following figures compare the recovery errors for all proposed algorithms, with random initialization for those that need initialization. The non-convex algorithm on the manifold of phases, Algorithm 2, runs the Riemannian trust-region method (RTR) [62] using the toolbox Manopt [63]. Algorithm 3 runs 15 iterations with warm-start using the same toolbox. For the phase unwrapping algorithm, Algorithm 6, we use an implementation of LLL available in the MILES package [70]. The SDP is solved with CVX [65]. The EM algorithm is implemented as explained in Section VI. We compared the algorithms with an oracle who knowns the random shifts $r_{j}$ and therefore simply averages out the Gaussian noise. Experiments are run on a computer with 30 CPUs available. These CPUs are used to compute the invariants in parallel (with Matlab’s parfor), while the EM algorithm benefits from parallelism to run the many thousands of FFTs it requires efficiently (built-in Matlab). The algorithms that need $\tilde{y}[0]$ and $\tilde{y}[1]$ are given the correct values.

Figures VII.1 and VII.2 present the recovery error and computation time of all algorithms as a function of the number of observations $M$ for fixed noise level $\sigma=1$. Of course, the oracle who knows the shifts of the observations is unbeatable. Algorithms 2 and 3 outperform all invariant approach methods. The inferior performance of the SDP might be explained by the fact that we are minimizing a smooth non-linear objective. This is in contrast to SDPs with linear or piecewise linear objectives which tend to promote “simple” (i.e., low rank) solutions [74], [75, Remark 6.2]. Additionally, while this is not depicted on the figure, we note that for $\sigma=0$ all methods get exact recovery up to machine precision. EM outperforms the best invariant features approaches by a factor of 3, at the cost of being significantly slower. For large $M$, the best invariant features approaches are faster than EM by a factor of 25. Note, however, that for $M$ up to about 300, EM is faster than the other algorithms. For invariant features approaches (aside from the SDP), almost all of the time is spent computing the bispectrum estimator, while inverting the bispectrum is relatively cheap.

Figures VII.3 and VII.4 show the recovery error and computation time as a function of the noise level $\sigma$ with $M=10,\!000$ observations. Surprisingly, for high noise level $\sigma\gtrsim 3$, the invariant features algorithms outperform EM.

VIII Conclusions and perspective

The goal of this paper is twofold. First, we have suggested a new approach for the MRA problem based on features that are invariant under translations. This technique enables us to deal with any noise level as long as we have access to enough measurements and particularly it achieves the sample complexity of MRA. The invariant features approach has low computational complexity and it requires less memory with respect to alternative methods, such as EM. If one wants to have a highly accurate solution, it can therefore be used to initialize EM.

A main ingredient of the invariant features approach is estimating the signal’s Fourier phases by inverting the bispectrum. Hence, the second goal of this paper was to study algorithms for bispectrum inversion. We have proposed a few algorithms for this task. In the presence of noise, the non-convex algorithms on the manifold of phases, namely, Algorithms 2 and 3, perform the best. Empirically, these algorithms have a remarkable property: despite their non-convex landscape, they appear to converge to the target signal from random initialization. We provide some analysis for Algorithm 2 but this phenomenon is not well understood.

Our chief motivation for this work comes from the more involved problem of cryo-EM. In cryo-EM, a 3D object is estimated from its 2D projections at unknown rotations in a low SNR environment. One line of research for the object recovery is based on first estimating the unknown rotations [76, 77, 78]. However, the rotation estimation is performed in a very noisy environment and therefore might be inaccurate. An interesting question is to examine whether the 3D object can be estimated directly from the acquired data using features that are invariant under the unknown viewing directions [79].

Acknowledgement

The authors are grateful to Afonso Bandeira, Roy Lederman, William Leeb, Nir Sharon and Susannah Shoemaker for many insightful discussions. We also thank the reviewers for their useful comments, and particularly the anonymous reviewer who proposed a simpler proof for Lemma IV.2.

-A Stable inversion of bispectrum

We first give a proof of Propositions II.2.

Let $U$ denote the set of signals of length $N$ (either real or complex) whose DFTs are non-vanishing. Let $G=\{R_{0},\ldots,R_{N-1}\}$ denote the discrete group of cyclic shifts over signals of length $N$. In MRA, the parameter space is the quotient space $U/G$, turned into a metric space with distance

[TABLE]

In the above notation, $[x]=\{R_{0}x,\ldots,R_{N-1}x\}$ is the equivalence class of signal $x$. Let $V\subset\mathbb{C}^{N\times N}$ denote the set of matrices whose entries are nonzero, equipped with the Frobenius norm distance.

We now construct a function $\psi\colon V\to U/G$. Crucially, $\psi$ is designed such that if $B\in V$ is the bispectrum of a signal $x$ in $U$, then $\psi(B)=[x]$: it inverts bispectra. Our purpose for building $\psi$ is to establish local Lipschitz continuity.

For a given input $B\in V$, for each $r$ in $\{0,\ldots,N-1\}$, we construct the moduli and phases of $y_{r}\in\mathbb{C}^{N}$ separately. The inverse DFTs $\{x_{0},\ldots,x_{N-1}\}$ form an equivalence class in $U/G$: this is the output $\psi(B)$. Specifically:

For all $r$, set $|y_{r}[0]|=\sqrt[3]{|B[0,0]|}$ and $\mathrm{phase}(y_{r}[0])=\mathrm{phase}(B[0,0])$. 2. 2.

For all $r$, the moduli $|y_{r}[k]|$ for $k=1\ldots N-1$ are:

[TABLE] 3. 3.

Assume we are given a reference $\theta_{0}\in\mathbb{R}$ (specified later, independent of input $B$). For a unit-modulus complex number $z$, we define the principal $N$th root of $z$ as $z^{1/N}=e^{i\theta/N}$ with $\theta\in(\theta_{0}-\pi,\theta_{0}+\pi]$ such that $z=e^{i\theta}$. 4. 4.

With this definition of $N$th root, the phase of the first non-DC component $\mathrm{phase}(y_{r}[1])$ is

[TABLE]

(Notice $B[1,2]$ appears twice.) Each value of $r$ assumes one of $N$ possible $N$th roots. 5. 5.

The phases of $y_{r}[k]$ for $k=2\ldots,N-1$ are defined recursively, separately for each $r$:

[TABLE] 6. 6.

Define $x_{r}$ as the inverse DFT of $y_{r}$; it is easily checked that $x_{r}=R_{r}x_{0}$, so that $\{x_{0},\ldots,x_{N-1}\}$ forms an equivalence class in $U/G$.

We now argue that $\psi$ is locally Lipschitz continuous. The construction above implicitly defines a function $x_{0}\colon V\to U$, parameterized by a reference $\theta_{0}$. For a given $B_{0}\in V$, pick the reference such that $e^{i\theta_{0}}=\mathrm{phase}\left(B_{0}[N-1,1]B_{0}[1,2]B_{0}[1,2]\cdots B_{0}[1,N-1]\right)$. There exists a neighborhood $W$ of $B_{0}$ such that all $B\in W$ obey

$|B[k,\ell]|\geq\frac{1}{2}|B_{0}[k,\ell]|$ for all $k,\ell$, and 2. 2.

$\Re\{e^{-i\theta_{0}}\mathrm{phase}\left(Z\right)\}\geq-\frac{1}{2}$,

with $Z=B[N-1,1]B[1,2]B[1,2]\cdots B[1,N-1]$.

In that neighborhood, $\mathrm{phase}(y_{r}[1])$ as defined in (-A.1) is a smooth function of $B$. All other operations involved in computing $\psi(B)$ are smooth in $W$ as well, since entries of $B$ are nonzero. Thus, $x_{0}\colon W\to U$ is smooth. In particular, there exists a neighborhood $\bar{W}\subset W$ of $B_{0}$ such that $x_{0}$ is Lipschitz in $\bar{W}$: there exists $L>0$ such that

[TABLE]

Thus, for all $B_{0}\in V$ there exists a neighborhood $W$ of $B_{0}$ in $V$ and $L>0$ such that

[TABLE]

confirming $\psi$ is locally Lipschitz continuous. In particular, if $B=B_{0}$ is the bispectrum of $x\in U$ so that $\psi(B)=[x]$, then there exists $\delta>0$ small enough and $L>0$ such that if $\|B^{\prime}-B\|_{\mathrm{F}}\leq\delta$, then $[\hat{x}]=\psi(B^{\prime})$ obeys

[TABLE]

which shows $x$ can be estimated (up to shift) with finite error from a sufficiently good bispectrum estimator.

We now give a proof of Proposition II.3.

The estimator $\hat{B}_{x}$ has expectation equal to $B_{x}$ and variances on its individual entries are bounded by $c\frac{\sigma^{2}+\sigma^{6}}{M}$ for some $c=c(x)$.

Chebyshev’s inequality states that for a random variable $X$ with mean $\mu$ and variance $\sigma^{2}$ the probability that $|X-\mu|$ exceeds $t\sigma$ is bounded by $1/t^{2}$. We have $N^{2}$ random variables. By a union bound (which, importantly, does not require independence), we have

[TABLE]

Pick $t=\frac{N}{\sqrt{1-p}}$ so that the probability is at least $p$. Now pick $M$ such that $t\sqrt{c\frac{\sigma^{2}+\sigma^{6}}{M}}\leq\delta^{\prime}$. Explicitly,

[TABLE]

This shows that for any $p$ and $\delta^{\prime}$ we can pick $C$ such that with $M\geq C(\sigma^{2}+\sigma^{6})$ we have that $\hat{B}_{x}$ approximates $B_{x}$ within $\delta^{\prime}$ (entry-wise) with probability at least $p$. Pick $\delta^{\prime}=\delta/N$ to conclude. (We stress that this expression for $C$ is suboptimal.)

-B Analysis of optimization over phases for real signals

As aforementioned, if the signal under consideration is real, then the phases $\tilde{y}$ of its Fourier transform obey symmetries. These should be exploited, and the estimator should satisfy the same symmetries. Specifically,

[TABLE]

This notably implies that $T(\tilde{y})$ is Hermitian, and as a result that $\tilde{B}$ is Hermitian. It is then also sensible to take $W$ symmetric. Furthermore, this implies that $\tilde{y}[0]$ is real, so that $\tilde{y}[0]=\pm 1$, and, if $N$ is even, that $\tilde{y}[N/2]=\pm 1$ as well. In the latter case, observe that $\tilde{y}[N/2]$ changes sign when the signal is time-shifted by one index, so that we may fix $\tilde{y}[N/2]=1$ without loss of generality, if need be—this is discussed more below. Without loss of generality, let us assume $\tilde{y}[0]=1$ (in the MRA framework, this would be estimated from $\hat{\mu}_{x}$).

If the signal is real, the conditions of Lemma V.1 can be slightly alleviated so that the global optima of (III.4) corresponds exactly to $\tilde{y}$ as follows:

Lemma -B.1**.**

For $N\geq 5$, let $x\in\mathbb{R}^{N}$ be a real signal whose DFT $y$ is nonzero for frequencies $k$ in $\{1,\ldots,K\}\cup\{N-1,\ldots,N-K\}$, possibly also for $k=0$, and zero otherwise. Up to integer time shifts, $x$ is determined exactly by its bispectrum $B$ provided $\frac{N}{3}\leq K\leq\frac{N-1}{2}$. Furthermore, the global optima of (III.4) (with conjugate-reflection constraints) correspond exactly to the relevant phases of $y$—up to the effects of integer time shifts—provided $W[k,\ell]$ is positive when $B[k,\ell]\neq 0$.

Proof.

See Appendix -K. ∎

In the following, it is helpful to introduce visual notation for the symmetries. Given any vector $u\in\mathbb{C}^{N}$, consider the following linear operations:

[TABLE]

and

[TABLE]

Using this notation, the symmetries of $z$, namely, $z[N-k]=\overline{z[k]}$ for all $k$, can be written $z_{\uparrow}=\overline{z_{\downarrow}}$. As a result, for odd $N$ the phase estimation problem lives on the following submanifold of $\mathcal{M}$:

[TABLE]

If $N$ is even, the submanifold takes a slightly different form:

[TABLE]

The corresponding optimization problem is

[TABLE]

To carry out the optimization of (-B.2), we follow the same protocol as for the complex case, namely, we obtain expressions for the Riemannian gradient and the Riemannian Hessian of the cost function $f$ restricted to the manifold $\mathcal{M}_{\mathbb{R}}$, at which point we will be in a position to use standard algorithms.

The first step is to obtain an expression for the orthogonal projector from $\mathbb{C}^{N}$ to the tangent spaces of $\mathcal{M}_{\mathbb{R}}$. These tangent spaces are readily obtained by differentiating the constraints:

[TABLE]

This implicitly forces $\dot{z}[0]=0$ and, if $N$ is even, $\dot{z}[N/2]=0$. The orthogonal projection from $\mathbb{C}^{N}$ to that linear space is

[TABLE]

One can check that, for any $u$, $\left(\operatorname{Proj}_{z}^{\mathbb{R}}u\right)[k]=0$ if $k=0$ and if $k=N/2$ for $N$ even, as expected. As for the complex case, these steps yield explicit expressions for the Riemannian gradient and Hessian of $f$ on $\mathcal{M}_{\mathbb{R}}$:

[TABLE]

and

[TABLE]

To reach the last equality, we used the fact that $\operatorname{Proj}_{z}^{\mathbb{R}}(\alpha\circ z)=0$ for any real vector $\alpha$, so that one of the terms vanished under the projection.

Interestingly, for $z$ restricted to $\mathcal{M}_{\mathbb{R}}$, the cost function and its derivatives simplify. Indeed, $T(z)$ is now Hermitian, so that $M(z)$ is Hermitian (using a symmetric $W$) and (V.5) becomes:

[TABLE]

The last simplification of the Hessian follows from this result, under symmetry assumptions which are valid in the real case.

Lemma -B.2**.**

Under the same symmetry conditions as in Lemma V.2, if additionally $W$ and $\tilde{B}$ are Hermitian, then, if $z\in\mathcal{M}_{\mathbb{R}}$ and $\dot{z}\in\mathrm{T}_{z}\mathcal{M}_{\mathbb{R}}$, it holds that $M(z)\dot{z}=M(\dot{z})z$.

Proof.

See Appendix -L. ∎

It remains to define a retraction on $\mathcal{M}_{\mathbb{R}}$. An obvious choice is:

[TABLE]

where $\operatorname{Retr}_{z}$ was defined in (V.6). Indeed, it is a simple exercise to check that $\operatorname{Retr}_{z}(\dot{z})$ is in $\mathcal{M}_{\mathbb{R}}$ for any $z\in\mathcal{M}_{\mathbb{R}}$ and $\dot{z}\in\mathrm{T}_{z}\mathcal{M}_{\mathbb{R}}$:

[TABLE]

Empirically, with exact data $\tilde{B}$ and unit weights $W$, we find that the parity of $N$ has an important effect on the nature of second-order critical points of (-B.2). If $N$ is odd and $z[0]=\tilde{y}[0]$, for $N=3,5,7,9$, we find empirically that $f\colon\mathcal{M}_{\mathbb{R}}\to\mathbb{R}$ has $N$ second-order critical points which are all global optima; they correspond to the phases of the DFT of all possible time-shifts of the unknown signal. Hence, any second-order critical point leads to exact recovery. On the other hand, if $N$ is even, we find that only $N/2$ second-order critical points are global optima. If $N=4,6,8$, there are no other second-order critical points. If $N=10,12$, there are $N/2$ additional second-order critical points. These are strict, non-global local optima. Interestingly, if one flips the sign of $z[N/2]$, one recovers the phases of one of the other $N/2$ possible time-shifts of the unknown signal. This is why, for even $N$, we recommend the following: (i) compute a second-order critical point $z$ of (-B.2) from some initial guess; (ii) let $z^{\prime}$ be $z$ with the sign of $z[N/2]$ flipped; (iii) if $f(z^{\prime})>f(z)$, run the optimization again with $z[N/2]$ flipped in $\mathcal{M}_{\mathbb{R}}$, this time starting at $z^{\prime}$. In the noiseless case, empirically, $z^{\prime}$ will already be optimal.

Of course, the parameterization of $\mathcal{M}_{\mathbb{R}}$ as proposed here is redundant. Given the symmetries, one could alternatively choose to work with only the phases $z[1]$ to $z[\lfloor\frac{N-1}{2}\rfloor]$, which are sufficient to determine all of $z\in\mathcal{M}_{\mathbb{R}}$. This is certainly computationally advantageous. The choice to work with a redundant parameterization above is motivated by the simpler exposition it allows. There is no conceptual or technical difficulty in implementing the above with a non-redundant parameterization.

-C Proof of Lemma III.1

Defining $z_{\theta}$ as $z_{\theta}[k]=z[k]e^{i\theta g(k)}$ with $g(k)=k\bmod N$, the statement to prove is equivalent to $f(z)=f(z_{\theta})$ for any $\theta$ which is an integer multiple of $\frac{2\pi}{N}$.

Let $\omega=e^{i\theta}$ and let $u\in\mathbb{C}^{N}$ be the vector defined by $u[k]=\omega^{g(k)}$. Then, $z_{\theta}=z\circ u=\mathrm{diag}(u)z$ so that

[TABLE]

Using (III.3) for the definition of $M$ we get

[TABLE]

Hence,

[TABLE]

All dependence on $\theta$ now resides in the matrix $\mathrm{diag}(u)^{*}\overline{T(u)}\mathrm{diag}(u)$. Its entries obey

[TABLE]

For $(k_{1},k_{2})$ ranging from 0 to $N-1$, the exponent $g(k_{2})-g(k_{1})-g(k_{2}-k_{1})$ evaluates to 0 if $k_{2}\geq k_{1}$ and to $N$ otherwise. In the first case, $\omega^{0}=1$. In the second case, $\omega^{N}=e^{i\theta N}=1$ owing to the assumption that $\theta$ is an integer multiple of $\frac{2\pi}{N}$. Consequently,

[TABLE]

and indeed $f(z)=f(z_{\theta})$.

-D Averaging over phases

Let $a_{1},a_{2},\dots,a_{K}$ be complex numbers with unit modulus. We define the average over the set of phases (i.e., the group SO(2)) as the solution of

[TABLE]

Expanding the square modulus, the problem is equivalent to

[TABLE]

The last expression is maximized if and only if $\overline{z}\sum_{k=1}^{K}a_{k}$ is real and non-negative, i.e., $z$ is the phase of $\sum_{k=1}^{K}a_{k}$. Explicitly, if $\sum_{k=1}^{K}a_{k}\neq 0$ then

[TABLE]

Otherwise, any unit-modulus $z$ is optimal.

-E Proof of Lemma IV.2

By expanding $u[1]$, we get

[TABLE]

where the last equality holds since $\hat{u}$ is real. Therefore, since $\hat{u}$ is non-negative, we get

[TABLE]

This equality holds if and only if $\hat{u}[k]=0$ whenever $\Re\left\{e^{2\pi ik/N}\right\}<1$. This implies in turn that $\hat{u}$ is a delta function and thus $u$ is a constant.

-F Proof of Theorem IV.3

In the noiseless case, plugging in the correct signal into (IV.2) satisfies the constraints and attains 0 for the objective value, which is clearly minimal. Thus, any global optimum $(Z,z)$ must have objective value 0, meaning $\tilde{B}\circ\overline{T(z)}=Z$.

Using Schur’s complement, the SDP constraint $Z\succeq zz^{*}$ can be written as

[TABLE]

Let

[TABLE]

where $\mbox{diag}\left(\tilde{y}\right)$ is a diagonal matrix with the entries of $\tilde{y}$. Observe that $Q\succeq 0$ implies $U^{*}QU\succeq 0$. Recall that each entry of $\tilde{y}$ has modulus one. Then, using $\tilde{B}=\tilde{y}\tilde{y}^{*}\circ T(\tilde{y})$, we get

[TABLE]

where $u:=\overline{\tilde{y}}\circ z$. This in turn implies that $T(\overline{u})\succeq 0$, hence, by properties of circulant matrices, $\hat{u}$ (the DFT of $u$) is non-negative. Furthermore, since $z[0]=\tilde{y}[0]$ and $z[1]=\tilde{y}[1]$, it follows that $u[0]=u[1]=1$. So, by Lemma IV.2 we conclude that $u[k]=1$ for all $k$ and $z=\tilde{y}$. This implies immediately $Z=zz^{*}=\tilde{y}\tilde{y}^{*}$.

-G Proof of Lemma V.1

We start by characterizing the global optima of (III.4). Recall that

[TABLE]

Defining $u=\overline{z}\circ\tilde{y}$, the objective function of (III.4) simplifies to

[TABLE]

Clearly, $f(\tilde{y})$ is maximal over $\mathcal{M}$. Any other global optimum must attain the same objective value. By our assumptions on $W$, this implies for all ${k,\ell,\ell-k\leavevmode\nobreak\ \in\{1,\ldots,K\}\bmod N}$ that

[TABLE]

Consider $k=1$ and $\ell=2,\ldots,K$; this easily leads to $u[k]=u[1]^{k}$ for $k\in\{1,\ldots,K\}$. Now set $k=K,\ell=1$. Using periodicity of indexing, since $\ell-k=1-K=N+1-K\bmod N$, and using the conditions $K\geq\frac{N+1}{2}$ and $K\leq N-1$, we have $N+1-K\in\{2,\ldots,K\}$, so that (-G.1) applies. On one hand, we find $u[\ell-k]=u[N+1-K]=u[1]^{N+1-K}$, and on the other hand we find $\overline{u[\ell-k]}=u[k]\overline{u[\ell]}=u[1]^{K-1}$. Thus, $u[1]^{N}=1$. If we let $u[1]=e^{i\theta}$ for some $\theta\in\mathbb{R}$, it follows that $\theta=m\frac{2\pi}{N}$ for some integer $m$. If $x$ has zero mean, $u[0]$ is irrelevant and can be set to $1$; otherwise, equation (-G.1) also holds with $k=\ell=0$, so that $u[0]=1$. Recalling $u=\overline{z}\circ\tilde{y}$, it follows that if $z$ is optimal, then $z[k]=\tilde{y}[k]e^{-2\pi im\frac{k}{N}}$ for relevant $k$’s. For $m$ in $\{0,\ldots,N-1\}$, these are exactly the phases of the DFTs of all $N$ circular time-shifts of $x$, which are indeed all optimal.

It remains to show that the amplitudes of the DFT of $x$ (the power spectrum of $x$) can be recovered from $B$. If $x$ has nonzero mean, the power spectrum can be read off the diagonal of $B$. If $y[0]=0$, the power spectrum can still be recovered. Indeed, consider $\log|B[k,\ell]|$:

[TABLE]

These provide linear equations in $\log|y[1]|,\ldots,\log|y[K]|$. It suffices to collect $K$ independent ones. Considering in order equations with $k=1,\ell=2,\ldots,K$

followed by $k=2,\ell=4$, we get a structured linear system. For example, with $K=5$:

[TABLE]

The determinant of the matrix is $(-1)^{K}6$, proving $|y[k]|$ can be recovered from $B$. Together with $y[0]=0$ and the phases, this is sufficient to recover $x$ up to global time shift.

-H Proof of identity (V.2)

Observe that $T$ can be expressed as

[TABLE]

where $T_{k}$ is a circulant matrix with ones in its $k$th (circular) diagonal and zero otherwise as defined in (V.3). Using (-H.1) and (III.3),

[TABLE]

Since this must hold for all $z$ and $X$, by identification on the last line,

[TABLE]

-I Proof of Lemma V.2

On one hand, note that

[TABLE]

On the other hand, we have from (V.2) and (V.3) that

[TABLE]

where we used the assumed symmetries of $\tilde{B}$ and $W$ in the last step only. The expressions match, concluding the proof.

-J Proof of Lemma V.5

Since $z$ is critical, $\nabla f(z)\circ\overline{z}$ is real. It remains to show that it is positive. For this, we will use the second-order condition. Under the symmetry assumptions of Lemma V.2 which hold a fortiori in the noiseless case, the term $\left\langle{\dot{z}},{\nabla^{2}f(z)[\dot{z}]}\right\rangle$ can be developed into

[TABLE]

where the last equality follows Lemma V.2. Hence,

[TABLE]

For some index $k$, consider the tangent vector $\dot{z}=(iz[k])e_{k}$ in $\mathrm{T}_{z}\mathcal{M}$, where $e_{k}\in\mathbb{R}^{N}$ is the $k$th canonical basis vector. Since $z$ is second-order critical, we have the inequality

[TABLE]

Plugging in the expression of $\dot{z}$, this is equivalent to

[TABLE]

The left hand side develops as follows:

[TABLE]

where in the last equality we substituted $\ell^{\prime}=\ell-k$. Using the definition of $\tilde{B}$ as in (II.5) we then get

[TABLE]

In particular, for $k=0$, this simplifies to $2\left\langle{\tilde{y}[0]},{z[0]}\right\rangle$. Then, the inequality is

[TABLE]

Let the sum in the right hand side be denoted by $t\in\mathbb{C}$. Clearly, $|t|\leq N$ and $2\left\langle{\tilde{y}[0]},{z[0]}\right\rangle\geq-2$, so

[TABLE]

Since $z$ is critical, the right hand side is real, so that it is equal to either $|2N+t|$ or $-|2N+t|$ (any real number is either its absolute value or the opposite). For contradiction, let us assume it is equal to $-|2N+t|$. Then,

[TABLE]

This is impossible if $N>2$. Hence, $(\tilde{y}[0]\overline{z[0]})\left(2N+t\right)=|2N+t|$, so that $\overline{\tilde{y}[0]}z[0]=\mathrm{phase}(2N+t)$. Using $|t|\leq N$, it is a simple exercise to determine that $\left\langle{\tilde{y}[0]},{z[0]}\right\rangle\geq\frac{\sqrt{3}}{2}$. Turning back to general $k$ and using (-J.1) and (-J.2), it follows that

[TABLE]

Manually checking the statement for $N=1,2$ concludes the proof.

-K Proof of Lemma -B.1

The proof is identical to that of Lemma V.1, up to a few differences we highlight. With the same definition of the vector $u$, the identity $u[k]\overline{u[\ell]}=\overline{u[\ell-k]}$ still holds for $k,\ell$ such that $k,\ell,\ell-k$ are in $\{1,\ldots,K\}\cup\{N-1,\ldots,N-K\}\bmod N$. Using real symmetry, $u[-k]=\overline{u[k]}$ for all $k$, hence the identity also reads $u[k]=u[\ell]u[k-\ell]$.

Using this rule for $k=2,\ldots,K$ and fixed $\ell=1$, we easily get $u[k]=u[1]^{k}$ for $k\in\{1,\ldots,K\}$. We now consider the rule for $k=N-K,\ell=K$. Since $K\geq\frac{1}{3}N$, it is clear that $k-\ell=N-2K\leq\frac{1}{3}N\leq K$. Since $K\leq\frac{N-1}{2}$, it also holds that $k-\ell\geq 1$. As a result, $u[N-K]=u[1]^{N-K}$. Likewise, using real symmetry and indexing modulo $N$, $u[N-K]=u[-K]=\overline{u[K]}=u[1]^{-K}$. Combining the two, it follows that $u[1]^{N}=1$. One can then conclude as in the proof of the previous lemma. The magnitudes of the DFT can also be recovered, following the same procedure as in Lemma V.1: if $y[0]\neq 0$, read the power spectrum off the diagonal; otherwise, obtain $|y[1]|,\ldots,|y[K]|$ via the same linear system and use $|y[N-k]|=|y[k]|$ for $k=1\ldots K$.

-L Proof of Lemma -B.2

For ease of notation, let $A=W^{(2)}\circ\tilde{B}$. By the assumptions of Lemma V.2, we know that $A[\ell-k,\ell]=A[k,\ell]$ for all $k,\ell$. Since $A$ is now also Hermitian, we have

[TABLE]

Hence, with the change of variable $\ell^{\prime}=k-\ell$ and using both $z[-k]=\overline{z[k]}$ and $\dot{z}[-k]=\overline{\dot{z}[k]}$ for any $k$ since $z\in\mathcal{M}_{\mathbb{R}}$ and $\dot{z}\in\mathrm{T}_{z}\mathcal{M}_{\mathbb{R}}$:

[TABLE]

This concludes the proof.