A Fast Approximate Maximum Likelihood Estimator for Low SNR Multi-Reference Alignment

Technical Summary: A Fast Approximate Maximum Likelihood Estimator for Low SNR Multi-Reference Alignment

Problem Statement

The paper addresses the Multi-Reference Alignment (MRA) problem, specifically over the special orthogonal groups SO(2) and SO(3). The objective is to recover an unknown signal from $n$ noisy observations, where each observation has been independently transformed by an unknown group action (rotation) and corrupted by additive Gaussian noise.

The primary challenge arises in the low signal-to-noise ratio (SNR) regime. In this setting, individual observations are so heavily corrupted that estimating the specific rotation applied to each observation is statistically infeasible. While the standard approach for computing the Maximum Likelihood Estimator (MLE) is the Expectation-Maximization (EM) algorithm, it suffers from two critical limitations in low-SNR settings:

1. Computational Prohibitive Cost: EM requires multiple iterations over the entire dataset. As SNR decreases, convergence becomes extremely slow, often requiring a number of iterations that grows rapidly.
2. Non-Convexity: The likelihood landscape is non-convex, making EM prone to converging to spurious local minima unless initialized with a high-quality estimate.

Alternative methods based on invariant features (e.g., moments, bispectra) avoid latent rotation estimation and offer single-pass consistency. However, in high-dimensional settings like SO(3), recovering these invariants requires solving large systems of polynomial equations, which becomes computationally expensive.

Methodology

The authors propose a fast approximate MLE derived from a low-SNR Taylor expansion of the marginalized log-likelihood function. The core insight is that in the low-SNR limit, the likelihood landscape simplifies, allowing for a closed-form, non-iterative solution.

The algorithm operates via a frequency marching strategy, recovering signal components sequentially from lower to higher frequencies (or angular degrees):

1. Preliminary Estimates:

   • SO(2): The zero-frequency coefficient is estimated by averaging. The magnitudes of higher-frequency Fourier coefficients are estimated using second-order empirical moments (subtracting noise variance).
   • SO(3): The degree-zero coefficients are estimated by averaging. The Gram matrices of the coefficient matrices for each angular degree $\ell$ are estimated via second-order moments. These Gram matrices determine the coefficient matrices up to a left-acting unitary matrix (analogous to phase in SO(2)).

2. Likelihood Expansion and Sequential Update:

   • Instead of optimizing the full likelihood, the method optimizes a restricted log-likelihood with respect to a single unresolved degree of freedom (a Fourier phase in SO(2) or a unitary matrix in SO(3)) while holding lower-frequency components fixed.
   • A low-SNR expansion of this restricted likelihood reveals that the leading-order term depends linearly on the unknown variable.
   • SO(2): The update for the phase of the $k$-th coefficient is a closed-form solution derived from a weighted average of the observations, effectively coupling the current frequency to previously recovered lower frequencies. This is shown to be equivalent to a third-order statistic (related to the bispectrum) but derived implicitly.
   • SO(3): The update for the unitary matrix $Q_\ell$ at degree $\ell$ is formulated as a Procrustes problem (maximizing the real part of the trace of a product of matrices). The solution is obtained via Singular Value Decomposition (SVD). The coupling between degrees is governed by Clebsch–Gordan coefficients, which arise naturally from the expansion of the likelihood over SO(3).

3. Algorithm Flow:

   • Initialize lower-frequency components.
   • Sweep from the lowest unknown frequency/degree to the highest.
   • At each step, compute a data-driven average (involving numerical integration over the rotation group for SO(3)) to form a target matrix.
   • Solve the closed-form update (phase rotation or Procrustes alignment) to recover the current component.
   • Higher-frequency components are initially set to zero during the forward sweep and can be refined in subsequent passes if necessary.

Key Contributions

1. Closed-Form Approximate MLE: The paper derives explicit, non-iterative estimators for MRA over SO(2) and SO(3) in the low-SNR regime. This avoids the iterative refinement of EM and the explicit construction/inversion of large invariant tensor systems required by bispectrum methods.
2. Statistical Consistency: The authors prove that the proposed estimator is consistent in the low-SNR asymptotic regime where the number of observations $n$ scales as $n = \omega(\sigma^6)$ (where $\sigma$ is the noise standard deviation). This matches the optimal estimation rate for MRA.
3. Computational Efficiency:
   • SO(2): The complexity is $O(n L N_\theta)$, where $N_\theta$ is the grid size for numerical integration. This is comparable to invariant methods and significantly faster than EM in low SNR, as it avoids the iteration count $T_{EM}$ which explodes as SNR drops.
   • SO(3): The dominant cost per sweep is $O(n(R + N_\Omega)L^3)$. This is more efficient than the bispectrum baseline, which scales as $O(n R^3 L^4)$, particularly because the proposed method avoids the cubic dependence on the number of radial degrees of freedom ($R$) inherent in explicit bispectrum tensor construction.
4. Initialization for EM: The method serves as a robust initialization for the standard EM algorithm. By providing a high-quality starting point, it mitigates the risk of EM getting trapped in local minima and accelerates convergence.

Results

Numerical experiments were conducted on synthetic data generated from molecular volumes (for SO(3)) and bandlimited signals (for SO(2)).

• Accuracy vs. Noise: In SO(2) experiments, as the noise level $\sigma$ increases, the reconstruction error of the proposed fast MLE converges toward that of the full EM algorithm, validating the low-SNR approximation.
• Runtime: The proposed method maintains an essentially constant runtime across varying noise levels. In contrast, EM runtime increases drastically as SNR decreases due to the growing number of required iterations.
• SO(3) Performance: On five molecular volumes (including TRPV1 and 80S ribosome) at SNR levels of $1/50$ and $1/500$:
  • The fast MLE produces meaningful low-frequency reconstructions in under an hour, significantly faster than both the bispectrum method and EM.
  • While the bispectrum method and EM (initialized with the fast MLE) achieve slightly higher Fourier Shell Correlation (FSC) values, the fast MLE provides a favorable accuracy-runtime trade-off.
  • EM initialized with random guesses performs poorly in low SNR, highlighting the necessity of the proposed initialization.

Significance and Claims

The paper positions this work as a foundational step establishing a concrete algorithmic framework based on low-SNR likelihood expansion. The authors claim that:

• The method recovers the same leading-order third-order rotational invariants as the bispectrum but does so implicitly through the likelihood expansion, avoiding the explicit construction and inversion of large bispectrum systems.
• It offers a computationally efficient alternative to EM that is statistically consistent and robust in the challenging low-SNR regimes typical of cryo-electron microscopy (cryo-EM) and cryo-electron tomography (cryo-ET).
• The approach is particularly valuable as an effective initialization for likelihood-based refinement pipelines, bridging the gap between fast invariant methods and statistically optimal but slow iterative methods.

The authors note that while the current work focuses on SO(2) and SO(3), the principles may extend to broader classes of compact group actions, and future work could apply this framework to the full cryo-EM model involving projection operators.