Expectation-maximization for low-SNR multi-reference alignment

TL;DR

This work analyzes the EM algorithm for multi-reference alignment in the challenging low-SNR regime, revealing a two-phase local convergence near the ground truth and a fundamental iteration bottleneck that scales as $T \gtrsim \mathrm{SNR}^{-2}$. It uncovers initialization-driven biases (Einstein from Noise) and finite-sample instabilities (Ghost of Newton), showing that Fourier phases can persist from initialization while magnitudes decay slowly, and that phase drift accumulates with sample size and iterations. The paper provides a detailed population-level Jacobian analysis in the low-SNR limit, a spectral decomposition in Fourier space, and a precise iteration-complexity bound, complemented by a thorough finite-sample analysis that identifies sample-size thresholds $n \gtrsim \mathrm{SNR}^{-3}$ and phase-drift phenomena $\propto 1/n$. Collectively, these results expose fundamental computational and initialization limitations of EM for MRA in practical, noisy settings and suggest mitigation directions, including mini-batching and potential second-order methods. The findings have implications for real-world applications in cryo-EM and related latent-group models where alignment under noise is essential yet challenging.

Abstract

We study the multi-reference alignment (MRA) problem of recovering a signal from noisy observations acted on by unknown random circular shifts. While the information-theoretic limits of MRA are well characterized in many settings, the algorithmic behavior at low signal-to-noise ratio (SNR), the regime of practical interest, remains poorly understood. In this paper, we analyze the expectation-maximization (EM) algorithm, a widely used method for MRA, and characterize its convergence dynamics and initialization dependence in the low-SNR limit. On the convergence side, we prove a two-phase phenomenon near the ground truth as $\mathrm{SNR}\to 0$: an initial contraction with error decaying as $\exp(-\, \mathrm{SNR} \cdot t)$ followed by a much slower phase scaling as $\exp(- \,\mathrm{SNR}^2 \cdot t)$, where $t$ is the iteration number. This yields an iteration-complexity lower bound $T \gtrsim \mathrm{SNR}^{-2}$ to reach a small fixed target accuracy, revealing a severe computational bottleneck at low SNR. We also identify a finite-sample instability, which we term \emph{Ghost of Newton}, in which EM initially approaches the ground truth but later diverges, degrading reconstruction quality. On the bias side, we analyze EM in the noise-only setting ($\mathrm{SNR}=0$), a regime referred to as Einstein from Noise, to highlight its pronounced sensitivity to initialization. We prove that the EM map preserves the Fourier phases of the initialization across all iterations, while the corresponding Fourier magnitudes contract toward zero at a slow rate of $(1+T)^{-1/2}$. Consequently, although the amplitudes vanish in the limit of $T \to \infty$ iterations, the reconstructed structure continues to reflect the geometry encoded by the template's Fourier phases. Together, these results expose fundamental computational and initialization-driven limitations of EM for MRA in the low-SNR regime.

Figures (11)

• Figure 1: Population EM dynamics in the low-SNR regime.(a) Two-phase convergence at positive $0 < \mathrm{SNR} \ll 1$. The plot shows the normalized population MSE $\|\hat{x}^{(t)}-x^\star\|_2^2/\|x\|_2^2$ versus iteration $t$ for several low SNR values (annotated on the curves). The population EM map exhibits a clear two-phase behavior: an initial moderate-rate geometric decay, $\mathrm{MSE}\propto \exp\left(-c \, \mathrm{SNR} \cdot t\right)$, followed by a slow phase where the error decreases much more slowly, $\mathrm{MSE}\propto\exp\left(- c' \, \mathrm{SNR}^2 \cdot t\right)$. (b) Fourier magnitudes contraction to 0 at $\mathrm{SNR}=0$. The Fourier magnitudes $|\mathsf{\hat{X}}^{(t)}[k]|^2$ are shown for a representative non-mean Fourier mode under the population EM operator. The dashed curve is the closed-form prediction from Theorem \ref{['thm:low-mag-rate']}, $|\mathsf{\hat{X}}^{(t)}[k]|^2/|\mathsf{\hat{X}}^{(0)}[k]|^2 = (1+2t|\mathsf{\hat{X}}^{(0)}[k]|^2)^{-1}$. The agreement corroborate that, in the regime of $\mathrm{SNR} = 0$, Fourier magnitudes decay toward zero at a rate $\asymp(1+t)^{-1}$, while the corresponding Fourier phases remain frozen along the population trajectory (not shown in the Figure). Simulation parameters: Population expectations defining the EM update are evaluated numerically via Gauss-Hermite quadrature over the Gaussian noise with a fixed signal dimension $d=5$.
• Figure 2: Empirical setup used in this work and reconstruction regimes.(a) Motivating application: Illustration of the MRA model \ref{['eqn:mainModel1D']} in a two-dimensional setting. The goal is to reconstruct a 2D Newton image from noisy, randomly shifted observations, starting from an initial template of Einstein's image. We study the expectation-maximization (EM) algorithm as the primary method analyzed in this work, and include a hard-assignment method as a baseline for comparison. (b) Reconstruction performance for both methods is quantified by the MSE, defined in \ref{['eqn:mseDef']}, and plotted as a function of SNR under cyclic 2D shifts. Four distinct regimes are observed: 1) Einstein from Noise: at very low SNR, the reconstruction resembles the initial template (Einstein); 2) Ghost of Newton: at moderately low SNR, the EM algorithm initially converges to the true structure (Newton), but in later iterations, the noise term causes the reconstruction to gradually deteriorate, leading to increasingly corrupted estimates (see Figure \ref{['fig:4']} for further illustration). 3) a transition region, 4) high-SNR regime, where alignment is accurate and performance approaches the baseline case, following the scaling law $\mathsf{MSE} \propto 1/\mathrm{SNR}$. All experiments use $d = 64 \times 64$ images, $n = 2 \times 10^4$ observations, $T = 200$ iterations, and average results over 30 Monte Carlo trials per data point.
• Figure 3: Empirical demonstration of the Jacobian's second-order expansion, block spectrum, and spectral radius.(a) Numerical support for the second-order expansion in \ref{['eq:Jhat_explicit_matrix_prop']}: the operator-norm discrepancy between the exact Jacobian $\widehat{J}(\beta)$ and its second-order approximation $\widehat{J}_2(\beta)$ (given by \ref{['eq:Jhat_explicit_matrix_prop']}) decays proportionally to $\beta^4$, consistent with a remainder term of order $O(\beta^4/\sigma^4)$. (b) Numerical support for the blockwise spectral characterization in Corollary \ref{['cor:block-spectral']}. In the $\{k,-k\}$ contracting blocks, the eigenvalues satisfy $\lambda_k(u_k)=1-2(\beta^2/\sigma^2)|\widehat{V}_k|^2+O(\beta^4/\sigma^4)$; correspondingly, the plotted deviation $\max_k\bigl|\lambda_k(u_k)-\bigl(1-2(\beta^2/\sigma^2)|\widehat{V}_k|^2\bigr)\bigr|$ exhibits slope $4$ on the log--log scale. Along the flat directions, $\lambda_k(w_k)=1-O(\beta^4/\sigma^4)$, and the quantity $\max_k|\lambda_k(w_k)-1|$ also scales as $\beta^4$, in agreement with the corollary. (c) Scaling of the spectral radius: the spectral gap $1-\rho({J}(\beta))$ follows $O(\beta^4/\sigma^4)$ and displays slope $4$ on the log--log scale, consistent with Corollary \ref{['cor:block-spectral']}(c). All curves are computed using Gauss--Hermite quadrature to evaluate the population Jacobian, with signal dimension $d=5$ and noise level $\sigma^2=1$.
• Figure 4: The Ghost of Newton phenomenon.(a) MSE, as defined in \ref{['eqn:mseDef']}, between the EM estimator at iteration $t$ and the true signal $x^\star$, plotted as a function of iteration $t$. Initially, the algorithm converges toward the true signal of Newton, but after several iterations, it begins to diverge, ultimately yielding a noisy reconstruction. (b) The EM update consists of two components: the aligned signal, $\hat{x}_{\mathrm{signal}}^{(t+1)}$, and the aligned noise, $\hat{x}_{\mathrm{noise}}^{(t+1)}$. In the Ghost of Newton regime, the aligned noise, despite originating from pure noise, resembles a noisy version of the Newton image. This misleading signal steers the algorithm away from the ground truth. (c) Visualization of the aligned signal $\hat{x}_{\mathrm{signal}}^{(t+1)}$, aligned noise $\hat{x}_{\mathrm{noise}}^{(t+1)}$, and the total estimate $\hat{x}^{(t+1)}$ at representative iterations. In the early stages, the aligned noise and the estimate resemble the initial template of Einstein. As the iterations progress, the Newton image gradually emerges. However, after approximately 50 iterations, the influence of the aligned noise dominates, steering the reconstruction away from Newton and resulting in a noisy and degraded estimate. The simulation corresponds to a fixed signal-to-noise ratio of $\mathrm{SNR} = 5 \times 10^{-3}$, with an image size of $d = 64 \times 64$ and $n = 2 \times 10^4$ observations.
• Figure 5: The Impact of the number of observations and SNR on the Ghost of Newton phenomenon.(a, b) Effect of increasing the number of observations at a fixed SNR. Parameters: image size $d = 64 \times 64$, SNR $= 5 \times 10^{-3}$. A larger number of observations reduces the impact of noise and mitigates the Ghost of Newton effect. (c, d) Effect of increasing the SNR at a fixed number of observations. Parameters: image size $d = 64 \times 64$, $n = 2 \times 10^4$ observations. Higher SNR leads to better alignment and convergence to the true signal. The $\text{MSE}$ plotted in both cases is defined in \ref{['eqn:mseDef']}, while the log-likelihood corresponds to \ref{['eqn:logLikelihoodMRA']}, with the constant final term omitted as it does not depend on the iteration number. In both scenarios, a sharp increase in the log-likelihood ($\mathcal{L}_{\mathsf{MRA}}$) is observed at the region where the estimator begins to diverge from Newton’s image---this transition is indicated by dashed lines. Notably, this pattern consistently appears across different numbers of observations and SNR levels.

Theorems & Definitions (103)

• Theorem 1.1: Informal: Iteration complexity of EM in MRA
• Proposition 2.1: Global convergence at high-SNR
• Definition 3.1: Jacobian and spectral radius
• Theorem 3.3: Local linear convergence of the population EM operator
• Proposition 3.4: Second-order Jacobian expansion in the low-SNR regime
• Corollary 3.5: Jacobian spectral decomposition
• Theorem 3.6: Two-phase convergence
• Remark 3.7: Basin of attraction
• Corollary 3.8: Tail iteration complexity of EM in MRA
• Proposition 3.9: Local expansion of Fourier magnitudes/phases around $x^\star$