{'fraction_non_alphanumeric': 0.047524934503837844, 'fraction_numerical': 0.020751941903755112, 'mean_word_length': 4.17912401809093, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 36, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Autocorrelation analysis for cryo-EM with sparsity constraints: Improved sample complexity and projection-based algorithms

The number of noisy images required for molecular reconstruction in single-particle cryo-electron microscopy (cryo-EM) is governed by the autocorrelations of the observed, randomly-oriented, noisy projection images. In this work, we consider the effect of imposing sparsity priors on the molecule. We use techniques from signal processing, optimization, and applied algebraic geometry to obtain new theoretical and computational contributions for this challenging non-linear inverse problem with sparsity constraints. We prove that molecular structures modeled as sums of idealized point masses are uniquely determined by the second-order autocorrelation of their projection images, implying that the sample complexity is proportional to the square of the variance of the noise. This theory improves upon the non-sparse case, where the third-order autocorrelation is required for uniformly-oriented particle images and the sample complexity scales with the cube of the noise variance. Furthermore, we build a computational framework to reconstruct molecular structures which are sparse in the wavelet basis. This method combines the sparse representation for the molecule with projection-based techniques used for phase retrieval in X-ray crystallography.


Introduction
Sparsity is a ubiquitous prior in many linear inverse problems, including regression [85,41], compressed sensing [25,17,28], and various image processing applications [27], to name a few. While sparse priors are also used for non-linear inverse problems, their applicability and theoretical foundations are limited to a few specific (usually linear and quadratic) models, e.g., [15,21,31,88]. Motivated by single-particle cryo-EM-an imaging technology for determining the 3-D structure of biological molecules-this paper uses modern techniques from signal processing, optimization, and applied algebraic geometry to provide new theoretical analysis and computational methods for a challenging non-linear inverse problem with sparsity constraints.
Cryo-EM has garnered increasing interest in the past decade due to a series of technological and algorithmic breakthroughs, driving a striking improvement in the obtainable resolution, up to the level where individual atoms can be distinguished. This has in turn opened new scientific horizons and led to many biological discoveries, e.g., [55,3,16].
In a cryo-EM experiment, a solution containing molecules to be imaged is rapidly frozen into a thin ice layer, which is then placed in an electron microscope. Next, the microscope acquires an image, called micrograph, which contains multiple 2-D tomographic projection images of the molecules. The 3-D orientations of individual projection images are unknown and random. To avoid damaging the samples, the electron dose must be kept low, resulting in a low signal-to-noise ratio (SNR). The cryo-EM computational problem is reconstructing the 3-D molecular structure from these projection images [36].
The renewed interest in cryo-EM led to a thorough investigation of its mathematical and statistical foundations [80,6]. In particular, a crucial challenge from a statistical perspective is understanding the sample complexity of cryo-EM, i.e., the number of images that are required to obtain accurate reconstructions. A remarkable result revealed an intimate connection between the sample complexity of cryo-EM (and related statistical models) and the method of moments in the low SNR regime. If the distribution of the 3-D rotations is uniform, the method of moments reduces to autocorrelation analysis. In particular, it was shown that if d is the lowest degree moment of the observations (i.e., the randomly oriented tomographic projections) that determines the molecular structure uniquely, a necessary condition for recovery is n = ω(σ 2d ), where σ 2 is the variance of the noise [67,1]. Specifically, if the distribution of rotations is uniform, then the third-order autocorrelation is the lowest order autocorrelation that determines a generic 3-D structure, implying a sample complexity of n = ω(σ 6 ) [4]. This agrees with long-standing empirical evidence [79].
Autocorrelation analysis was first introduced to cryo-EM by Zvi Kam more than 40 years ago [50]. Kam showed that the second-order autocorrelation of the projection images (which can be estimated with n = ω(σ 4 ) observations) determines the 3-D structure up to a set of orthogonal matrices, under the assumption that the rotations are drawn from a uniform distribution. A few methods have been proposed to resolve the missing orthogonal matrices based on typically unavailable side information, such as homologous models with known structure, or a few clean projections images, in order to construct ab initio models [12,59,46]. These ab initio models can then be refined using expectation-maximization: the prevalent algorithmic framework for the cryo-EM reconstruction problem that aims to maximize the non-convex posterior distribution [79,74,69]. In addition, it was recently shown that if the distribution of rotations is non-uniform, then there is at most a finite list of structures that agree with the second moment of the observations [76]. Techniques that are inspired by Kam's method were also proposed as a solution to the molecular reconstruction problem in X-ray free-electron lasers (XFEL), which, akin to the reconstruction problem in cryo-EM, involves recovering a 3-D structure from its randomly oriented diffraction images [71,56,24,57]. In contrast to cryo-EM, in XFEL the rotations are always uniformly distributed and the reconstruction problem is more involved since the measurements consist of the magnitudes of Fourier coefficients without their phases.
The first contribution of this paper is proving that if the sought-for molecular structure can be described as a sparse combination of ideal point masses, then the structure can be determined uniquely from the second-order autocorrelation of the observations, even if the rotations are distributed uniformly. This eliminates the orthogonal matrix ambiguity in Kam's original paper [50]. This result is the first theoretical guarantee for unique recovery from the second moment. It combines the sparsity assumption with proof techniques from real algebraic geometry, substantially reducing the sample complexity of the cryo-EM reconstruction problem from n = ω(σ 6 ) to n = ω(σ 4 ). The argument is constructive in the sense that it provides a polynomial-time recovery algorithm. However, said algorithm is tailored to the specific model of point mass functions. It is not well-suited to data discretized into pixels or voxels because it hinges on the ability to cluster points in the support of the second moment into certain distinct components, see Remark 11. (Accurate clustering becomes difficult when data is discretized and the components are close to each other.) The second contribution of this paper is a practical algorithm, fusing the secondorder autocorrelation of the projections with a sparse representation of the molecular structure, and requiring only n = ω(σ 4 ) observations. The algorithm builds on the realization that a typical 3-D structure can be represented by only a few coefficients in a suitable basis. Similar sparsity assumptions have been leveraged in a wide variety of scientific and engineering applications, including compressed sensing [25,17,28], image processing [27], and phase retrieval [15,21,47,31,88], to name a few. In cryo-EM, there have been several attempts to represent the 3-D structure as either a sparse mixture of Gaussians [87,49,48,52,89,19,70] or using alternative bases [86]. Yet, the sparsity property has still not been fully harnessed to represent and recover 3-D molecular structures, and it is not part of the standard computational pipeline of cryo-EM.
The technique is based on a new connection between Kam's theory for cryo-EM and the crystallographic phase retrieval problem-recovering a sparse signal from its Fourier magnitudes. In particular, we adapt projection-based algorithms that were designed for the crystallographic phase retrieval problem to the cryo-EM setting. These algorithms were extensively validated on experimental X-ray crystallography datasets by prior researchers, see for example [34,29,60,32,31]. Here, we demonstrate on simulated data that they are also useful in constructing ab initio models in cryo-EM. They can be used to mitigate computational and model bias issues associated with the non-convexity of the cryo-EM reconstruction problem [82,69,43]. This new computational approach opens the door to merging more aspects of the phase retrieval and cryo-EM fields in future work.
The rest of the paper is organized as follows. In Section 2, we provide background on the reconstruction problem in cryo-EM, the method of moments, Kam's theory, and the crystallographic phase retrieval problem. In Section 3, we prove that a structure composed of an ensemble of ideal point masses subject to uniform rotations can be recovered from the second order autocorrelation, implying a sample complexity of n = ω(σ 4 ). Section 4 outlines the practical computational framework and presents numerical results. Section 5 concludes the paper and discusses potential theoretical and computational extensions.

Preliminaries
2.1. The cryo-EM problem. Cryo-EM reconstruction seeks to determine a 3-D molecular structure Φ from its 2-D noisy tomographic projections, taken at random viewing angles. In this work, we focus on the case of uniformly random rotations. Uniformity is often taken as a baseline model, was the setting in Kam's paper [50], and is harder than the case of a non-uniform distribution of rotations in the sense that it requires asymptotically more images (without sparsity priors) [4,76]. Formally, let µ be the Haar probability measure on the compact group SO(3) of 3-D rotations, representing the uniform distribution. Assuming that we observe i.i.d. 2-D images of Φ after it has been randomly rotated according to µ and then tomographically projected to the plane, each projection image is modeled as: where ε(x, y) is a white Gaussian noise with known variance σ 2 , and R · Φ denotes the action of the rotation R on Φ. Here, typically, the variance of the noise σ 2 is much greater than the magnitude of the clean projection. The cryo-EM problem is to estimate the molecular structure Φ from n realizations of (1), i.e., from the 2-D observations I R1 , I R2 , . . . , I Rn . In Section 5, we discuss how the proposed framework can be extended to account for additional aspects in the generative model for cryo-EM images.
2.2. The method of moments. Our theoretical and computational contributions are based on the method of moments-a basic statistical inference technique tracing back to the seminal paper of Karl Pearson in the end of the 19th century. Specifically, we use the second moment of the observations (1), and relate it to the sought-for 3-D structure.
The (debiased) second observable moment is given by where B(σ 2 ) is a bias term that depends only on the noise variance. For large enough n, we have M 2 ((x 1 , y 1 ), (x 2 , y 2 )) ≈ M 2 ((x 1 , y 1 ), (x 2 , y 2 )), where (4) M 2 ((x 1 , y 1 ), (x 2 , y 2 )) := denotes the (debiased) population second moment, which is a function of Φ through (1). More precisely, for n = ω(σ 4 ) it holds M 2 = M 2 + o(1) with high probability. The idea of the method of moments is to find a structure Φ which matches the observable moments. It is an alternative to other standard statistical estimation methods, e.g., maximum likelihood estimation (MLE). That said, a recent paper suggests that in the low SNR regime, the method of moments approximates the MLE [51].

2.3.
Kam's method. We detail a specific approach to autocorrelation analysis in cryo-EM, introduced by Kam [50]. To this end, we need to introduce a convenient basis for representing a 3-D structure Φ, the spherical Bessel basis [2]. An expansion of maximum degree L is defined by first expanding the Fourier transform F(Φ) of Φ in spherical harmonics as where k denotes the radial frequency and Y m (θ, ϕ) are the spherical harmonics basis functions. In addition, the spherical harmonics coefficients A m (k) are expanded by spherical Bessel functions, up to degree S , as The functions j s (k) are the normalized spherical Bessel functions. By allowing L and S to grow unboundedly, any structure of a biological molecule can be represented in this basis. However, when expanding 3-D molecular structures from discretized projection images, the Nyquist criterion applied to the projection images limits the amount of extractable information. This determines bounds on the maximally allowable truncation parameters L and S ; see Appendix C for a detailed description. We aim to recover the coefficients a ms (up to rotation and reflection in R 3 ). As will be shown next, it is convenient to gather the coefficients into matrices A , of size S × (2 + 1), via A (s, m) := a ms , for each = 0, . . . , L.
Provided that the distribution of viewing angles is uniform, Kam [50] showed that the second moment of the Fourier transform of the projection images provides estimates for the following S × S matrices: Applying the Cholesky decomposition to each C in (7) and imposing Φ to be realvalued, knowledge of (7) identifies each matrix of coefficient A up to an unknown real, orthogonal transformation, provided S ≥ 2 + 1. That is, we can compute A O for some unknown orthogonal matrix O in the group O(2 + 1). Therefore, the second moment determines Φ up to a set of orthogonal matrices. To recover these matrices, and thus the 3-D structure, additional information is required. In this paper, we suggest using a sparsity assumption.
2.4. Crystallographic phase retrieval. One of the contributions of this paper is to relate the cryo-EM reconstruction problem to crystallographic phase retrieval. Phase retrieval is the main computational challenge in X-ray crystallography, which is still a leading method for elucidating the atomic structure of molecules. The prevalence of crystallography is witnessed by the remarkable fact that 25 Nobel Prizes have been awarded for work directly or indirectly involving crystallography [38]. Although there exist additional important phase retrieval applications (see for example [77,24,7,5]), X-ray crystallography is by far the most widely investigated application. The crystallographic phase retrieval problem entails recovering a sparse signal x from its periodic autocorrelation (or, equivalently, from its Fourier transform magnitudes, namely, its power spectrum). While simply stated, and despite its importance, the theoretical foundations of this problem have been investigated only recently. In particular, it was conjectured that a generic sparse signal can be recovered from its periodic autocorrelation if the number of non-zero entries is smaller than half the signal's length [9]. This conjecture was verified for a few cases. The relation of the crystallographic phase retrieval problem with the beltway problem from combinatorial optimization is explored in [39]. Our theoretical reconstruction guarantees in the following section can be viewed as analogous results in the setting of cryo-EM.
The standard algorithms for crystallographic phase retrieval build on two projection operators: one onto the measured data (the power spectrum) and the second onto the space of sparse signals. While simple algorithms that alternate between these two projections tend to quickly stagnate, a more sophisticated family of algorithms, based on reflections, shows excellent performance, though their running time is exponential in the sparsity level [31,11]. These algorithms are tightly related to splitting methods, such as Douglas-Rachford and the alternating direction method of multipliers (ADMM), and have been applied to a wide variety of problems [32]. A main contribution of this paper is a modification of these algorithms to autocorrelation analysis for cryo-EM. In particular, we focus on one such algorithm, called relaxed-reflect-reflect (RRR), but alternative algorithms, such as Fienup's hybrid input-output algorithm [34], the difference map algorithm [29], and the relaxed averaged alternating reflections algorithm [60], can be adapted to cryo-EM by the same strategy. Importantly, if the model is accurate (e.g., no noise) RRR iterations halt only when they find a solution that satisfies both constraints (defined by the projection operators). Thus, RRR does not suffer from local minima as gradient-based algorithms do.

Superior sample complexity: The second moment suffices for sums of point masses
This section presents our main theoretical result: the second moment suffices to recover an idealized sparse volume, i.e., a volume given as a weighted sum of point masses. We deduce that the second moment also suffices for a pixelated and blurred variant of the model. Our theorems imply an associated sample complexity of n = ω(σ 4 ). This stands in contrast to previous results which do not assume sparsity. There, the third moment is required for recovery, and the associated sample complexity is n = ω(σ 6 ) [4, 33].
3.1. Models and main theoretical results. We use an atomistic representation of a molecule. In our first idealized model, an atom is specified by a weighted Dirac delta function, and a molecule is a sum of such point masses. Let a 1 , . . . , a p ∈ R 3 be the 3-D points representing atom locations, and w 1 , . . . , w p be positive weights corresponding to the scattering potentials of the individual atoms. Then is the molecule composed of the atoms (a i , w i ). Relabeling if necessary, we assume that the 2 -norms a i are in descending order. We model each projection image I R as a mixture of p Dirac delta functions on R 2 plus noise: Here π : R 3 → R 2 denotes coordinate projection onto the first two coordinates, and ε is a white Gaussian noise with (known) variance σ 2 . Given n 2-D images as in (9), the reconstruction problem is to recover the atoms (a i , w i ) in (8) up to a global rotation and reflection. (The reflection ambiguity exists because a molecule and its reflection in the microscope's image plane are indistinguishable given cryo-EM data, see e.g. [36].) See Figure 1 for an illustration of the setup. Under this model the (debiased) second population moment, obtained by substituting (9) into (4), reads: where µ is the uniform distribution. Note that M 2 is a measure on R 2 × R 2 .
We introduce certain assumptions on the atom locations: A1. The vectors a i are pairwise linearly independent; A2. The norms a i are distinct.
We remark that conditions A1 and A2 are quite restrictive, e.g., ruling out molecules with nontrivial point-group symmetries.
Our first main theoretical result is stated as follows.
Theorem 1. Consider the model given by (8)- (9). Assume that conditions A1-A2 hold. Then the support of the second moment M 2 uniquely determines the set . . , p}. Therefore, M 2 (Eq. (10)) uniquely determines the set of atom locations {a i : i = 1, . . . , p} up to a rotation and reflection in R 3 .
Theorem 1 is proven in Subsection 3.3, after auxiliary results are given in Subsection 3.2.
Building on the uniqueness in Theorem 1, we obtain the following constructive result as well.
Theorem 3. Consider the model given by (11). Fix an integer p ≥ 3 and real numbers r > 0 and w + > w − > 0. Assume that A1-A2 hold, and for each i = 1, . . . , p we have a i ≤ r and w − ≤ w i ≤ w + . Then there exists m = m (p, r, w + , w − ) with the following property. Whenever m ≥ m and 2 m × 2 m pixels are used in (11), then the second moment M As the details are technical, we prove Theorem 3 in the appendices. We only use two properties of the Gaussian kernel k: that it is real-analytic and that its Fourier transform does not vanish. Corollary 4. Assume the setting of Theorem 3 with m ≥ m . Then, the sample complexity for generic unique recovery (in the sense of [4]) is n = ω(σ 4 ) as σ → ∞.
The rest of the section provides the proofs of Theorems 1 and 2, with Theorem 3, Corollary 4 and supporting results shown in the appendices. We emphasize that Algorithm 1 is a theoretical algorithm, not intended for use in practice due to its noise sensitivity as explained in Remark 11. By contrast, Algorithm 2 in the subsequent section is built for practical situations.
3.2. Support of M 2 . To recover the atoms from M 2 , the main information that we use is actually qualitative. Specifically, we rely on the particular structure of the support of the second moment M 2 in R 2 × R 2 . To describe this, we need to first understand the possible images of one pair of atoms.
The next lemma characterizes S ij ⊆ R 2 × R 2 as the solution set to a system of polynomial equations and inequalities. This will enable proof techniques from real algebraic geometry.
Lemma 7. S ij is connected, compact and semialgebraic. Letting ((x 1 , y 1 ), (x 2 , y 2 )) be variables on R 2 ×R 2 , then S ij is cut out by one quartic equation and two quadratic inequalities: It has dimension 3 as a semialgebraic set if condition A1 holds.
A few different examples of the sets S ij are illustrated in Figure 2, for varying values of a i , a j and a i , a j . There we show the projection of S ij to R 3 when y 2 is dropped.
The next result is immediate from Definitions 5 and 6.
3.3. Information-theoretic uniqueness: Proof of Theorem 1. We begin by proving Theorem 1. The key is a converse to Lemma 7. While Lemma 7 implies the quartic equation in (13) (plus the quadratic inequalities there) determine the set S ij , we need that S ij determines the quartic. The proof of this converse uses results from real algebraic geometry [14]. Figure 2. Illustration of the semialgebraic sets S ij described in Lemma 7, projected onto (x 1 , y 1 , x 2 )-space. Upper row from left to right: a i = a j = 1 and a i , a j = cos π 2 , cos π 4 , cos π 6 , cos π 20 , respectively. Lower row from left to right: a i , a j = a j cos π 3 , a i = 1 and a j = 1, 0.9, 0.8, 0.7, respectively. By Lemma 13, these sets make up the support of M 2 , which is the basis of our theory.
Lemma 9. Assume that condition A1 holds. Let i = j. Then the ideal of the real Zariski closure of S ij in R 2 ×R 2 is principal and generated by the quartic polynomial Further, q ij is irreducible over R.
Corollary 10. Assume that conditions A1-A2 hold. Then, the irredundant irreducible decomposition of the Zariski closure of the support of M 2 is We can now prove the information-theoretic uniqueness.
Proof of Theorem 1. The support of M 2 determines the real radical prime ideal of each of top-dimensional irreducible component of its Zariski closure. By A1-A2, Corollary 10 and Lemma 9, these ideals are q ij = {q ij h : h ∈ R[x 1 , y 1 , x 2 , y 2 ]} for i = j. The ideal q ij uniquely determines q ij , since the coefficient of x 2 1 y 2 2 in (15) is 1. Extracting the coefficients of x 2 2 , x 2 1 , x 1 x 2 in (15), q ij determines the triple ( a i 2 , a j 2 , a i , a j ). Ranging over i, j, we have proven that the support of M 2 fixes the set: (17) {( a i 2 , a j 2 , a i , a j ) : i = j}.
By A2 and our assumption the norms a i are descending, knowledge of (17) lets us fill in the Gram matrix: However G determines A up to left multiplication by a 3 × 3 orthogonal matrix. Indeed considering a truncated rank-3 eigendecomposition, we write where Q ∈ R p×3 has orthonormal columns and D ∈ R 3×3 is diagonal and positivesemidefinite. Then, . Therefore, the atoms' locations a i are determined up to a global rotation and reflection.

3.4.
Recovery algorithm: Proof of Theorem 2. Now, we move forward and prove Theorem 2. We present Algorithm 1 for efficiently recovering the atoms (a i , w i ) from M 2 . The algorithm is theoretical in that it relies on oracle access to the following information.
Assumption 1. We assume oracle access to: consists of four or more Zariskigeneric points on S ij ; O2. the value of the measure M 2 on the set S ij for all i = j.
Remark 11. In principle, O1 and O2 can be estimated from the sample moment M 2 (2) if n = ω(σ 4 ). It would require the ability to identify points in the support of M 2 and cluster them according to the components S ij . However this encounters difficulty when dealing with noisy moments discretized in pixels. The next section is dedicated to a different computational framework, better-suited for practical settings.
Proceeding, Algorithm 1 interpolates sample(S ij ) to recover q ij from Lemma 12.
Lemma 12. Assume that condition A1 holds, and O1 is known. Let sample(S ij ) = {((x 1k , y 1k ), (x 2k , y 2k )) : k = 1, . . . , |sample(S ij )|}. Consider the matrix Then it has rank 4, with kernel spanned by By Lemma 12, we compute the triples {( a i 2 , a j 2 , a i , a j ) : i = j} in O(p 2 ) time, by forming the matrices (22) and computing their kernels. We then fill in the Gram matrix (18) as in the proof of Theorem 1. The atoms' locations a i are recovered from the truncated eigendecomposition (21).
The calculation of the weights w i is based on the following.
Lemma 13. Assume that conditions A1-A2 hold. Then, for each i = j, the measure of S ij with respect to M 2 is Therefore, O2 tells us all off-diagonal entries of ww ∈ R p×p . We complete this uniquely to a rank-1 matrix by using where j, j are any indices such that i, j, j are all distinct. (This step requires p ≥ 3.) The weights w i are lastly recovered either by computing the leading eigenvector/eigenvalue pair of ww or as the square root of the diagonal of ww , using the fact that the w are non-negative.

Remark 14.
We note that (25) is a particular case of the problem of recovering a low-rank matrix with corrupted diagonal entries; see, e.g., [65,72,73] for more on that problem.
We summarize the procedure of this section in Algorithm 1.

Kam's method with sparsity constraints
This section introduces a computational framework to leverage sparsity in recovering the underlying molecular structure. The goal is to devise a principled way to compute ab initio approximations of the underlying structures, that can then be improved further in a refinement step (which is typically performed using expectation-maximization [74,69]). In this section, we maintain the assumption that the distribution µ of viewing angles is uniformly distributed. If µ is a nonuniform distribution, it is known that there is at most a finite list of structures that are consistent with the observed second order moment [76]; employing sparsity to aid in the recovery problem with non-uniform distribution will be considered in future work.
We use projection-based optimization techniques from the related problem of crystallographic phase retrieval, coupled with information extracted from the second moment of the projection images. Without imposing the underlying sparsity, the second moment of the projection images determines the structure up to an ambiguity encoded by a set of unknown orthogonal matrices. The key idea of the algorithm is to alternatingly project the molecular structure onto constraints encoded by the sparsity and by the projection image moments, respectively.
Analogously to (10), the moments of the projection images furnish information about the underlying 3-D structure. Unlike our theoretical results, however, we consider a general 3-D structure Φ expanded in a spherical Bessel basis as in Section 2.3, and assume it can be represented by only a few wavelet coefficients. The next section introduces wavelet bases, and later we provide the details of the projection-based algorithm.

Wavelet bases.
We encode sparsity of a 3-D molecular structure Φ by a sparse expansion in wavelets [22]-a popular choice of sparsifying, localized bases in a wide range of applications [63]. Our algorithm can easily be adapted to any specific wavelet basis, and, more generally, to any choice of basis, for instance, sparsifying bases learned through data.
We denote the multilevel wavelet basis by f m,n , where m = 1, . . . , m max denotes the level of the wavelet and n = 1, . . . , n max (m) the index of the function within the level. As a shorthand, we define W : R M ×M ×M → R mmaxnmax(m) as the map sending a 3-D structure to its vector of coefficients when expanded in the wavelet basis, i.e., The inverse mapping SB −1 : L =0 C S ×(2 +1) → R M ×M ×M then expands a set of coefficients in the spherical Bessel basis into its corresponding 3-D structure: As discussed in Section 2, Kam's method identifies matrices A O , for O an unknown orthogonal matrix, for each with S ≥ 2 + 1. By possibly reducing the value of L to the largest index with this property, we will for ease of notation assume that this property holds for = 0, . . . , L. Therefore, at the onset of the algorithm, we have access to a set of coefficient matrices B = (B 0 , . . . , B L ) satisfying for unknown orthogonal matrices O . Our algorithm aims to recover an approximation of these unknown orthogonal matrices, which leads to an approximation of Φ. This orthogonal matrix retrieval problem is an analogue to the problem of the missing phases in the phase retrieval problem [10]. We therefore adapt a popular algorithm from the phase retrieval literature into the problem of cryo-EM. The algorithm repeatedly utilizes two projections onto the set of structures with a given sparsity level and a set determined by the projection images. These two projections are the main pillars of the algorithm and can be used in different ways, as explained next. But first, we define the two projection operators.
4.2.1. First projection: Moment constraint. We begin by defining the first projection operator, denoted by ρ 1 , as the projection onto the set defined by the C matrices in (7). Let SB(Φ) = (A 0 , . . . , A L ) in L =0 C S ×(2 +1) be the ordered collection of matrices of coefficients in the spherical Bessel basis. Define ρ 1 (Φ) as the projection where the matrices D are defined by with C from (7). (32) is an instance of the Orthogonal Procrustes problem. Although it is a non-convex optimization problem, it can be solved in closed form in terms of the singular value decomposition of B T A , see e.g., [75]. In the implementation, the matrices defining the operations SB and SB −1 are precomputed. The computational complexity of subsequently solving an instance of (31) is then

4.2.2.
Second projection: Sparsity constraint. The projection ρ 2 promotes sparsity in a given local wavelet basis. For a structure Φ and an integer K, define ρ 2 (Φ, K) as the structure with wavelet coefficients obtained by retaining the K largest components of W (Φ) and replacing the remaining elements by zero, i.e., and α m,n = 1 if c m,n has magnitude among the K largest magnitudes of the c m,n , and zero otherwise. The computational complexity of this step is O(M 3 ). We again emphasize that, generally, any localized basis or frame can be used to define ρ 2 , and we fix a wavelet basis for the sake of definiteness.

Algorithm.
A straightforward algorithm to attempt to recover Φ is through alternating projections. This procedure is described by fixing a sparsity level K and iterating the two projections ρ 1 and ρ 2 in turn. The use of the two projections in an alternating fashion is intended to promote convergence to an intersection point of the two sets. In the case of projecting onto convex sets, convergence results are known [20], but convergence is not guaranteed for the non-convex projections in (31) and (33). Indeed, for non-convex sets, alternating projection schemes frequently suffer from convergence to local minima, and a method to escape the local minima is required. To achieve this, the phase retrieval literature details different iteration schemes combining the two projections ρ 1 and ρ 2 in different ways, for instance using the Relaxed-Reflect-Reflect (RRR) algorithm [30,31]. In terms of the projection operators, this iterative scheme can be written out as where β ∈ (0, 2) is a scalar hyperparameter. The algorithm is summarized in Algorithm 2. As aforementioned, other phase retrieval algorithms which are based on two projection operators, such as the difference map algorithm and the relaxed averaged alternating reflections algorithm, can be adapted to cryo-EM similarly.

Algorithm 2 Recovering Φ
Input: Projection images, sparsity level K, maximum number of iterations N , hyperparameter β. Output: Estimated structure Φ (1) Form the matrices C from (7) (2) Compute the Cholesky factorizations of the C to produce matrices B in (30)

4.3.
Simulation results. We apply Algorithm 2 to two structures from the online EM data bank [58], EMD-0409 [44] and EMD-25892 [61], as well as the Shepp-Logan phantom [78]. For each structure, we run Algorithm 2 for a given value of L and K to obtain a reconstruction. To measure the reconstruction quality, we follow the standard procedure in the cryo-EM community, and compute the Fourier shell correlation (FSC) between the estimated structure and the ground truth. Specifically, the FSC of two structures Φ 1 and Φ 2 is defined by where one structure is the estimated structure, and the second is the ground truth. The FSC is real-valued because of symmetry of the summation. The resolution is determined when the FSC curve drops below 0.5. EMD-0409 has dimensions 128×128×128, with each voxel having physical length of 1.117Å. EMD-25892 has dimensions 320 × 320 × 320, and voxel size 1.68Å. The volumes were downsampled by a factor of 2 and 5, respectively, to give structures of size 64 × 64 × 64. The ground truth matrices were generated exactly and the matrices B in (30) were generated using O chosen uniformly at random. To fix the units for the Shepp-Logan phantom, we assume the voxels to have side length 1Å. The simulations used Haar wavelets to define W .
The result of applying Algorithm 2 to each structure is shown in Figure 3. For all three example structures, during a run of the algorithm, the resolution initially rapidly improves. Afterwards, the improvement slows down and exhibits an exploratory and oscillating behavior. This is typical for RRR-based algorithms [31].
During the run of the algorithm, the optimal resolutions obtained for the three structures were 17.2Å, 20.52Å, and 4.59Å, respectively, and the resolutions at the initialization of the algorithm were 21.3Å, 329.8Å, and 8.4Å, respectively. As a comparison, the resolutions between the ground truth structures and their truncation into the spherical Bessel bases with the chosen values of L are 6.0Å, 18.50Å and 2.20Å, respectively; these resolutions are bounds on the optimally obtainable resolutions. Figure 4 also shows a comparison of EMD-0409 with its reconstruction, truncated to different values of L. Visually, the reconstructed element captures the relevant features of the ground truth for each value of L. The resolution increases with L.
Movie 1 visualizes the reconstructed volume as a function of the iteration number. Note that the reconstruction at each step is visually similar to the ground truth, although the computed resolution noticeably improves during the run of the algorithm. This implies that even knowledge of the coefficients A O with the wrong rotation matrices O provides some information about the ground truth.

Discussion
The contribution of this paper is twofold. As the first contribution, our theoretical results imply that a sparse mixture of point masses can be uniquely recovered from the second order moment, even in the case of a uniform distribution of viewing angles, whereas previous work has only proven recovery using the third order moment. Thus, fewer images are required for reconstruction. This has a number of potential experimental implications. Firstly, since microscope time is expensive, this may greatly reduce the cost of the experimental part of the cryo-EM pipeline. This might be especially important for XFEL, where throughput is a major bottleneck and viewing directions are uniformly distributed [84,54]. Secondly, it may enable reconstruction of structures where a limited number of projection images can be captured. This might be the case, for example, when the molecule may appear in several conformational states, and a limited number of images will be available for each conformation.
The second contribution is a new algorithm for ab initio modeling. The computational framework introduced in this article opens the door to incorporating a number of promising techniques from crystallographic phase retrieval into cryo-EM algorithms. There is, for instance, flexibility in choosing the projection operators ρ 1 and ρ 2 . It may include biologically-oriented priors, such as minimum atom-atom distance or Wilson statistics [81,40], or data-driven priors based on previously resolved structures [53]. A systematic study of adapting these techniques will be initiated in coming work. Additional future work includes extending the use of sparsifying priors in other parts of the cryo-EM reconstruction pipeline, for instance in existing approaches to iterative refinement [74] or in autocorrelation analysis using micrographs without particle picking [8]. Yet another important direction is to incorporate the sparsity prior into reconstruction by the method of moments when there is a non-uniform distribution of viewing directions [76].
Appendix A. Proofs of auxiliary results for Theorems 1 and 2 A.1. Proof of Lemma 7. Note that S ij is connected and compact, since S ij = θ ij (SO(3)) and SO (3) is compact and connected, while θ ij is continuous. It also is semialgebraic, as SO (3) is a real algebraic variety and θ ij is a polynomial map (see the Tarski-Seidenberg theorem [14]).
For the converse, take ((x 1 , y 1 ), (x 2 , y 2 )) ∈ T ij . Let where ε = sign ( a i , a j − x 1 x 2 − y 1 y 2 ). Put b i = (x 1 , y 1 , z 1 ) and b j = (x 2 , y 2 , z 2 ) in R 3 . By the choice of z 1 and z 2 , Also, from the equality constraint in (13), it holds z 2 This with the choice of ε implies From (36) and (37), there exists R ∈ SO(3) such b i = Ra i and b j = Ra j . Hence ((x 1 , y 1 ), (x 2 , y 2 )) ∈ S ij , whence T ij ⊆ S ij . We conclude T ij = S ij . The dimension of S ij as a semialgebraic set is the maximal dimension of a cell in any cylindrical algebraic decomposition of it [14,Cor. 2.8.9]. This agrees with the maximal rank attained by the differential of θ ij : where T denotes tangent space. We recall that the tangent space to rotation matrices is parameterized by skew-symmetric matrices. Specifically, T R (SO (3) Putting (x 1 , y 1 , z 1 ) := Ra i and (x 2 , y 2 , z 2 ) := Ra j , we rewrite Thus, dim(S ij ) = 3, unless W (R) is rank-deficient for all R ∈ SO(3). We claim it is rank-deficient for specific R if and only if (x 1 , y 1 , z 1 ) and (x 2 , y 2 , z 2 ) are linearly dependent or z 1 = z 2 = 0. This is proven using a computer algebra system, e.g. [42]. Indeed, if I is ideal in the ring Q[x 1 , y 1 , z 1 , x 2 , y 2 , z 2 ] generated by the 3 × 3 minors of W (R), the claim follows from calculating the primary decomposition [26]: Given the claim, rank(W (R)) < 3 for all R ∈ SO(3) if and only if a i and a j are linearly dependent. In other words, S ij has dimension 3 if and only if a i and a j are linearly independent. The proof of Lemma 7 is complete.
A.2. Proof of Lemma 8. As mentioned in the main text, this is immediate from Definitions 5 and 6 and (10).
A.3. Proof of Lemma 9. If X is a subset of a real Euclidean space R k , we write X for the Zariski closure in R k , and I(X) for the real radical ideal of X.
First, note that S ij is irreducible. This is because S ij = θ ij (SO (3)), θ ij is polynomial and SO (3) is an irreducible algebraic variety. So, I(S ij ) is prime [14,Thm. 2.8.3(ii)]. Also, S ij has dimension 3 as an algebraic variety. This is by [14,Prop. 2.8.2], and Lemma 7 which states S ij has dimension 3 as a semialgebraic set. So, I(S ij ) has height 1 [14, Def. 2.8.1]. Here, every prime ideal with height 1 is principal, as R[x 1 , y 1 , x 2 , y 2 ] is a unique factorization domain. It follows that (38) I(S ij ) = f , for some irreducible polynomial f ∈ R[x 1 , y 1 , x 2 , y 2 ], where angle brackets indicate ideal generation. By Lemma 7, we know S ij ⊆ Z(q ij ), where Z denotes the zero set in R 2 × R 2 . Taking closures, S ij ⊆ Z(q ij ). Equivalently, q ij ∈ I(S ij ). By (38), this means f evenly divides q ij , say, (39) q ij = f g, for some g ∈ R[x 1 , y 1 , x 2 , y 2 ]. To conclude the proof, it suffices to prove that g is a nonzero scalar. Then q ij = f = I(S ij ), and q ij is irreducible because f is. For a contradiction, assume that f has positive degree. Then, (39) implies where the subscript indicates the top total degree part of the polynomial. Here, Thus, (40), the assumption that f has positive degree and unique factorization, we deduce that (possibly after multiplying by nonzero scalars) for some α, β, γ, δ, ε, ζ, η, θ, ι, κ ∈ R. Now we insert (42) and (41) into (39). Equating the constants and the coefficients of x 2 1 , y 2 1 , x 2 2 , y 2 2 gives All the left-hand sides in (43) are nonzero because a i and a j are linearly independent. Therefore, α, . . . , κ are all nonzero. Next, we equate the coefficients of x 1 , y 1 , x 2 , y 2 in (39). The result is that whence the 4 × 2 matrix has rank 1. So all its 2 × 2 minors vanish, in particular Finally, we equate the coefficients of x 1 y 1 in (39): Eq. (44) and (45) imply But this contradicts the earlier finding that α, . . . , κ are all nonzero. So the assumption that f has positive degree is false, and f is a nonzero scalar. This proves Lemma 9.
A.4. Proof of Corollary 10. By Lemma 8, the support of M 2 is ∪ p i,j=1 S ij . This has Zariski closure ∪ p i,j=1 S ij . We claim its irredundant irreducible decomposition is The claim follows from several facts. First, for all i = j, S ij is irreducible. It has dimension 3 and defining equation q ij (41), by A1 and Lemma 9. When i = j, i = j and (i, j) = (i , j ), then q ij and q i j are not scalar multiples of each other by A2 (cf. their coefficients on x 2 1 , x 2 2 , x 2 1 y 2 1 ). Hence S ij = S i j . Next, because when we substitute x 1 = x 2 and y 1 = y 2 into (41) we get a nonzero result as the constant term does not vanish. All together, (47) is the claimed irredundant irreducible decomposition as wanted.
A.6. Proof of Lemma 13. From (14), we have Here µ(θ −1 ij (S ij )) = µ(SO(3)) = 1 by definition of S ij . On the other hand, for all i , j = 1, . . . , p with (i , j ) = (i, j), S ij ∩ S i j is a semialgebraic set with positive codimension in S i j by the fact that (47)  Step 0: We state a general fact about real analytic functions that we will use: Let H(y, z) : R m × R n → R be real analytic jointly in (y, z). Let ν be a compactly supported and absolutely continuous measure on R n . Then H(y, z)dν(z) is realanalytic in y. This can be justified by appropriately differentiating under the integral sign, see [35].
Step 1: We introduce notation. First, define where B(0, r) is the 2 -ball of radius r in R 3 centered at 0. X is the space of possible molecules.
We define I ∀m ≥ m .
Thus the corresponding zero sets in X ×2 stabilize too: Equivalently, for all {(a i , w i )} ∈ X we have: (50) the fiber (M Step 3: We deduce an equality in which there is no pixelation.