Method of moments for 3-D single particle ab initio modeling with non-uniform distribution of viewing angles

TL;DR

This paper extends the method of moments for cryo-EM ab initio modeling to handle unknown, non-uniform viewing angle distributions, simplifying the statistical requirements and proposing algorithms for structure reconstruction.

Contribution

It removes the uniformity assumption in Kam's method, demonstrating that first and second moments suffice for non-uniform distributions and developing algorithms for both known and unknown distributions.

Findings

First and second moments are sufficient for non-uniform distributions.

Efficient algorithms are proposed for known distribution cases.

Non-convex optimization methods can recover structure and distribution in unknown cases.

Abstract

Single-particle reconstruction in cryo-electron microscopy (cryo-EM) is an increasingly popular technique for determining the 3-D structure of a molecule from several noisy 2-D projections images taken at unknown viewing angles. Most reconstruction algorithms require a low-resolution initialization for the 3-D structure, which is the goal of ab initio modeling. Suggested by Zvi Kam in 1980, the method of moments (MoM) offers one approach, wherein low-order statistics of the 2-D images are computed and a 3-D structure is estimated by solving a system of polynomial equations. Unfortunately, Kam's method suffers from restrictive assumptions, most notably that viewing angles should be distributed uniformly. Often unrealistic, uniformity entails the computation of higher-order correlations, as in this case first and second moments fail to determine the 3-D structure. In the present paper, we remove this hypothesis, by permitting an unknown, non-uniform distribution of viewing angles in MoM. Perhaps surprisingly, we show that this case is statistically easier than the uniform case, as now first and second moments generically suffice to determine low-resolution expansions of the molecule. In the idealized setting of a known, non-uniform distribution, we find an efficient provable algorithm inverting first and second moments. For unknown, non-uniform distributions, we use non-convex optimization methods to solve for both the molecule and distribution.

Tables

Table 1. Table 1 . Uniqueness for inverting the first two moments in the case of a known, in-plane uniform ρ 𝜌 \rho , according to band limits. Generically finitely many solutions for A 𝐴 A is denoted by ✓ , infinitely many solutions for A 𝐴 A is denoted by ✗ , and indecisive numerics is denoted by ? .

	$L=2$	$L=3$	$L=4$	$L=5$	$L=6$	$L=7$	$L=8$	$L=9$	$L=10$
$P=0$	✗	✗	✗	✗	✗	✗	✗	✗	✗
$P=1$	✗	✗	?	✓	✓	✓	✓	✓	✓
$P=2$	✗	✓	✓	✓	✓	✓	✓	✓	✓
$P=3$	✗	✓	✓	✓	✓	✓	✓	✓	✓
$P=4$	✗	✓	✓	✓	✓	✓	✓	✓	✓

Table 2. Table 2 . Uniqueness for inverting the first two moments in the case of an unknown, totally non-uniform ρ 𝜌 \rho , according to band limits. Generically finitely many solutions for ( A , B ) ​  mod SO(3) 𝐴 𝐵  mod SO(3) (A,B)\textup{ mod SO(3)} is denoted by ✓ , finitely many solutions for A ​  mod SO(3) 𝐴  mod SO(3) A\textup{ mod SO(3)} but infinitely many solutions for B ​  mod SO(3) 𝐵  mod SO(3) B\textup{ mod SO(3)} is denoted by ∼ ∼ \thicksim , infinitely many solutions for A ​  mod SO(3) 𝐴  mod SO(3) A\textup{ mod SO(3)} is denoted by ✗ , and indecisive numerics is denoted by ? .

	$L=2$	$L=3$	$L=4$	$L=5$	$L=6$	$L=7$	$L=8$	$L=9$	$L=10$
$P=0$	✗	✗	✗	✗	✗	✗	✗	✗	✗
$P=1$	✓	✓	✓	✓	✓	✓	✓	✓	?
$P=2$	✓	✓	✓	✓	✓	✓	✓	✓	✓
$P=3$	✓	✓	✓	✓	✓	✓	✓	✓	✓
$P=4$	$\sim$	✓	✓	✓	✓	✓	✓	✓	✓

Table 3. Table 3 . Uniqueness for inverting the first two moments in the case of an unknown, in-plane uniform ρ 𝜌 \rho , according to band limits. Generically finitely many solutions for ( A , B ) ​  mod SO(3) 𝐴 𝐵  mod SO(3) (A,B)\textup{ mod SO(3)} is denoted by ✓ , infinitely many solutions for A ​  mod SO(3) 𝐴  mod SO(3) A\textup{ mod SO(3)} and B ​  mod SO(3) 𝐵  mod SO(3) B\textup{ mod SO(3)} is denoted by ✗ , and indecisive numerics is denoted by ? .

	$L=2$	$L=3$	$L=4$	$L=5$	$L=6$	$L=7$	$L=8$	$L=9$	$L=10$
$P=0$	✗	✗	✗	✗	✗	✗	✗	✗	✗
$P=1$	✗	✗	✗	✗	✗	✗	✗	✗	✗
$P=2$	✗	✓	✓	✓	✓	✓	?	?	?
$P=3$	✗	✓	✓	✓	✓	✓	✓	?	?
$P=4$	✗	✓	✓	✓	✓	✓	✓	?	?

Table 4. Table 4 . Vanilla, projected and specialized Jacobian tests: these show that a system of polynomial equations generically has only finitely solutions. Notation: F : ℂ N → ℂ M : 𝐹 superscript ℂ 𝑁 → superscript ℂ 𝑀 F\mathrel{\mathop{\mathchar 58\relax}}\mathbb{C}^{N}\rightarrow\mathbb{C}^{M} is a polynomial map, ℂ N = ℂ N ′ ⊕ ℂ N ′′ superscript ℂ 𝑁 direct-sum superscript ℂ superscript 𝑁 ′ superscript ℂ superscript 𝑁 ′′ \mathbb{C}^{N}=\mathbb{C}^{N^{\prime}}\oplus\mathbb{C}^{N^{\prime\prime}} where π ′ superscript 𝜋 ′ \pi^{\prime} , π ′′ superscript 𝜋 ′′ \pi^{\prime\prime} are orthogonal projections onto the factors, and d 1 , d 2 , d 3 subscript 𝑑 1 subscript 𝑑 2 subscript 𝑑 3 d_{1},d_{2},d_{3} are the dimension bounds in Theorem  6 .

	polynomial map	arbitrary choices	linearization	rank check
vanilla	${ℂ}^N\stackrel{𝐹}{⟶}{ℂ}^M$	$q\in {ℂ}^N$	${ℂ}^N\stackrel{dF\left(q\right)}{⟶}{ℂ}^M$	$\text{rank}\left(dF\left(q\right)\right)=N-d_1$
projected	${ℂ}^{N^{\prime }}\oplus {ℂ}^{N^{\prime \prime }}=$ ${ℂ}^N\stackrel{𝐹}{⟶}{ℂ}^M$	$q\in {ℂ}^N$	${ℂ}^{N^{\prime }}\subset {ℂ}^N$ $\stackrel{dF\left(q\right)}{⟶}{ℂ}^M$	$dim\left({\pi }^{\prime }\left(\text{ker}dF\left(q\right)\right)\right)=d_2$
specialized	${ℂ}^{N^{\prime }}\oplus {ℂ}^{N^{\prime \prime }}=$ ${ℂ}^N\stackrel{𝐹}{⟶}{ℂ}^M$	$q^{\prime }\in {ℂ}^{N^{\prime }}$ $q^{\prime \prime }\in {ℂ}^{N^{\prime \prime }}$	${ℂ}^{N^{\prime }}\oplus q^{\prime \prime }$ $\stackrel{dF\left(q^{\prime }\oplus q^{\prime \prime }\right)}{⟶}{ℂ}^M$	$\text{rank }d\left({F|}_{{ℂ}^{N^{\prime }}\oplus q^{\prime \prime }}\right)\left(q^{\prime }\right)$ $=N^{\prime }-d_3$

Equations

$${\tilde{m}}_1(r,\phi )=\frac{1}{n}\sum _{j=1}^n{\hat{I}}_j(r,\phi ),\text{ and }{\tilde{m}}_2(r,\phi ,r^{\prime },{\phi }^{\prime })=\frac{1}{n}\sum _{j=1}^n{\hat{I}}_j(r,\phi ){\hat{I}}_j(r^{\prime },{\phi }^{\prime }),$$

$$\mathcal{P}\varphi (x_1,x_2)\text{mathchar}58={\int }_{-\infty }^{\infty }\varphi (x_1,x_2,x_3)d{x}_3.$$

$$I_j=h_j\ast \mathcal{P}\left(R_j^T\cdot \varphi \right)+{\epsilon }_j,R_j\in \text{SO}\left(3\right),j=1,\dots ,n,$$

$$R_j^T\cdot \varphi (x_1,x_2,x_3)\text{mathchar}58=\varphi (R_j{\left[\begin{array}{ccc}{x}_1& x_2& x_3\end{array}\right]}^T).$$

$$I_j=\mathcal{P}\left(R_j^T\cdot \varphi \right)+{\epsilon }_j,j=1,\dots ,n.$$

$$\widehat{I_j}(x_1,x_2)=\widehat{\mathcal{P}(R_j^T\cdot \varphi )}(x_1,x_2)+\widehat{{\epsilon }_j}=(R_j^T\cdot \hat{\varphi })(x_1,x_2,x_3){|}_{x_3=0}+\widehat{{\epsilon }_j}.$$

$$m_1={\mathbb{E}}_{R\sim \rho }[\widehat{\mathcal{P}\left(R^T\cdot \varphi \right)}],\text{ and }{m}_2={\mathbb{E}}_{R\sim \rho }[\widehat{\mathcal{P}\left(R^T\cdot \varphi \right)}\otimes \widehat{\mathcal{P}\left(R^T\cdot \varphi \right)}].$$

$$\hat{\varphi }\left(\kappa ,\theta ,\phi \right)=\sum _{\ell =0}^L\sum _{m=-\ell }^{\ell }\sum _{s=1}^{S(\ell )}{A}_{\ell ,m,s}{F}_{\ell ,s}(\kappa )Y_{\ell }^m(\theta ,\phi ).$$

$$Y_{\ell }^m(\theta ,\phi )=\sqrt{\frac{(2\ell +1)}{4\pi }\frac{(\ell -m)!}{(\ell +m)!}}{P}_{\ell }^m(\cos \theta )e^{im\phi }$$

$$P_{\ell }^m(x)=\frac{(-1{)}^m}{2^{\ell }\ell !}(1-x^2{)}^{m/2}\frac{d^{\ell +m}}{d{x}^{\ell +m}}(x^2-1{)}^{\ell }.$$

$$R^T\cdot Y_{\ell }^m(x)=Y_{\ell }^m\left(Rx\right)=\sum _{m^{\prime }=-\ell }^{\ell }{U}_{m,m^{\prime }}^{\ell }(R)Y_{\ell }^{m^{\prime }}(x),x\in S^2.$$

$$\rho (R)=\sum _{p=0}^P\sum _{u,v=-p}^p{B}_{p,u,v}{U}_{u,v}^p(R),R\in \text{SO}\left(3\right).$$

$$dR=\frac{1}{8{\pi }^2}\sin (\beta )d\alpha d\beta d\gamma ,$$

$${\hat{I}}_j(\kappa ,\phi )=\sum _{q=-Q}^Q\sum _{t=1}^{T(q)}{a}_{q,t}^j{f}_{q,t}(\kappa )e^{iq\phi }.$$

$$\widehat{I_R}(\kappa ,\phi )=\sum _{q=-Q}^Q\sum _{t=1}^{T(q)}{a}_{q,t}^R{f}_{q,t}(\kappa )e^{iq\phi }.$$

$$\widehat{I_R}(\kappa ,\phi )$$

$$=\sum _{\ell =0}^L\sum _{s=1}^{S(\ell )}\sum _{m=-\ell }^{\ell }{A}_{\ell ,m,s}{F}_{\ell ,s}(\kappa )R^T\cdot Y_{\ell }^m(\frac{\pi }{2},\phi )$$

$$=\sum _{\ell =0}^L\sum _{s=1}^{S(\ell )}\sum _{m=-\ell }^{\ell }\sum _{m^{\prime }=-\ell }^{\ell }{A}_{\ell ,m,s}{F}_{\ell ,s}(\kappa )U_{m,m^{\prime }}^{\ell }(R)Y_{\ell }^{m^{\prime }}(\frac{\pi }{2},\phi ).$$

$$\frac{1}{2\pi }{\int }_0^{\infty }{\int }_0^{2\pi }{f}_{q_1,t_1}(\kappa )e^{i{q}_1\phi }{f}_{q_2,t_2}(\kappa )e^{-i{q}_2\phi }d\phi \kappa d\kappa =\text{mathbbm}{1}_{q_1=q_2}\text{mathbbm}{1}_{t_1=t_2}$$

$$a_{q,t}^R=\sum _{\ell =|q|}^L\sum _{s=1}^{S(\ell )}\sum _{m=-\ell }^{\ell }{A}_{\ell ,m,s}{U}_{m,q}^{\ell }(R){\gamma }_{\ell ,s}^{q,t},$$

$${\gamma }_{\ell ,s}^{q,t}\text{mathchar}58$$

$$=\frac{1}{2\pi }\sqrt{\frac{(2\ell +1)}{4\pi }\frac{(\ell -q)!}{(\ell +q)!}}{P}_{\ell }^q(0){\int }_0^{\infty }{F}_{\ell ,s}(\kappa )f_{q,t}(\kappa )\kappa d\kappa .$$

$${\mathbb{E}}_R[a_{q,t}^R]$$

$$=\sum _{\ell =|q|}^L\sum _{s=1}^{S(\ell )}\sum _{m=-\ell }^{\ell }\sum _{p=0}^P\sum _{u,v=-p}^p{A}_{\ell ,m,s}{B}_{p,u,v}{\gamma }_{\ell ,s}^{q,t}\int U_{m,q}^{\ell }(R)U_{u,v}^p(R)dR$$

$$=\sum _{\ell =|q|}^{\min (L,P)}\sum _{s=1}^{S(\ell )}\sum _{m=-\ell }^{\ell }{A}_{\ell ,m,s}{B}_{\ell ,-m,-q}{\gamma }_{\ell ,s}^{q,t}\frac{(-1{)}^{m+q}}{2\ell +1}.$$

$${\int }_{\text{SO}(3)}\overline{U_{m,n}^{\ell }(R)}{U}_{u,v}^p(R)dR=\frac{1}{2\ell +1}{\text{mathbbm}1}_{\ell =p}{\text{mathbbm}1}_{u=m}{\text{mathbbm}1}_{v=n},$$

$$\overline{U_{u,v}^p(R)}=(-1{)}^{u+v}{U}_{-u,-v}^p(R).$$

$$m_1(q)\text{mathchar}58={\left(\mathbb{E}[a_{q,t}^R]\right)}_{t=1,\dots ,T(q)}=\sum _{\begin{array}{c}\ell \text{mathchar}58|q|\le \ell \le L\\ \ell \equiv q(\text{mod }2)\end{array}}{\Gamma }_{\ell }^q{\mathcal{A}}_{\ell }{\beta }_{\ell }^q.$$

$${\mathbb{E}}_R\left[a_{q_1,t_1}^R{a}_{q_2,t_2}^R\right]=\sum _{{\ell }_1=|q_1|}^L\sum _{s_1=1}^{S({\ell }_1)}\sum _{m_1=-{\ell }_1}^{{\ell }_1}\sum _{{\ell }_2=|q_2|}^L\sum _{s_2=1}^{S({\ell }_2)}\sum _{m_2=-{\ell }_2}^{{\ell }_2}{A}_{{\ell }_1,m_1,s_1}{\gamma }_{{\ell }_1,s_1}^{q_1,t_1}$$

$$\times A_{{\ell }_2,m_2,s_2}{\gamma }_{{\ell }_2,s_2}^{q_2,t_2}\int U_{m_1,q_1}^{{\ell }_1}(R)U_{m_2,q_2}^{{\ell }_2}(R)\rho (R)dR$$

$$\int U_{m_1,q_1}^{{\ell }_1}(R)U_{m_2,q_2}^{{\ell }_2}(R)\rho (R)dR=\sum _{p=0}^P\sum _{u,v=-p}^p{B}_{p,u,v}\int U_{m_1,q_1}^{{\ell }_1}(R)U_{m_2,q_2}^{{\ell }_2}(R)U_{u,v}^p(R)dR.$$

$$U_{m_1,q_1}^{{\ell }_1}(R)U_{m_2,q_2}^{{\ell }_2}(R)=\sum _{{\ell }_3=|{\ell }_2-{\ell }_1|}^{{\ell }_1+{\ell }_2}{\mathcal{C}}_{{\ell }_3}({\ell }_1,{\ell }_2,m_1,m_2,q_1,q_2)U_{m_1+m_2,n_1+n_2}^{{\ell }_3}(R),$$

$${\mathcal{C}}_{{\ell }_3}({\ell }_1,{\ell }_2,m_1,m_2,q_1,q_2)=C({\ell }_1,m_1;{\ell }_2,m_2|{\ell }_3,m_1+m_2)C({\ell }_1,q_1;{\ell }_2,q_2|{\ell }_3,q_1+q_2),$$

Full text

∗Corresponding author: Nir Sharon.

*†*The first two authors contributed equally.

Nir Sharon*∗,†*

School of Mathematical Sciences, Tel Aviv University

Joe Kileel*†*

Program in Applied and Computational Mathematics, Princeton University

Yuehaw Khoo

Department of Statistics and the College, The University of Chicago

Boris Landa

Applied Mathematics Program, Yale University

Amit Singer

Program in Applied and Computational Mathematics and Department of Mathematics, Princeton University

1. Introduction

Single-particle cryo-electron microscopy (cryo-EM) is an imaging method for determining the high-resolution 3-D structure of biological macromolecules without crystallization [25, 35]. The reconstruction process in cryo-EM determines the 3-D structure of a molecule from its noisy 2-D tomographic projection images. By virtue of the experimental setup, each projection image is taken at an unknown viewing direction and has a very high level of noise, due to the small electron dose one can apply to the specimen before inflicting severe radiation damage, e.g., [12, 24, 41]. The computational pipeline that leads from the raw data, given many large unsegmented micrographs of projections, to the 3-D model consists of the following stages. The first step is particle picking, in which 2-D projection images are selected from micrographs. The selected particle images typically undergo 2-D classification to assess data quality and further improve particle picking. At this point, the 3-D reconstruction process begins, where often it is divided into two substeps of low-resolution modeling and 3-D refinement. In this paper, we focus on the mathematical aspects of the former, namely the modeling part. In particular, we suggest using the method of moments (MoM) for ab initio modeling. We illustrate this workflow with an overview given in Figure 1.

The last step in the reconstruction, also known as the refinement step, aims to improve the resolution as much as possible. This refinement process is typically a variant of the expectation-maximization (EM) algorithm which seeks the maximum likelihood estimator (MLE) via an efficient implementation, e.g., [52]. As such, 3-D refinement requires an initial structure that is close to the correct target structure [28, 51]. Serving this purpose, an ab initio model is the result of a reconstruction process which depends solely on the data at hand with no a priori assumptions about the 3-D structure of the molecule [49]. We remark that the two primary challenges for cryo-EM reconstruction are the high level of noise and the unknown viewing directions. Mathematically, without the presence of noise, the unknown viewing directions could be recovered using common lines [61, 62]. Then, the 3-D structure follows, for example, by tomographic inversion, see, e.g., [2]. Reliable detection of common lines is limited however to high signal-to-noise (SNR) ratio. As a result, the application of common lines based approaches is often limited to 2-D class averages rather than the original raw images [56]. Other alternatives such as frequency marching [7] and optimization using stochastic gradient have been suggested [48]. As optimization processes are designed to minimize highly non-convex cost functions, methods like SGD are not guaranteed to succeed. In addition, as in the case of EM, it is not a priori clear how many images are required.

Approximately forty years ago, Zvi Kam proposed a method for 3-D ab initio reconstruction based on computing the mean and covariance of the 2-D noisy images [33]. In order to uniquely determine the volume, the third moment (triple correlation) is also used besides the mean and covariance. In this approach, known as Kam’s method, the 3-D volume is reconstructed without estimating the viewing directions. In this sense, Kam’s method is strikingly different from common lines based approaches and maximum likelihood and other optimization methods that rely on orientation estimation for each image. Crucially, Kam’s method is effective at arbitrary levels of noise, given sufficiently many picked particles for accurate estimation of the moment statistics. Additionally, Kam’s method does not require any starting model, and it requires only one pass through the data to compute moments (contrary to other approaches needing access to the measurements multiple times). Despite the aforementioned advantages, Kam’s method relies on the restrictive assumption that the viewing directions for the images are distributed uniformly over the sphere. This hypothesis, alongside other technical issues, has so far prevented a direct application of Kam’s method to experimental cryo-EM data, for which viewing angles are typically non-uniform [4, 26, 44, 59]. This situation motivates us to explore generalizations of Kam’s method better suited to cryo-EM data.111We remark that Kam’s method, assuming uniform rotations, is of significant current interest in X-ray free electron laser (XFEL) single molecule imaging, where the assumption of uniformity more closely matches experimental reality [21, 45, 65].

In this paper, we generalize Kam’s theory to the case of non-uniform distribution of viewing directions. We regard Kam’s original approach with uniform distribution of viewing angles as a degenerate instance of MoM. In our formulation, we estimate both the 3-D structure and the unknown distribution of viewing angles jointly from the first two moments of the Fourier transformed images. More precisely, for $n$ images $I_{j},j=1,\ldots,n$, the first and second empirical moments of the Fourier transformed images, given in polar coordinates, $\widehat{I}_{j}(r,\varphi),j=1,\ldots,n$, are

[TABLE]

which upon the above discretization become 2-D and 4-D tensors, respectively. Our basic rationale for trying to obtain the volume from the first two moments is as follows. Supposing the distribution of rotations of the image plane to be uniform, then in the limit $n\rightarrow\infty$ the first moment is radially symmetric, that is, it is only a function of $r$ but is independent of $\varphi$. Therefore, $\widetilde{m}_{1}$ may be regarded as a 1-D vector. Similarly, the second moment is a 3-D tensor (rather than 4-D) since it will only depend on $\varphi$ and $\varphi^{\prime}$ through $\varphi-\varphi^{\prime}$ as $n\rightarrow\infty$. Also $I_{j}(r^{\prime},\varphi^{\prime})$ is linearly related to the molecule’s volume via a tomographic projection. Thus, for images of size $N\times N$ pixels, the first and second moments should give rise to $\mathcal{O}(N^{3})$ polynomial equations for the unknown volume and distribution. Assuming the volume is of size $N\times N\times N$ (and the distribution is of lower dimensionality), then the first and second moments have “just” the right number of equations (in terms of leading order) to determine the unknowns. Unfortunately, when the distribution of viewing directions is uniform, as noted by Kam [33], the information encoded in the second moment is algebraically redundant; essentially it is the autocorrelation function (or equivalently, the power spectrum), and this information is insufficient for determining the structure of the molecule. As we will see, a non-uniform distribution of viewing directions introduces additional terms in both the first and second moments, and extends the number of independent equations beyond the autocorrelation case. In particular, we will show that non-uniformity guarantees uniqueness from the analytical counterparts of $\widetilde{m}_{1}$ and $\widetilde{m}_{2}$ in cases of a known distribution, and it guarantees finitely many solutions in other, more realistic, cases of an unknown distribution.

Our work is inspired by several earlier studies on simplified models in a setting called Multi-Reference Alignment (MRA). In MRA, a given group of transformations acts on a vector space of signals [5]. For example, the group ${\text{SO}\negmedspace\left(2\right)}$ acts on the space of band-limited signals over the unit circle by rotating them counterclockwise (as a 1-D analog of cryo-EM). The task then is to estimate a ground truth signal from multiple noisy samples, corresponding to unknown group elements of a finite cyclic subgroup of ${\text{SO}\negmedspace\left(2\right)}$ acting on the signal. The papers [6, 9] show that for a uniform distribution over the group, the signal can be estimated from the third moment, and the number of samples required scales like the third power of the noise variance. On the other hand, for a non-uniform and also aperiodic distributions over the group, the signal can be estimated from the second moment, and the required number of samples scales quadratically with the noise variance [1].

Despite the success of signal recovery in MRA from the first two moments under the action of the cyclic group, it is not apparent that such a strategy is still applicable in the case of cryo-EM. First, in cryo-EM, each image is obtained from the ground truth volume not just by applying a rotation in ${\text{SO}\negmedspace\left(3\right)}$, but also a tomographic projection. Moreover, the studies mentioned above (of MRA) consider finite abelian groups, whereas, in the case of cryo-EM, the group under consideration is the continuous non-commutative group ${\text{SO}\negmedspace\left(3\right)}$. The goal of this paper is then to investigate whether the first and second moment of the images is also sufficient for solving the inverse problem of structure determination in the cryo-EM setting.

1.1. Our contribution

We formulate the reconstruction problem in cryo-EM as an inverse problem of determining the volume and the distribution of viewing directions from the first two moments of the images. Assuming the volume and distribution are band-limited functions, they are discretized by Prolate Spheroidal Wave Functions (PSWFs) and Wigner matrices, respectively. The moments give rise to a polynomial system in which the unknowns are the coefficients of the volume and the distribution. Using computational algebraic geometry techniques [20, 23, 58], we exhibit a range of band limits for the volume and the distribution such that the polynomial system has only finitely many solutions, pointing to the possibility of exact recovery in these regimes. Additionally, we comment on numerical stability issues, by providing condition number formulas for moment inversion. In the setting where the rotational distribution is known, we prove that the number of solutions is generically 1 and present an efficient algorithm for recovering the volume using ideas from tensor decomposition [31]. For the practical case of an unknown distribution, we rely on methods from non-convex optimization and demonstrate, with synthetic data, successful ab initio model recovery of a molecule from the first two moments.

1.2. Organization

The paper is organized as follows. In Section 2, we present discretizations for the volume and distribution and derive the polynomial system obtained from the first two moments. In Section 3, we demonstrate that there exists a range of band limits where the polynomial system for the unknown molecule and distribution has only finitely many solutions. In Section 4, we discuss some implementation details on how the system is solved and present numerical and visual results. Proofs and background material are provided in appendices. For research reproducibility, MATLAB code is publicly available at GitHub.com.222The full address of the GitHub repository is https://github.com/nirsharon/nonuniformMoM.

2. Method of moments

We begin by introducing the image formation model. Then, convenient basis for discretizing various continuous objects, namely the images and the volume (in the Fourier domain) as well as the distribution for orientations, are introduced. From these, relationships between the moments of the 2-D images and the 3-D molecular volume can be derived, enabling us to fit the molecular structure to the empirical moments of the images.

2.1. Image formation model and the 3-D reconstruction problem

In cryo-EM, data is acquired by projecting particles embedded in ice along the direction of the beaming electrons, resulting in tomographic images of the particles. The particles orient themselves randomly with respect to the projection direction. More formally, let $\phi\colon\mathbb{R}^{3}\to\mathbb{R}$ be the Coulomb potential of the 3-D volume, and the projection operator be denoted by $\mathcal{P}\colon\mathbb{R}^{3}\rightarrow\mathbb{R}^{2}$, where

[TABLE]

Assuming the $j$-th particle comes from the same volume $\phi$ but rotated by $R_{j}\in\textup{SO}(3)$, the image formation model is [10, 25]

[TABLE]

where $\varepsilon_{j}$ is a random field modeling the noise term and $h_{j}$ is a point spread function, whose Fourier transform is known as the contrast transfer function (CTF)  [42, 50, 60]. Each image is assumed to lie within the box $[-1,1]\times[-1,1]$. For size $N\times N$ discretized images, we assume the random field $\varepsilon_{j}\sim\mathcal{N}(0,\sigma^{2}I_{N^{2}}),\ j=1,\ldots,n$. Here $R_{j}$ denotes an element in the group of $3\times 3$ rotations ${\text{SO}\negmedspace\left(3\right)}$, and we define the group action by333Here we prefer to write the action of $R^{T}$ and correspondingly later we use Wigner $U$-matrices, instead of $R$ and Wigner $D$-matrices. While simply notational, these conventions allow us to cite identities from [19] verbatim, which are in terms of Wigner $U$-matrices and not Wigner $D$-matrices.

[TABLE]

The rotations $R_{j}$’s are not known since the molecules can take any orientation with respect to projection direction. For the purpose of simplifying the exposition, we shall henceforth disregard the CTF, by assuming

[TABLE]

The presence of CTF is not expected to have a major impact on our main results, and we will incorporate the CTF in a future work. Typically, it is convenient to consider Fourier transform of the images, since by projection slice theorem, the Fourier transform $\widehat{I_{j}}$ of $I_{j}$ gives a slice of the Fourier coefficients $\widehat{\phi}$ of the volume $\phi$:

[TABLE]

The goal of cryo-EM is to recover $\widehat{\phi}$ from the Fourier coefficients of the projections $\widehat{I_{j}}(x_{1},x_{2})$. While reconstructing $\widehat{\phi}$ given estimated $R_{j}$’s amounts to solving a standard computed tomography problem, we wish to reconstruct $\widehat{\phi}$ directly from the noisy images without estimating the rotations, for reasons detailed above. To this end, we assume the rotations are sampled from a distribution $\rho$ on ${\text{SO}\negmedspace\left(3\right)}$, where $\rho\colon{\text{SO}\negmedspace\left(3\right)}\rightarrow\mathbb{R}$ is a smooth band-limited function. Then from the empirical moments of the images $\{\widehat{I_{j}}\}_{j=1}^{n}$, we jointly estimate the volume $\widehat{\phi}$ and the distribution $\rho$.

2.2. Representation of the volume, the distribution of rotations and the images

As mentioned previously, the proposed method of moments consists of fitting the analytical moments

[TABLE]

to their empirical counterparts $\widetilde{m}_{1}$ and $\widetilde{m}_{2}$ as appears in (1) after debiasing.444By the law of large numbers, $\widetilde{m}_{1}\rightarrow m_{1}$ and $\widetilde{m}_{2}\rightarrow m_{2}+\sigma^{2}I$ almost surely as $n\rightarrow\infty$, so $m_{1}$ is fitted to $\widetilde{m}_{1}$ and $m_{2}$ to $\widetilde{m}_{2}-\sigma^{2}I$. For notational convenience, we drop $\sigma^{2}I$ in what follows, either assuming $\tilde{m}_{2}$ has been appropriately debiased already or $\sigma=0$. Through fitting to the empirical moments, we seek to determine the Fourier volume $\hat{\phi}$ and also the distribution $\rho$. In this section, we present discretizations of $\hat{\phi}$ and $\rho$ by expanding them using convenient bases.

2.2.1. Basis for the Fourier volume $\widehat{\phi}$

Since the image formation model involves rotations of the Fourier volume $\widehat{\phi}$, it is convenient to represent $\widehat{\phi}$ as an element of a function space closed under rotations; in fact, this is the same as representing $\widehat{\phi}$ using spherical harmonics (see the Peter-Weyl theorem [19]):

[TABLE]

Here $Y_{\ell}^{m}$ are the (complex) spherical harmonics:

[TABLE]

with associated Legendre polynomials $P_{\ell}^{m}$ defined by:

[TABLE]

In Cartesian coordinates, spherical harmonics are polynomials of degree $\ell$. Without loss of generality, the radial frequency functions $F_{\ell,s}$ should form an orthonormal family (for each fixed $\ell$) with respect to $\kappa^{2}d\kappa$, where $s=1,\ldots,S(\ell)$ is referred to as the radial index. Choices of radial functions suitable for molecular densities include spherical Bessel functions [3], which are eigenfunctions of the Laplacian on a closed ball with Dirichlet boundary condition, as well as the radial components of 3-D prolate spheroidal wave functions  [57].

We assume the volume is band-limited with Fourier coefficients supported within a radius of size $\pi N/2$, i.e., the Nyquist cutoff frequency for the images $I_{j}$’s discretized on a grid of size $N\times N$ (over the square $[-1,1]\times[-1,1])$. Under this assumption, the maximum degree and radial indices $L$ and $S(\ell)$ in (8) are essentially finite. Further details on the particular basis functions $F_{\ell,s}$ and cutoffs $L$ and $S(\ell)$ that we choose to use are deferred to Section A in the appendix. Note that in practice, as we target low-resolution modeling, one can choose to decrease either the cutoff or the grid size to obtain more compact settings. The coefficients $A_{\ell,m,s}\in\mathbb{C}$ furnish our representation of $\widehat{\phi}$ using spherical harmonics. Note that since $\phi$ is real valued, its Fourier transform is conjugate-symmetric, which imposes restrictions on the coefficients $A_{\ell,m,s}$. The specific constraints are presented in Section 4.1.

The advantage of expanding $\widehat{\phi}$ in terms of spherical harmonics is that the space of degree $\ell$ spherical harmonics is closed under rotation; in group-theoretic language, this space forms a linear representation of $\textup{SO}(3)$.555In fact, this is an irreducible representation of $\textup{SO}(3)$ and varying $\ell$ these give all irreps. Thus the action of a rotation on $\widehat{\phi}$ amounts to a linear transformation on the expansion coefficients $A_{\ell,m,s}$ (with a block structure according to $\ell$ and $s$). More precisely, fixing the vector space spanned by $\{Y^{m}_{\ell}(\theta,\varphi)\}_{m=-\ell}^{\ell}$ for a specific $\ell$, the action of a rotation $R$ on this vector space is represented by the Wigner matrix $U^{\ell}(R)\in\mathbb{C}^{(2\ell+1)\times(2\ell+1)}$ (see [19, p. 343]) so that:

[TABLE]

In particular, the matrix $U^{\ell}(R)$ is unitary, with entries degree $\ell$ polynomials in the entries of $R$ [19]. For all $R_{1},R_{2}\in\textup{SO}(3)$ and $\ell$, the group homomorphism property reads $U^{\ell}(R_{1}R_{2})=U^{\ell}(R_{1})U^{\ell}(R_{2})$. In light of (11), 3-D bases of the form $\{F_{\ell,s}(\kappa)Y_{\ell}^{m}(\theta,\phi)\}_{\ell,m,s}$ have been called steerable bases.

2.2.2. Basis for the probability distribution of rotations $\rho$

As we shall see, when expanding the volume in terms of spherical harmonics, the analytical moments (7) involve integrating different monomials of $\{U^{l}(R)\}_{\ell=0}^{L}$ with respect to the measure $\rho(R)dR$. To this end, we assume the probability density $\rho$ over ${\text{SO}\negmedspace\left(3\right)}$ is a smooth band-limited function (and in a function space closed under rotation) by expanding

[TABLE]

By Peter-Weyl, these form an orthonormal basis for $L^{2}(\textup{SO}(3))$, and for higher $p$ they are increasingly oscillatory functions on $\textup{SO}(3)$. Thus, expansion (12) is analogous to using spherical harmonics to expand a smooth function on the sphere, or using Fourier modes for a function on the circle. The cutoff $P\in\mathbb{N}$ is the band limit of the distribution $\rho$; we shall see in the next section that since we use only first and second moments it makes sense to assume $P\leq 2L$. Note that in the special case of a uniform distribution, the only nonzero coefficient is $B_{0,0,0}=1$. Also, $dR$ denotes the Haar measure, which is the unique volume form on the group of total mass one that is invariant under left action. Using the Euler angles parameterization of ${\text{SO}\negmedspace\left(3\right)}$, the Haar measure is of the form

[TABLE]

where the normalizing constant ensures $\int_{\textup{SO}(3)}dR=\int_{\alpha=0}^{2\pi}\int_{\beta=0}^{\pi}\int_{\gamma=0}^{2\pi}dR=1$.

2.2.3. Basis for the 2-D images

At this point, we discuss convenient representations for the images after Fourier transform, $\widehat{I}_{j}$. Similarly to volumes, it is desirable to represent images using a function space closed under in-plane rotations, i.e., $\textup{SO}(2)$. By the Peter-Weyl theorem, this is the same as expanding using Fourier modes, in a 2-D steerable basis:

[TABLE]

Here the radial frequency functions $f_{q,t}$ (for fixed $q$) are taken to be an orthonormal basis with respect to $\kappa d\kappa$, with $\kappa$ referred to as the radial frequency. Comparing to expansion (8) (see Section 2.2), it makes most sense to set $Q=L$. Again, owing to the Nyquist frequency for the discretized images $I_{j}$, we may bound the cutoffs $T(q)$. Typical choices for $f_{q,t}$ for representing tomographic images include Fourier-Bessel functions [66] and the radial components of 2-D prolate spheroidal wave functions  [57]. Details on our specific choices are given in Section A.2 in the appendix.

2.2.4. Choice of radial functions

For the finite expansions in (8) and (14) to accurately represent the Fourier transforms of the electric potential and its slices, one should carefully choose the radial functions $F_{\ell,s}$ and $f_{q,t}$, together with the truncation-related quantities $L$, $S(\ell)$, $Q$, and $T(q)$. In this work, we consider $F_{\ell,s}$ and $f_{q,t}$ to be the radial parts of the three-dimensional and two-dimensional PSWFs [57], respectively. In Appendix A, we describe some of the key properties of the PSWFs, and propose upper bounds for setting $L$, $S(\ell)$, $Q$, and $T(q)$. In practice, band limits would be selected by balancing these expressivity considerations together with the well-posedness and conditioning considerations of Section 3.

2.3. Low-order moments

In this section, we derive the analytical relationship between the first two moments for the observed images $\{a^{j}_{q,t}\}_{j,q,t}$, and the coefficients $\{A_{\ell,m,s}\}_{\ell,m,s}$ and $\{B_{p,u,v}\}_{p,u,v}$ of the volume and distribution of rotations. These relationships will be used to determine $\{A_{\ell,m,s}\}_{\ell,m,s}$ and $\{B_{p,u,v}\}_{p,u,v}$ via solving a nonlinear least-squares problem.

To this end, we first register a crucial relationship between the coefficients of the 2-D images and the 3-D volume. By indexing the images in terms of $R\in{\text{SO}\negmedspace\left(3\right)}$ (instead of $j$ in (14)), we have:

[TABLE]

On the other hand, using the Fourier slice theorem and (11):

[TABLE]

Multiplying (15) and (16) by $f_{q,t}(\kappa)e^{-iq\varphi}$ and integrating against $\frac{1}{2\pi}\kappa d\kappa d\varphi$, then combining the orthogonality relation

[TABLE]

with $Y^{m^{\prime}}_{\ell}(\frac{\pi}{2},\varphi)\propto e^{im^{\prime}\varphi}$, tells us

[TABLE]

where $\gamma_{\ell,s}^{q,t}$ are constants depending on the radial functions:

[TABLE]

From the term $P^{q}_{\ell}(0)$, we see $\gamma_{\ell,s}^{q,t}=0$ if $q\not\equiv\ell$ (mod $2$) (and if $|q|>\ell$ then $\gamma_{\ell,s}^{q,t}\mathrel{\mathop{\mathchar 58\relax}}=0$). Also one may check $\gamma_{\ell,s}^{-q,t}=(-1)^{q}\gamma_{\ell,s}^{q,t}$. Equation (19) connects 2-D image coefficients with 3-D volume coefficients. We note we may as well choose $Q=L$ in (15), since if $|q|>L$ then $a^{R}_{q,t}=0$. In practice, the coefficients $\gamma_{\ell,s}^{q,t}$ are calculated via numerical integration over a closed segment, according to the localization property of the PSWFs, see Appendix A and [39].

2.3.1. The first moment

In this section, from (19) the relationship between the first moment of the images and the volume is derived. Taking the expectation over $R$, and using the distribution expansion (12), we get

[TABLE]

The last equation follows from the orthogonality of the Wigner matrix entries [11, p. 68]

[TABLE]

and

[TABLE]

The first moment gives a set of bilinear forms in the unknowns $\{A_{\ell,m,s}\}_{\ell,m,s}$ and $\{B_{p,u,v}\}_{p,u,v}$, as seen in (24) for each $(q,t)$ with $|q|\leq\min(L,P)$ and $1\leq t\leq T(q)$.

It is convenient to provide compact notation for the first moment formula. To this end, we introduce:

1. (1)

$\mathcal{A}_{\ell}$, a matrix of size $S(\ell)\times(2\ell+1)$ given by $(\mathcal{A}_{\ell})_{s,m}=A_{\ell,m,s}$ 2. (2)

$\beta^{q}_{\ell}$, a vector of size $2\ell+1$ given by $(\beta^{q}_{\ell})_{m}=\tfrac{(-1)^{m}}{2\ell+1}B_{\ell,-m,-q}$ 3. (3)

$\Gamma^{q}_{\ell}$, a matrix of size $T(q)\times S(\ell)$ given by $(\Gamma^{q}_{\ell})_{t,s}=(-1)^{q}\gamma_{\ell,s}^{q,t}$.

Item 2 is zero if $\ell<\lvert q\rvert$ and item 3 is zero if either $\ell<\lvert q\rvert$ or $\ell\not\equiv q\,(\textup{mod }2)$. In this notation, the first moment formula (24) (with fixed $q$ and varying $t$) reads:

[TABLE]

Here $m_{1}(q)\in\mathbb{C}^{T(q)}$ is nonzero only if $\lvert q\rvert\leq\min(L,P)$.

2.3.2. The second moment

Higher moments require higher powers of the image coefficients, and so in the case of the second moment and for $|q_{1}|,|q_{2}|\leq L$, we have

[TABLE]

where

[TABLE]

The product of two Wigner matrix entries is expressed as a linear combination of Wigner matrix entries [19, p. 351],

[TABLE]

where

[TABLE]

is the product of two Clebsch-Gordan coefficients. This product is nonzero only if $(\ell_{1},\ell_{2},\ell_{3})$ satisfy the triangle inequalities. Substituting (31) into (30), and invoking (25) and (26), we obtain:

[TABLE]

where the sum is over $p$ satisfying $\max(|\ell_{1}-\ell_{2}|,|m_{1}+m_{2}|,|q_{1}+q_{2}|)\leq p\leq\min(\ell_{1}+\ell_{2},P)$. Now substituting into (28) gives:

[TABLE]

where the first sum has the range of (28) and the second sum has range of (33). The second moment thus gives a set of polynomials in unknowns $\{A_{\ell,m,s}\}_{\ell,m,s}$ and $\{B_{p,u,v}\}_{p,u,v}$, quadratic in the volume coefficients and linear in the distribution coefficients, namely, the expression in (34) for each $(q_{1},t_{1},q_{2},t_{2})$ with $|q_{1}|\leq L$, $|q_{2}|\leq L$, $\lvert q_{1}+q_{2}\rvert\leq P$, $1\leq t_{1}\leq T(q_{1})$ and $1\leq t_{2}\leq T(q_{2})$. Also, it may be assumed that $P\leq 2L$, since $B_{p,u,v}$ with $p>2L$ does not contribute in either (34) or (24).

As for the first moment, it will be convenient to rewrite the second moment in compact notation. Let us further introduce:

1. (4)

$\mathcal{B}^{q_{1},q_{2}}_{\ell_{1},\ell_{2}}$, a matrix of size $(2\ell_{1}+1)\times(2\ell_{2}+1)$ given by

[TABLE]

where the sum is over $\max(\lvert\ell_{1}-\ell_{2}\rvert,\lvert m_{1}+m_{2}\rvert,\lvert q_{1}+q_{2}\rvert)\leq p\leq\min(\ell_{1}+\ell_{2},P)$ and $C_{p}$ denotes the product Clebsch-Gordan coefficients in (32).

Item 4 is zero if either $\ell_{1}<\lvert q_{1}\rvert$ or $\ell_{2}<\lvert q_{2}\rvert$ or $\max(\lvert\ell_{1}-\ell_{2}\rvert,\lvert q_{1}+q_{2}\rvert)>P$. Now for fixed $q_{1},q_{2}$ and varying $t_{1},t_{2}$, the second moment (34) neatly reads:

[TABLE]

Here $m_{2}(q_{1},q_{2})\in\mathbb{C}^{T(q_{1})\times T(q_{2})}$ is nonzero only if $\lvert q_{1}\rvert,\lvert q_{2}\rvert\leq L$ and $\lvert q_{1}+q_{2}\rvert\leq P$.

3. Uniqueness Guarantees and Conditioning

Here, we derive uniqueness guarantees and comment on intrinsic conditioning for the polynomial system defined by the first and second moments, (27) and (35).

Analysis comes in four cases, according to assumptions on the distribution $\rho$: whether $\rho$ is known or unknown; and if $\rho$ is invariant to in-plane rotations, i.e., $\rho$ depends only on the viewing directions up to rotations that retain the $z$-axis. This invariance restricts $\rho$ to be a non-uniform distribution function over $S^{2}$, see subsection 4.2. If $\rho$ is not invariant to in-plane rotations, we say $\rho$ is totally non-uniform as a distribution on the entire ${\text{SO}\negmedspace\left(3\right)}$. Throughout, our general finding is well-posedness, i.e., the molecule is uniquely determined by first and second moments up to finitely many solutions, under genericity assumptions, for a range of band limits $L$ and $P$. In the case of a known totally non-uniform distribution, we prove the number of solutions is $1$, and give an efficient, explicit algorithm to solve for $\{A_{\ell,m,s}\}$. For all cases, sensitivity of the solution to errors in the moments is quantified by condition number formulas.

3.1. Known, totally non-uniform $\rho$

For this case, we have a provable algorithm that recovers $\{A_{\ell,m,s}\}$ from (27) and (35) (up to the satisfaction of technical genericity and band limit conditions). Remarkably, while the polynomial system is nonlinear (consisting of both quadratic and linear equations), our method is based only on linear algebra. The main technical idea is simultaneous diagonalization borrowed from Jennrich’s well-known algorithm for third-order tensor decomposition [31], that was also used recently for signal recovery in MRA [46].

Theorem 1**.**

The molecule $\{A_{\ell,m,s}\}$ is uniquely determined by the analytical first and second moments, (27) and (35), in the case the distribution $\{B_{p,u,v}\}$ is totally non-uniform, known and $P\geq 2L$, provided it also holds:

1. (i)

The matrices $B_{1}\mathrel{\mathop{\mathchar 58\relax}}=\mathcal{B}^{L,L}_{L,L}$ and $B_{2}=\mathcal{B}^{L,-L}_{L,L}$ of size $(2L+1)\times(2L+1)$ both have full rank, and $B_{1}B_{2}^{-1}$ has distinct eigenvalues. Likewise $B_{3}\mathrel{\mathop{\mathchar 58\relax}}=\mathcal{B}^{L-1,L-1}_{L-1,L-1}$ and $B_{4}=\mathcal{B}^{L-1,1-L}_{L-1,L-1}$ of size $(2L-1)\times(2L-1)$ both have full rank, and $B_{3}B_{4}^{-1}$ has distinct eigenvalues. 2. (ii)

Writing $B_{1}B_{2}^{-1}=\mathrel{\mathop{\mathchar 58\relax}}Q_{12}D_{12}Q_{12}^{-1}$ and $B_{3}B_{4}^{-1}=\mathrel{\mathop{\mathchar 58\relax}}Q_{34}D_{34}Q_{34}^{-1}$ for eigendecompositions, the vectors $b_{12}\mathrel{\mathop{\mathchar 58\relax}}=Q_{12}^{-1}\beta^{L}_{L}$ of size $2L+1$ and $b_{34}\mathrel{\mathop{\mathchar 58\relax}}=Q_{34}^{-1}\beta^{L-1}_{L-1}$ of size $2L-1$ both have no zero entries. 3. (iii)

For $\ell\leq L-2$, the matrix $\mathcal{B}_{\ell,L}^{\ell,L}$ of size $(2\ell+1)\times(2L+1)$ has full row rank. 4. (iv)

For all $\ell$, the matrix $\mathcal{A}_{\ell}$ of size $S(\ell)\times(2\ell+1)$ has full column rank. 5. (v)

For $\ell\geq\lvert q\rvert$ with $\ell\equiv q\,\,(\textup{mod }2)$, the matrix $\Gamma^{q}_{\ell}$ of size $T(q)\times S(\ell)$ has full column rank.

Moreover in this case, there is a provable algorithm inverting (27) and (35) to get $\{A_{\ell,m,s}\}$ in time $\mathcal{O}\left(L^{2}\cdot T^{3}\right)$, where $T\mathrel{\mathop{\mathchar 58\relax}}=\max_{q}T(q)$.

Proof.

For this proof, we need some general properties of the Moore-Penrose pseudo-inverse, denoted by $\dagger$, as in [8]. In particular, if $Y\in\mathbb{C}^{n_{1}\times n_{2}}$ has full column rank and $Z\in\mathbb{C}^{n_{2}\times n_{3}}$ has full row rank, then $Y^{\dagger}Y=I_{n_{2}}$, $ZZ^{\dagger}=I_{n_{2}}$, $(YZ)^{\dagger}=Z^{\dagger}Y^{\dagger}$, and also, pseudo-inversion and transposition commute.

Proceeding, the second moment with $q_{1}=L,q_{2}=L$ tells us:

[TABLE]

and with $q_{1}=L,q_{2}=-L$:

[TABLE]

where $\Gamma_{L}^{L}=(-1)^{L}\Gamma_{L}^{-L}$. We compute $(-1)^{L}$ times the Moore-Penrose psuedo-inverse of (37) and then multiply this on the right of (36). Because $\Gamma_{L}^{L}$ and $\mathcal{A}_{L}$ are each tall with full column rank by assumptions ($v$) and ($iv$), respectively, and $B_{2}$ is invertible by ($i$), properties of the pseudo-inverse imply:

[TABLE]

where we have substituted in an eigendecomposition $B_{1}B_{2}^{-1}=Q_{12}D_{12}Q_{12}^{-1}$. As $B_{1}B_{2}^{-1}$ has distinct eigenvalues by condition ($i$), we see that the eigenvectors of $(-1)^{L}m_{2}(L,L)m_{2}(L,-L)^{\dagger}$ are unique up to scale and given as the columns of $\Gamma^{L}_{L}\mathcal{A}_{L}Q_{12}$. Thus, $\Gamma^{L}_{L}\mathcal{A}_{L}Q_{12}=X\Lambda$, where $X$ consists of eigenvectors of (3.1) and $\Lambda$ is an unknown (as yet) diagonal matrix.

To disambiguate the scales $\Lambda$, we compare with the first moment for $q=L$:

[TABLE]

Multiplying on the left by $X^{\dagger}$ gives $X^{\dagger}m_{1}(L)=\Lambda b_{12}$, an equality of matrix-vector products in which the only unknown is the diagonal matrix $\Lambda$. By the full support of $b_{12}$ (assumption ($ii$)), this determines $\Lambda$. Substituting into $X\Lambda$, we now know $\Gamma^{L}_{L}\mathcal{A}_{L}Q_{12}$. Multiplying on the left by $\Gamma^{L\dagger}_{L}$ and on the right by $Q_{12}^{-1}$ tells us $\mathcal{A}_{L}$.

Backward marching, the second moment with $q_{1}=L-2$ and $q_{2}=L$ reads:

[TABLE]

At this point, we know the first term, and thus the second term gives us $\mathcal{A}_{L-2}$ by appropriately multiplying by pseudo-inverses ($B^{L-2,L}_{L,L}$ is right-invertible by ($iii$)).

Then, we may look at the second moments with $q_{1}=L-4$ and $q_{2}=L$ to similarly determine $\mathcal{A}_{L-4}$, and so on, to $\mathcal{A}_{0}$ or $\mathcal{A}_{1}$ (depending on the parity of $L$). Analogous reasoning and usage of the assumptions gives $\mathcal{A}_{L-1},\mathcal{A}_{L-3},\ldots$

We have provided an algorithm to solve for each $\mathcal{A}_{\ell}$, which proves uniqueness of $\mathcal{A}_{\ell}$ as a byproduct. The time complexity of the algorithm is $\mathcal{O}(L^{2}T^{3})$ since it involves $\mathcal{O}(L^{2})$ matrix operations –matrix multiplications, pseudo-inversions or eigendecompositions – of matrices whose dimensions are all bounded by $T$. (Note that back-substituting to solve for $\mathcal{A}_{\ell}$ involves $\mathcal{O}(L-\ell)$ such matrix operations.) ∎

We remark that condition (iv), which just involves the choice of radial bases, appears to always hold for PSWFs using the cutoffs proposed in Appendix A. Conditions (i), ($ii$) and (iii) just involve the distribution, and are full-rank, spectral and non-vanishing hypotheses. Condition (iv) just involves the molecule and in particular requires $S(L)\geq 2L+1$, which limits $L$ to be less than the Nyquist frequency where $S(L_{\textup{Nyquist}})=1$.

Our algorithm goes by reverse666Reverse frequency marching is natural given the sparsity structure of (35): only $\mathcal{A}_{\ell_{1}}$ and $\mathcal{A}_{\ell_{2}}$ with $\ell_{1}\geq|q_{1}|$, $\ell_{1}\equiv q_{1}\,\,(\textup{mod }2)$ and $\ell_{2}\geq|q_{2}|$, $\ell_{2}\equiv q_{2}\,\,(\textup{mod }2)$ appear in the moments $m_{2}(q_{1},q_{2})$. frequency marching, as we solve for top-frequency coefficients from the second moment (35) where $q_{1},q_{2}=\pm L,\pm(L-1)$ via eigenvectors (similar to simultaneous diagonalization in Jennrich’s algorithm), and then solving for lower-frequency coefficients via linear systems. While our conditions in Theorem 1 are certainly not necessary, fortunately for generic777This means generic with respect to the Zariski topology [30]. Equivalently, there is a non-zero polynomial $p$ in $A,B$ such that $p(A,B)\neq 0$ implies the conditions in Theorem 1 are met. $(A,B)$, those conditions are satisfied, so that the method applies:

Lemma 2**.**

Condition (ii) in Theorem 1 holds for Zariski-generic $\{B_{p,u,v}\}$. If $S(L)\geq 2L+1$, then condition (iii) holds for Zariski-generic $\{A_{\ell,m,s}\}$. At least for $L\leq 100$, conditions (i) and ($iii$) hold for Zariski-generic $\{B_{p,u,v}\}$.

Proof of Lemma 2.

Conditions ($i$)-($iv$) are all Zariski-open, i.e., their failure implies $\{A_{\ell,m,s}\}$ or $\{B_{p,u,v}\}$ obey polynomial equations. As such, to conclude genericity, it suffices to exhibit a single point $\{A_{\ell,m,s}\}$ or $\{B_{p,u,v}\}$, where the conditions are met. For conditions ($i$), ($iii$), we verified the conditions hold at randomly selected points on computer up to $L\leq 100$. Conditions ($ii$) and ($iv$) are obviously generic. ∎

By uniqueness, $A$ is a well-defined function of the first and second moments $m_{1}$ and $m_{2}$ almost everywhere. It is useful to quantify the “sensitivity” of $A$ to errors in $m_{1},m_{2}$, as, e.g., in practice one can access only empirical estimates $\widetilde{m}_{1}$ and $\widetilde{m}_{2}$. An a posteriori (absolute) condition number for $A$ is given by the reciprocal of the least singular value of the Jacobian matrix of the algebraic map:

[TABLE]

Throughout this section, all condition formulas are in the sense of [16, Section 14.3], for which the domain and image of our moment maps are viewed as Riemannian manifolds. To this end, when $\rho$ is unknown, dense open subsets of the orbit spaces $\{(A,B)\,\,\textup{mod SO}(3)\}$, $\{A\,\,\textup{mod SO}(3)\}$, $\{B\,\,\textup{mod SO}(3)\}$ naturally identify as Riemmannian manifolds (for the construction, see [15]).

3.2. Known, in-plane uniform $\rho$

For this case, given a particular image size (and other image parameters), together with band limits $L$ and $P$, we have code888Available in GitHub: https://github.com/nirsharon/nonuniformMoM/JacobianTest. which decides if, for generic $A$ and $B$, the molecule $A$ is determined by (27) and (35), up to finitely many solutions. The basis for this code is the so-called Jacobian test for algebraic maps, see Appendix B. Below is an illustrative computation.

Computational Result 3**.**

Consider $43\times 43$ pixel images, and the following parameters for prolates (representative values): a bandlimit $c$ (see Appendix A) chosen as the Nyquist frequency, 2-D prescribed accuracy (94) set to $\epsilon=10^{-3}$ and 3-D truncation parameter (74) to be $\delta=0.99$ 999The value of $\delta$ means we allow only $1\%$ of the energy to be outside the ball, and is chosen to best model a molecule structure which is assumed to be mostly supported inside a ball.. We varied band limits $L$ in (8) and $P$ in (12), and randomly fixed (12) to give a known in-plane uniform distribution. For each $(L,P)$, we computed the numerical rank of the Jacobian matrix of the polynomial map $m_{B}$ of (41) at a randomly chosen $A$, with random $B$. The Jacobian was convincingly of full numerical rank for a variety of band limits, as seen in Table 1. Cases where the gap between the two least singular values of the Jacobian matrix exceeds a threshold of $10^{6}$ are set as indecisive numerics, and appears in the table as ?. Note that if the rank was calculated in exact arithmetic, this gives a proof that for generic $(A,B)$ generic fibers of the map $m_{B}$ consist of finitely many $A$; i.e., first and second moments (with known in-plane uniform distribution) determine the molecule up to finitely many solutions. For fibers and related definitions, see Appendix B.

Again, the sensitivity of $A$ as a locally defined function of (27) and (35) is quantified by the reciprocal of the least singular value of the Jacobian matrix of $m_{B}$.

3.3. Unknown, totally non-uniform $\rho$

In this case, it is important to note that solutions come in symmetry classes. If $(A,B)$ have specified moments, then so too for $(R\cdot A,R\cdot B)$ for all $R\in\textup{SO}(3)$, that is, we may jointly rotate the molecule and probability distribution and the moments are left invariant. So, solutions come in $3$-dimensional equivalence classes, and we are interested in solutions modulo $\textup{SO}(3)$.

That said, we have code which accepts a particular image size (and other image parameters), together with band limits $L$ and $P$. The code then numerically decides which of the following situations occur: i) for generic $(A,B)$, both $A$ and $B$ are determined by (27) and (35) up to finitely many solutions modulo $\textup{SO}(3)$; ii) for generic $(A,B)$, the molecule $A$ is determined by (27) and (35) up to finitely many solutions modulo $\textup{SO}(3)$, whereas the distribution $B$ has infinitely many solutions; iii) for generic $(A,B)$, both $A$ and $B$ have infinitely many solutions modulo $\textup{SO}(3)$. Note these cases are (essentially) exhaustive, since if $B$ is determined so is $A$ in the regime of Theorem 1. Moreover, we noticed the case ii) really does arise, e.g., this seems to happen when $P=2L$.

Computational Result 4**.**

We keep the running example of $43\times 43$ pixel images, and the prolates parameters of a bandlimit $c$ chosen as the Nyquist frequency, 2-D prescribed accuracy (94) set to $\epsilon=10^{-3}$ and 3-D truncation parameter (74) of $\delta=0.99$. We varied band limits $L$ in (8) and $P$ in (12). For each $(L,P)$, we computed the numerical rank of the Jacobian matrix of the polynomial map

[TABLE]

at a randomly chosen point in the domain. The numerical rank of the Jacobian convincingly equaled three less (that is $d_{1}=3$, see Appendix B) than full column rank for a variety of band limits, see Table 2. Cases where the gap between the third and fourth least singular values of the Jacobian matrix exceeds a threshold of $10^{6}$ are set as indecisive numerics, and appears in the table as ?. If the rank were calculated in exact arithmetic, this furnishes a proof that generic fibers of the map $m$ consist of finitely many $\textup{SO}(3)$-orbits; that is, first and second moments determine both the molecule and the totally non-uniform distribution up to finitely many solutions (modulo global rotation).

For band limits $L$ and $P$ such that generically there are only finitely many solutions for $(A,B)$ mod $\textup{SO}(3)$, the sensitivity of $(A,B)\textup{ mod SO}(3)$ as a (locally defined) function of (27) and (35) is quantified by the reciprocal of the fourth least singular of $m$. For band limits such that generically there are only finitely many solutions for $A$ mod $\textup{SO}(3)$, the sensitivity of $A\textup{ mod SO}(3)$ as a locally defined of (27) and (35) is quantified by the reciprocal of the fourth least singular value of

[TABLE]

where $\dagger$ denotes pseudo-inverse and $\mathcal{P}_{A}$ is the differential of $(A,B)\mapsto A\textup{ mod SO}(3)$. We compute (43) by analytically differentiating (27) and (35), evaluating at $(A,B)$ and place as diagonal blocks of a matrix, and finally applying pseudo-inverse which is SVD-based.

3.4. Unknown, in-plane uniform $\rho$

Again in this case, solutions come in $3$-symmetry classes, orbits under the action of global rotation, so we are interested in solutions modulo $\textup{SO}(3)$. We have code which accepts a particular image size (and other image parameters), together with band limits $L$ and $P$, and numerically decides if for generic $(A,B)$, both $A$ and $B$ are determined by (27) and (35) up to finitely many solutions modulo $\textup{SO}(3)$, or if there are infinitely many solutions. We did not find parameters giving a “mixed” result as in case ii) above.

Computational Result 5**.**

For $43\times 43$ pixel images, and the parameters for prolates (representative values): a bandlimit $c$ chosen as the Nyquist frequency, 2-D prescribed accuracy (94) set to $\epsilon=10^{-3}$ and 3-D truncation parameter (74) of $\delta=0.99$. We varied band limits $L$ in (8) and $P$ in (12), restricting (12) to an in-plane uniform distribution. For each $(L,P)$, we computed the numerical rank of the Jacobian matrix of the polynomial map:

[TABLE]

at a randomly chosen point in the domain. The numerical rank of the Jacobian convincingly equaled three less than full column rank for a variety of band limits, see Table 3. Cases where the gap between the third and fourth least singular values of the Jacobian matrix exceeds a threshold of $10^{6}$ are set as indecisive numerics, and appears in the table as ?. If the rank was calculated in exact arithmetic, this furnishes a proof that generic fibers of the map $m$ consist of finitely many $\textup{SO}(3)$-orbits; that is, first and second moments determine both the molecule and the in-plane uniform distribution up to finitely many solutions (modulo global rotation).

For band limits $L$ and $P$ such that generically there are only finitely many solutions for $(A,B)$ mod $\textup{SO}(3)$, the sensitivity of $(A,B)\textup{ mod SO}(3)$ as a function of moments is quantified by the reciprocal of the fourth least singular of $m$. For example, in the $P=2$ row of Table 3, when evaluating at random $(A,B)$, this worked out to:

[TABLE]

Further, in the $L=4$ column of Table 1, evaluating at random $(A,B)$ gave:

[TABLE]

In practice, we run this refined Jacobian test (takes $<1$ minute on a standard laptop) to identify well-conditioned band limits $L$ and $P$ before we attempt non-convex optimization.

4. Numerical Optimization and First Visual Examples

After studying the theoretical properties of the polynomial system which is defined by the first two moments, we discuss in this section aspects of numerically inverting the polynomial map via optimization.

4.1. Incorporating natural constraints in optimization

When determining the coefficients $A=\{A_{\ell,m,s}\}_{\ell,m,s}$ and $B=\{B_{p,u,v}\}_{p,u,v}$, the search space has to be restricted in order to ensure the coefficients stem from some physical volume and density.

4.1.1. Constraints on the volume

To ensure the volume $\phi\colon\mathbb{R}^{3}\rightarrow\mathbb{R}$ is a real-valued function, one has to ensure its Fourier transformation $\widehat{\phi}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{3}\rightarrow\mathbb{C}$ satisfies conjugate symmetry $\hat{\phi}(\kappa,\theta,\varphi)=\bar{\hat{\phi}}(\kappa,\pi-\theta,\pi+\varphi)$. That is, in spherical coordinates,

[TABLE]

Assuming the basis $\{F_{\ell,s}\}$ is a set of real-valued functions, along with the facts that $\overline{Y^{m}_{\ell}(\theta,\varphi)}=(-1)^{m}Y^{-m}_{\ell}(\theta,\varphi)$ and $Y^{m}_{\ell}(\pi-\theta,\pi+\phi)=(-1)^{\ell}Y^{m}_{\ell}(\theta,\phi)$, we get

[TABLE]

This further implies

[TABLE]

Having such relationships, $\{A_{\ell,m,s}\}_{\ell,m,s}$ can thus be written in terms of some real coefficients $\{\alpha_{\ell,m,s}\}_{\ell,m,s}$ as:

[TABLE]

The latter means that instead of solving a complex optimization problem in terms of the coefficients $A_{\ell,m,s}$, one can work with the real coefficients $\alpha_{\ell,m,s}$ of  (46). Otherwise, the equality constraints (45) are required.

4.1.2. Constraints on the density

Similarly, to ensure the density $\rho$ being a real-valued function, we need to ensure

[TABLE]

The fact that $\overline{U^{p}_{u,v}(R)}=(-1)^{v-u}U^{p}_{-u,-v}(R)$ leads to

[TABLE]

Again, from such relationships, it can be shown that an alternative to (48) can be written in terms of real coefficients $\beta_{p,u,v}$:

[TABLE]

Here, $\prec_{\textup{lex}}$ is the lexicographical order, that is $(u_{1},v_{1})\prec_{\textup{lex}}(u_{2},v_{2})$ iff $u_{1}<u_{2}$ or both $u_{1}=u_{2}$ and $v_{1}<v_{2}$.

Two additional constraints are required. First, the integral of any density function is one. To ensure such a correct normalization, we simply let

[TABLE]

which means it is no longer considered as unknown. Finally, the nonnegativity of the density is ensured via a collocation method, that is requiring

[TABLE]

for $R_{i}$’s on a near uniform, refined grid on ${\text{SO}\negmedspace\left(3\right)}$. While (51) does not prevent the density from becoming negative off the ${\text{SO}\negmedspace\left(3\right)}$ grid, requiring the density to be non-negative entirely on ${\text{SO}\negmedspace\left(3\right)}$ leads to an optimization problem that is much more costly to solve in practice. Note that we do not enforce positivity of $\rho$ by requiring it to be a sum-of-squares, as, e.g., already in the case of an in-plane uniform distribution on the sphere $S^{2}\subset\mathbb{R}^{3}$, not all nonnegative polynomials may be written as a sum-of-squares, see Motzkin’s example when $P=6$ [43].

4.2. Accommodating invariance to in-plane rotations

While molecules typically exhibit preferred orientations, there is no physical reason why molecules should have preferred in-plane orientations. In this section, we focus on the case of non-uniform rotational distributions invariant to in-plane rotations since these distributions better model real cryo-EM data sets.

For simplicity, we fix the image plane as perpendicular to the $z$-axis. We add the prior that the density for drawing $R$ equals the density for drawing $Rz(\alpha)$, for all $R\in{\text{SO}\negmedspace\left(3\right)}$ and all rotations $z(\alpha)$ of $\alpha\in\mathbb{R}$ radians about the $z$-axis. This assumption reads

[TABLE]

Therefore,

[TABLE]

Here we used the group representation property of $U^{p}$. Checking explicitly the action of $z(\alpha)$ on degree $p$ spherical harmonics,

[TABLE]

So continuing the above,

[TABLE]

This is equivalent to $B_{p,u,v}=0$ for $v\neq 0$ where $v$ ranges over $-p,-p+1,\ldots,p$. To sum, we have found that in-plane invariance is captured by:

[TABLE]

For a sanity check, a distribution with in-plane invariance should sample a rotation with density only depending on which point maps to the north pole. Namely, $\rho(R)$ should only depend on the last column of $R$, that is, $R(\mathrel{\mathop{\mathchar 58\relax}},3)=R_{\bullet 3}$. Indeed, this holds as $U^{p}_{u0}(R)=(-1)^{u}\sqrt{\frac{4\pi}{2l+1}}\overline{Y^{u}_{p}}(R_{\bullet 3})$ [19, Eqn. 9.44, Pg. 342].

Restricting the expansion of $\rho$ as above, we easily see the first moment is independent of $\varphi$. It is now merely a linear combination of basis functions $F_{\ell,s}(\kappa)$. Likewise, for the second moment, angular dependency is only on the difference $\varphi_{1}-\varphi_{2}$, meaning it is a linear combination of basis functions $e^{im(\varphi_{1}-\varphi_{2})}F_{\ell_{1},s_{1}}(\kappa_{1})F_{\ell_{2},s_{2}}(\kappa_{2})$. Thus, in subsection (3.4), we have the following polynomial map, now with fewer $B$ variables and fewer invariants than in subsection (3.3)

[TABLE]

4.3. Direct method – known totally non-uniform distribution

For the “easy” case of a known, totally non-uniform distribution, we have implemented the provable algorithm in Theorem 1. The method’s performance is illustrated by way of an example. As the ground truth volume, we use EMD-0409, that is, the catalytic subunit of protein kinase A bound to ATP and IP20 [32], as presented at the online cryo-EM data-bank [38]. The volumetric array’s original dimension is $128$ voxels in each direction, which we downsampled by a factor of three to $43$. The volume was expanded using PSWFs with a band limit $c$ chosen to be the Nyquist frequency and 3-D truncation parameter (74) of $\delta=0.99$. Before downsampling, the full expansion consists of degree $L=40$; with downsampling and proper truncation, we aim to recover the terms up to degree $L=7$. For the known totally non-uniform distribution, we took $P=14$ (per Theorem 1), and then formed a particular distribution using a sums-of-squares. Precisely, we formed a random linear combination of Wigner entries up to degree $7$, multiplied this by its complex conjugate, invoked (26) and (32) to rewrite the result as a linear combination of Wigner entries up to degree $14$, repeated for a second square, added, and finally normalized to satisfy (50). Then, with the distribution known as such, the volume contributes $1080$ unknowns (without discounting for (45)). Providing the algorithm with $m_{1}$ and $m_{2}$, our method took $0.24$ seconds on a standard laptop, and recovered the unknowns $A$ up to a relative error in $\mathcal{L}^{2}$ norm of $5.4\times 10^{-11}$. Visual results are in Figure 2.

4.4. Setting up a least-squares formulation

For the cases where we lack a direct method, we formulate the problem in terms of minimizing a least-squares cost function. First, we define the unknowns of our optimization process to be the coefficients of the volume $A=\{A_{l,m,s}\}$ and distribution $B=\{B_{p,u,v}\}$. The explicit formulas (27) and (35) provide means to write the low-order moments (7) as functions of our unknown coefficients, that is $m_{1}=m_{1}(A,B)$ and $m_{2}=m_{2}(A,B)$.

In practice, given data images, one estimates the low-order statistics using the empirical moments $\widetilde{m}_{1}$ and $\widetilde{m}_{2}$ of (1), but now given in PSWFs coordinates

[TABLE]

The connection between the empirical moments and their analytical formulas as functions of our unknowns gives rise to a nonlinear least-squares

[TABLE]

where $\lambda$ is a parameter chosen to balance the errors from both terms. In particular, two main considerations determine the value of $\lambda$. First is the number of elements in each summand. Namely, the second moment includes many more entries than the first moment. Therefore, without the effect of noise, $\lambda$ is set to be the ratio between the number of entries in first moment and the second moment. The second factor to balance is the different convergence rates of the empirical moments, see also [1]. The nonlinear least-squares (60) may be adjusted to incorporate the constraints on $\{A_{l,m,s}\}$ and $\{B_{p,u,v}\}$ that ensure $\phi$ is a real-valued volume and $\rho$ a probability density.

We remark that it is interesting to consider pre-conditioners, or more intricate weighings, in the formation of the nonlinear least-squares cost (60). Such might alleviate high condition numbers observed in Section 3, and potentially accelerate optimization algorithms. While we have not tested a pre-conditioner in optimization experiments yet, one possibility would be to consider the following normalized cost:

[TABLE]

Effectively, (61) scales each polynomial in $(A,B)$ given by $m_{1}$ and $m_{2}$ to take value $1$.

4.5. Complexity analysis of inverting the moments via gradient-based optimization

Before moving forward to further numerical examples, we state the computational load of minimizing the least-squares cost function (60). It is worth noting that in many modern ab initio algorithms, like SGD [48] and EM [52], the runtime of each iteration is measured with respect to the size of the set of data images, which can be huge. In our approach, we only carry out one pass over the data to collect the low-order statistics. In here, we assume the empirical moments are already given, and so the complexity of each iteration is merely a function of the size of the moments or equivalently depends on the size and resolution of the data images, as reflected by their PSWF representations.

Many possible algorithms exist to minimize the least squares problem (60), for example direct gradient descent methods, such as trust-region [47], or alternating approaches, including alternating stochastic gradient descent. Here, we present the complexity of evaluating the cost function and its gradient, regardless of the specific algorithm or implementation one wishes to exploit.

For simplicity, denote by $S$ and $T$ two bounds for the radial indices $S(\ell)$ and $T(q)$ of the 3-D and 2-D PSWF expansions, respectively. Typically, it is sufficient to take $S=S(0)$ and $T=T(0)$, as radial degree decreases as overall degree ($\ell$) increases.

Starting from the first moment (27): with a fixed $\ell$ we have to apply two matrix-vector products in a row which requires an order of $\mathcal{O}\left(S\ell+TS\right)$ arithmetic operations. The variable $\ell$ increases up to $L$, which sums up to a total of $L\cdot\mathcal{O}\left(S\ell+TS\right)=\mathcal{O}\left(LS(L+T)\right)$. The gradient uses the precomputed remainder $m_{1}(A,B)-\left(\widetilde{m}_{1}\right)$ and is calculated by two terms with similar complexity as the above. Namely, the cost of both evaluation and gradient calculations is again $\mathcal{O}\left(LS(L+T)\right)$.

For the second moment, we follow (35): establishing $\Gamma^{q}_{\ell}\,\mathcal{A}_{\ell}$ is done in $\mathcal{O}\left(TSL\right)$ and applying the product in $\mathcal{O}\left(TL^{2}\right)$. Overall, the evaluation is bounded by

[TABLE]

The gradient is a bit more complicated, in short, there are two terms for the volume derivatives and one term for the distribution part, with the precomputed remainder $m_{2}(A,B)-\left(\widetilde{m}_{2}\right)$ we get an overall complexity of $\mathcal{O}\left(L^{2}S(L^{2}+T^{2}+TL)\right)$. In summary, the first moment requires third-order complexity with respect to the different parameters where the second moment requires a total power of five.

Finally, the parameters $T$, $S$, and $L$ can be described by the PSWF representation: the length $L$ of the 3-D PSWF expansion and the bound on the radial indices $S$ are related to the parameter $c$ of sampling rate, and are bounded according to (A.1.1). Additional bound, now on the radial 2-D expansion $T$, uses the accuracy parameter $\epsilon$ of the 2-D images and the above $L$ as given in (94). For more details on those parameters, see Appendix A.

4.6. Remark on using semidefinite programming (SDP) relaxation

Solving the nonlinear least-squares problem in Eq. (60) could suffer from slow convergence because the cost function is a polynomial of degree 6. We remark that in principle, it is possible to apply a semidefinite programming relaxation to facilitate the optimization. For convenience, let the second moments $m_{2}(A,B)_{q_{1},t_{1},q_{2},t_{2}}$ be summarized as

[TABLE]

where $G_{q_{1},t_{1},q_{2},t_{2}}(\cdot)$ is a linear operator that captures the RHS of Eq. (34). If we define

[TABLE]

the optimization problem can be written as

[TABLE]

To deal with the non-convex constraint $\bar{A}=AA^{T}$, we propose the following relaxed constraint

[TABLE]

which gives the following non-linear least squares problem

[TABLE]

Comparing with (60), although (65) is still a non-convex problem, the degree of the polynomial in the cost function of (65) is 4 (instead of 6). Furthermore, one can solve (65) efficiently by minimizing $(A,\bar{A})$ and $B$ in an alternating fashion. Therefore if at the optimum $\bar{A}\approx AA^{T}$ in spite of the relaxation (64), solving (65) can be advantageous.

We remark on the special case when the density coefficient $B$ is given. In this situation, one can consider an SDP relaxation

[TABLE]

The nuclear norm minimization strategy as in matrix completion [17] is used to promote $\bar{A}$ to be of rank-1. We test the SDP in (4.6) when given a fixed $B_{0}$. We generate $B_{0}$ for a non-uniform distribution from a 6-th degree nonnegative polynomial over the rotation group, i.e. letting $P=6$. We generate a volume with random coefficients $A_{0}$ with $L=3$. Noise is added to the moments in the following manner:

[TABLE]

Where

[TABLE]

and

[TABLE]

In this case, we set in (4.6),

[TABLE]

The stability results in recovering $A_{0}$ are shown in Figure 3. We ran five simulations for every $\epsilon$ and average the relative error

[TABLE]

Results show an exact recovery in the noiseless case and slowly increasing in relative error as $\epsilon$ grows.

4.7. Volume from moments – non-uniform vs. uniform

As a first numerical example, we present a recovery comparison between the cases of uniform and non-uniform distributions of rotations. In this example, we use as a ground truth a low degree approximation of a mixture of six Gaussians, given in a non-symmetric conformation. The approximation, which we ultimately use as our reference, is attained by discretizing the initial volume to $23\times 23\times 23$ and truncating the PSWFs expansion to $L=4$. This expansion consists of $118$ coefficients in total. The other PSWFs parameters that we use are a band limit $c$ that corresponds to the Nyquist frequency and 3-D truncation parameter (74) of $\delta=0.99$. The original volume and its approximation appear in Figure 4.

We divide the example into two scenarios of different distributions, uniform and non-uniform. In each case, we start from the analytic moments (7), calculated with respect to 2-D prescribed accuracy (94) of $\epsilon=10^{-3}$, and obtain an estimation based on minimizing the least squares cost function (60). The optimization is carried with a gradient-based method, specifically we use an implementation of the trust-region algorithm, see e.g., [47]. In the first case, we use as the distribution of rotations a quadratic expansion $P=2$ which is in-plane uniform. Based on the in-plane invariance, we present this distribution as a function on the sphere in Figure 5. For the second case, we use a uniform distribution of rotations.

In both cases, we let the optimization reach numerical convergence, where the progress in minimization is minor. In this example, it is usually at about $100-150$ iterations. In the case of non-uniform distribution, we observe that choosing a random initial guess can have an effect on the speed of convergence but has almost no influence on the resulted volume. In other words, we gain numerical evidence for uniqueness. The estimated volume, in this case, is depicted on the left side of Figure 6.

On the other hand, in the case of a uniform distribution, while convergence was typically quicker than in the non-uniform case, the results vary between different initial guesses, indicating the richness of the space of possible solutions. One such solution appears on the right side of Figure 6. This behavior of the optimization solver agrees with our previous knowledge on the ill-posedness of Kam’s method and also with the Jacobian test which shows degree deficiency of the polynomial system defined by the first and second moment under the uniform distribution.

4.8. Comparing volumes using FSC

A commonly used cryo-EM resolution measure is the Fourier shell correlation (FSC) [29]. The FSC measures cross-correlation coefficient between two 3-D volumes over each corresponding shell. That is, given two volumes $\phi_{1}$ and $\phi_{2}$, the FSC in a shell $\kappa$ is calculated using all voxels $\bm{\kappa}$ on this $\kappa$-th shell:

[TABLE]

Customary, the resolution is determined by a cutoff value. The threshold question is discussed in [64], where in our case since we wish to compare a reconstructed volume against its ground truth, we use the $0.5$ threshold. Since we focus on ab initio modeling, we aim to estimate a low-resolution version of the molecule from the first two moments. Thus, we expect the cutoff to reach a value which ensures a good starting point for a refinement procedure.

4.9. Visual example and the effect of non-uniformity

We next introduce an example for the most realistic scenario of an unknown, in-plane uniform distribution, by inverting the moment map of a real-world structure through minimization of a least-squares cost function (60). In this example, we once again illustrate the feasibility of numerically approaching the solution, without any prior assumption on the volume.

The example is constructed as follows. As the ground truth volume, we once again use EMD-0409, the catalytic subunit of protein kinase A bound to ATP and IP20 [32], as presented at the online cryo-EM data-bank [38]. The map original dimension is $128\times 128\times 128$ voxels. Since we aim to recover a low-resolution model, we reduce complexity and downsample it by a factor of three to $43$. We firstly expand this volume using PSWFs with a band limit $c$ chosen as the Nyquist frequency and 3-D truncation parameter (74) of $\delta=0.99$. The full expansion consists of degree $L=40$, and after truncation to maximize conditioning, as done in Section 3.3, we aim to recover the low degree counterpart up to degree $L=6$. The moments were calculated with respect to 2-D prescribed accuracy (94) of $\epsilon=10^{-3}$ and in the absence of noise. The volume contributes $657$ unknowns to be optimized.

As the ground truth distribution, we choose three different functions: uniform, highly non-uniform and a non-uniform case in-between. The two non-uniform cases are cubic spherical harmonics expansions ($P=3$) and satisfy in-plane invariance and so we present them in Figure 7 as functions on the sphere, together with a histogram to compare and illustrate their “non-uniformness”. The non-uniform distributions add extra $15$ unknowns which means that, in total, we optimize $672$ unknowns in the cases of non-uniform distribution and only $657$ unknowns in the case of uniform distribution.

In the optimization process, we use the limit of the empirical moments (59) ($n\to\infty$) as our input moments. As before, we use a trust-region algorithm, see e.g., [47], which is a gradient-based method. To fix the initialization between the different cases, we start the search with the zero volume. In cases of non-uniform distribution, we provide a random non-uniform distribution to start with. Our method is implemented in MATLAB R2017b, and we calculated the example on a laptop with a 2.9 GHz Intel Core i5 processor and 16 GB 2133 MHz memory.

The result we present next is obtained after $60$ iterations of trust-region, each iteration usually uses up to $30$ inner iterations to estimate the most accurate step size. The runtime of this example is about $55$ minutes for each model, where at this point, our naive implementation does not support any parallelization which potentially can lead to a significant improvement in the total runtime. For example, the evaluation of the second moment and its associated gradient part are related to the leading complexity term as described in Section 4.5. Their implementation is based upon matrix product as seen in the form (35). This part can remarkably benefit from parallel execution. Note that evaluating the PSWF functions, as well as the product Clebsch-Gordan coefficients (which appears in the moments), are all calculated offline as a preprocessing step.

We present a comparison between the different FSC curves for the three cases. As implied by Figure 8, the resolution increases (lower FSC cut) as the non-uniformity becomes more significant. Specifically, with the uniform distribution we obtain merely $39.1\textup{\AA}$, where for the two other non-uniform cases we get $22.5\textup{\AA}$ and $19.0\textup{\AA}$ as the non-uniformity increases in the examples of Figure 8.

A visual demonstration of the output of the optimization is presented in Figure 9, where we plot side by side the ground truth and three models, from the uniform to the most non-uniform one.

4.10. Recovery from noisy images

We conclude this section with an example of recovering a volume from its noisy projection images. The volume is a mixture of six Gaussians, synthetically designed to have no spatial symmetry. The volume’s size is $15\times 15\times 15$ and its full PSWF expansion is of length $L=13$, with band limit $c$ chosen as the Nyquist frequency and 3-D truncation parameter (74) of $\delta=0.99$. We use an in-plane uniform distribution of rotations, very localized on a $45$ degree cone, represented with an expansion length of $P=3$. The distribution function is shown on Figure 10 and can model a realistic scenario of highly anisotropic viewing directions (see, e.g., [4]). Using the distribution, we generated projection images to obtain $200,000$ observations. These images were then contaminated with noise. The SNR of an image $I_{j}$ with the noise term $\varepsilon_{j}$ is $\operatorname{SNR}_{j}=\mathinner{\!\left\lVert I_{j}-\varepsilon_{j}\right\rVert}^{2}/\mathinner{\!\left\lVert\varepsilon_{j}\right\rVert}^{2}$, using the Frobenius norm. The noise was chosen to achieve an average SNR value of $1/3$. Three examples of clean images and their noisy versions are depicted in Figure 11. As seen in Figure 11, the projections are hardly noticed in the noisy images for the naked eye.

We expand the noisy images using a 2-D PSWF basis, as appears in (14). Then, the coefficients and their double-products are averaged to estimate the first and second moments as in (59). The reconstruction uses the empirical moments to estimate the volume and distribution. For the volume, our gradient-based least-squares algorithm targets its full expansion, which consists of $192$ unknowns. The unknown distribution includes $8$ unknowns spherical harmonics coefficients. We reached the result we present next very quickly, starting from a random initial guess. It took about $15$ iterations of trust-region; each iteration could use up to $30$ inner iterations to estimate the most accurate step size. The runtime of this example is less than $10$ minutes.

A visual demonstration of the estimated volume is provided in Figure 12. We present the estimation, side by side, to the original volume. As seen in the various pictures, the reconstruction, while not perfect, captures most features and the general shape of the structure. This encouraging result indicates that inverting the moments is possible also from noisy moments and that the mapping has some robustness to small perturbations.

5. Discussion and Conclusion

The method of moments offers an attractive approach for modeling volumes in cryo-EM. This statistical method completely bypasses the estimation of viewing directions by treating them directly as nuisance parameters. The assumption of a non-uniform distribution of viewing angles enables in many cases volume estimation using only the first and second moments of the data. This phenomenon opens the door for fast, single-pass reconstruction algorithms, based on inverting the map from the volume and distribution to the low-order statistics of the projection images.

This paper extended Zvi Kam’s original method of moments for cryo-EM to the setting of a non-uniform distribution of viewing directions. We formulated the reconstruction problem using appropriate discretizations for the images, the volume, and the distribution. Then, we derived moment formulas using properties of the spherical harmonic functions and Wigner matrix entries. Computational algebra was employed to analyze the resulting large-scale system of polynomial equations. The analysis shows the seeming complication of an unknown, non-uniform distribution renders 3-D reconstruction easier than in the uniform case, as now only first and second moments are required to determine a low-resolution expansion of the molecule, up to finitely many solutions. Intermediate cases were treated; remarkably, when the distribution is known and totally non-uniform over ${\text{SO}\negmedspace\left(3\right)}$, there is an efficient, provable algorithm to invert the first and second moments non-linearly. Additionally, our work addressed several numerical and computational aspects of the method of moments. An implementation of a trust-region method was presented and used to illustrate the advantages of our approach over Kam’s classical approach by numerical experiments involving synthetic volumes.

We regard our work as a definite, albeit initial step toward developing the method of moments for ab initio modeling from experimental datasets. Firstly, even in the synthetic cases considered here, further work on the optimization side is warranted. Variations on our nonlinear cost function that incorporate a pre-conditioner, e.g., (61), could be considered. Secondly, other techniques for large-scale nonlinear least squares optimization should be tried, such as Levenberg-Marquardt [40] or Variable Projection [18], where in the latter one can exploit the linearity in the moments with respect to the distribution, by eliminating out the distribution. Thirdly, to get our method working on images, further effects, such as the CTF and imperfect centering of picked particles, should be incorporated into the moment formulas. Fourthly, accurate covariance estimation in high dimensions requires eigenvalue shrinkage [22], the theory for which may call for a modification in the non-uniform setting.

To simplify our exposition, we have stuck to the asymmetric and homogeneous cases here, although both of these can be relaxed in the method of moments. Specifically, as already noted in Kam’s original paper [33], point symmetries of molecules are reflected in the vanishing of certain expansion coefficients, see also [63]. Therefore, MoM can be reformulated using fewer coefficients for symmetric molecules. This fact, alongside with further improvement of the representation of the distribution, may pave the way for recovery under practical cases of very restricted viewing angles, as reported in literature [4, 26, 44, 59]. At the same time, heterogeneity, at least if it is finite and discrete, can be expressed using a mixture of volumes and a corresponding mixture of moments, see [14, 5]. In future work, computational algebra should be applied to these cases to check whether the first and second moments remain sufficient for unique recovery.

To conclude, we raise one further possibility, in some sense at odds with the message of this paper. In the non-uniform case, we have determined that the first and second moments are sufficient information-theoretically for volume recovery. Nonetheless, the resulting optimization landscape is potentially challenging, due to non-convexity or ill-conditioning. Thus, despite the increased statistical cost of estimating the third moment, it seems worthwhile to ask what can be gained computationally by reprising the third moment in MoM (or at least, using a carefully chosen slice of the third moment). Specifically, we would like to answer this question: can the third moment facilitate more efficient modeling at higher resolution?

Acknowledgments

The authors thank Nicolas Boumal, Peter Bürgisser, Eitan Levin, Dilano Saldin and Yoel Shkolnisky for stimulating conversations, and the anonymous referees for their valuable comments.

This research was supported in parts by Award Number R01GM090200 from the NIGMS, FA9550-17-1-0291 from AFOSR, Simons Foundation Math$+$X Investigator Award, the Simons Collaboration on Algorithms and Geometry, the Moore Foundation Data-Driven Discovery Investigator Award, and NSF BIGDATA Award IIS-1837992. BL’s research was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement 723991 - CRYOMATH).

Appendix A Prolate Spheroidal Wave Functions

Here we describe key properties of the PSWFs, and propose a method for setting the expansion parameters $L$, $S(\ell)$, $Q$, and $T(q)$. We begin with the three-dimensional PSWFs, where we describe important properties established in the literature [57, 34, 53], and outline our choice for setting $L$ and $S(\ell)$, accordingly. Then, we proceed with a short analogous description for the two-dimensional PSWFs (summarizing results of [57]), and derive a method for choosing $Q$ and $T(q)$ by directly exploiting the fact that the images to be expanded are tomographic projections of a bandlimited and localized volume function (employing our previous representation for the volume function).

A.1. Volume function representation with three-dimensional PSWFs

Let $\Phi\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{3}\rightarrow\mathbb{R}$ be a square integrable (volume) function on $\mathbb{R}^{3}$, representing the true underlying electric potential of the molecule, and denote by $\hat{\Phi}$ its three-dimensional Fourier transform. It is common practice to assume that $\Phi(x)$ is bandlimited (i.e., $\hat{\Phi}$ is restricted to a ball) while being localized in space. Functions satisfying this property are naturally represented by three-dimensional PSWFs, as detailed next.

We say that the function $\Phi(x)$ as $c$-bandlimited if $\hat{\Phi}(\omega)$ vanishes outside a ball of radius $c$. That is, $\Phi$ is $c$-bandlimited if

[TABLE]

where $\mathbf{B}$ is the unit ball. Among all $c$-bandlimited functions, the three-dimensional PSWFs on $\mathbf{B}$ [57] are the most energy concentrated in $\mathbf{B}$, while constituting an orthonormal system over ${\mathcal{L}^{2}(\mathbf{B})}$. Namely, they satisfy

[TABLE]

for $i=1,2,\ldots$, i.e., $\Psi_{1}$ is the most energy concentrated $c$-bandlimited function, $\Psi_{2}$ is the most energy concentrated $c$-bandlimited function orthogonal to $\Psi_{1}$, and so on. Three-dimensional PSWFs can be obtained as the solutions to the integral equation

[TABLE]

where we denote the solutions (the PSWFs with bandlimit $c$) as $\Psi_{\ell,m,s}^{c}(x)$ and their corresponding eigenvalues as $\alpha_{\ell,m,s}^{c}$, where the enumeration over $i$ in (68) is replaced with an enumeration over the triplet $\ell,m,s$ described below, and the eigenvalues appear in non-increasing ordering with respect to the original enumerate $i$. $\Psi_{\ell,m,s}^{c}(x)$ and $\alpha_{\ell,m,s}^{c}$ together form the eigenfunctions and eigenvalues of (69), with $m\in\mathbb{Z}$, $\ell\in\mathbb{N}\cup\left\{0\right\}$, and $s\in\mathbb{N}$. Furthermore, the functions $\Psi_{\ell,m,s}^{c}(x)$ are orthogonal on both $\mathbf{B}$ and $\mathbb{R}^{3}$ using the standard $\mathcal{L}^{2}$ inner products on $\mathbf{B}$ and $\mathbb{R}^{3}$, respectively, and are dense in both the class of $\mathcal{L}^{2}(\mathbf{B)}$ functions and in the class of $c$-bandlimited functions on $\mathbb{R}^{3}$. In spherical coordinates, the functions $\Psi_{\ell,m,s}^{c}(x)$ agree with the form in the right-hand side of (8), and can be expressed as

[TABLE]

where $Y^{m}_{\ell}(\theta,\varphi)$ are the spherical harmonics (see (9)). Numerical evaluation of the three-dimensional PSWFs (in particular of the radial part $F_{\ell,s}^{c}$) was considered in [39].

From the properties of the three-dimensional PSWFs mentioned above, any volume function $\Phi(x)\in\mathcal{L}^{2}(\mathbb{R}^{3})$ can expanded in $\mathbf{B}$ as

[TABLE]

where $\overline{(\cdot)}$ denotes complex conjugation. Next, we consider the truncation of the expansion in (71), where it is convenient to bound the resulting truncation error in terms of the assumed spatial localization of $\Phi(x)$. Towards this end, we say that the function $\Phi(x)$ is $\varepsilon$-concentrated if

[TABLE]

Additionally, we define the normalized eigenvalues

[TABLE]

where we mention that $0\leq\lambda^{c}_{\ell,m,s}\leq 1$, $\lambda^{c}_{\ell,m,s}=\lambda^{c}_{\ell,0,s}$ for all triplets $(\ell,m,s)$, and $\lambda^{c}_{\ell,m,s}\underset{s\rightarrow\infty}{\longrightarrow}0$ for every $\ell$. Now, we propose to set $S(\ell)$ according to

[TABLE]

where $\delta\in(0,1)$ is some constant, and set $L$ to be the largest $\ell$ for which $S(\ell)$ is defined (i.e., such that the set $\left\{s\mathrel{\mathop{\mathchar 58\relax}}\;\lambda_{\ell,0,s}^{c}\geq\delta\right\}$ is non-empty). Correspondingly, the volume function resulting from the truncating the expansion in (71), according to the chosen $S(\ell)$ and $L$, is

[TABLE]

The following proposition bounds the error of approximating $\Phi(x)$ by $\phi(x)$.

Proposition 1**.**

Let $\Phi(x)$ be $c$-bandlimited with a unit $\mathcal{L}^{2}\left(\mathbf{B}\right)$ norm and assume it is $\varepsilon$-concentrated. Then,

[TABLE]

The proof follows immediately from Theorem 5 in [34] and from our choices of $S(\ell)$ and $L$. It is evident that the approximation error in the right-hand side of (76) can be made arbitrarily small by taking $\delta$ sufficiently small. Furthermore, in the case that $\Phi(x)$ is localized in space, i.e., $\varepsilon\ll 1$, we can take $\delta$ to be large, possibly even close to $1$, and still get approximation errors sufficiently small for our purposes.

A.1.1. Length of the expansion

Clearly, the number of basis functions taking part in the expansion (75), which is given explicitly by $\sum_{\ell=0}^{L}\sum_{m=-\ell}^{\ell}S(\ell)$, depends on the number of normalized eigenvalues $\lambda^{c}_{\ell,m,s}$ exceeding $\delta$. In this respect, the normalized eigenvalues $\lambda^{c}_{\ell,m,s}$ are known to admit the following three distinct regions of behavior (when sorted in descending order). The first is called the “flat region”, where $\lambda^{c}_{\ell,m,s}$ take values very close to $1$, the second is called the “transitional region”, where $\lambda^{c}_{\ell,m,s}$ shift rapidly from values close to $1$ to values close to [math], and the third is called the “decay region”, where $\lambda^{c}_{\ell,m,s}$ are very close to [math] and exhibit a super-exponential decay rate. As for the number of basis functions chosen according to (74), the following holds [53]:

[TABLE]

where the first, second, and third terms on the right-hand side of (A.1.1) correspond to the number of normalized eigenvalues $\lambda^{c}_{\ell,m,s}$ exceeding $\delta$ from the flat region, the transitional region, and the decay region of the eigenvalues, respectively. Clearly, the asymptotically dominant term is $\mathcal{O}(c^{3})$, which corresponds to the number of terms in the expansion chosen from the flat region. Additionally, we need an extra $\mathcal{O}(c^{2}\log{(c)})$ terms if we take $\delta$ to be small (note that the second term in the right hand-side of (A.1.1) is negative for $\delta>0.5$, meaning that asymptotically we need less than ${c^{3}}/{4.5\pi}$ terms for values of $\delta$ close to $1$). The remaining $o(c^{2}\log{(c)})$ terms from the decay region are negligible compared to the leading asymptotic terms.

A.1.2. Fourier domain representation

Up to this point, we have shown that three-dimensional PSWFs are naturally adapted for expanding a volume function $\Phi(x)$ which is bandlimited and localized in space, where we provided an appropriate error bound (76). However, note that in (8) we actually expand the Fourier transform of the molecular potential. We now connect our previous expansion of $\Phi(x)$ with the expansion of its Fourier transform, and show that in fact (and uniquely for PSWFs) the two coincide, in the sense that expanding a function in three-dimensional PSWFs is equivalent to expanding its Fourier transform in three-dimensional PSWFs (after an appropriate scaling and dilation). Let $\hat{\Psi}_{\ell,m,s}$ denote the three-dimensional Fourier transform of $\Psi_{\ell,m,s}$, then by (69) it is easy to verify that

[TABLE]

where $\mathbf{1}_{c\mathbf{B}}(\omega)$ is the indicator function on $c\mathbf{B}$. It is evident that the Fourier transform of each three-dimensional PSWF is equal to itself up to a constant factor, a dilation by $c$, and a restriction to a ball of radius $c$. Consequently, by taking the Fourier transform of (75) we have

[TABLE]

where

[TABLE]

We conclude this part as follows. Given a bandlimit $c$ (typically chosen as the Nyquist frequency corresponding to the projection images’ resolution), we take the radial part $F_{\ell,m,s}(k)$ of (8) as $F_{\ell,m,s}^{c}(k/c)/{c^{3/2}}\cdot\mathbf{1}_{c}(k)$, where $\mathbf{1}_{c}(k)$ is the indicator function on $[0,c]$, and $F_{\ell,m,s}^{c}(r)$ is the radial part of the three-dimensional PSWFs on $\mathbf{B}$ (the factor $1/{c^{3/2}}$ ensures that $F_{\ell,m,s}(k)$ are orthonormal over $[0,\infty)$ w.r.t the measure $k^{2}dk$). Then, setting $S(\ell)$ according to (74) for a given parameter $\delta$ allows for the controlled approximation error (76).

A.2. Projection image representation with two-dimensional PSWFs

In the sequel, we are interested in providing a suitable representation for the projection images of the rotated copies of $\phi(x)$. By the Fourier slice theorem, the two-dimensional Fourier transforms of such projections are equal to slices from the three-dimensional Fourier transform of $\phi(x)$ (i.e., of $\hat{\phi}(\omega)$). Therefore, if $\phi(x)$ is $c$-bandlimted, then the projection images are bandlimited to a disk of radius $c$. Additionally, we expect the projection images to be localized in the unit disk if $\phi(x)$ is sufficiently localized in the unit ball. For such projection images, two-dimensional PSWFs are expected to provide a natural representation (see [36]).

We briefly summarize properties of the two-dimensional PSWFs which are used in our context. In essence, the properties of the two-dimensional PSWFs are analogous to those of the three-dimensional PSWFs when replacing the unit ball $\mathbf{B}$ with the unit disk $\mathbf{D}$. Let $P\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{2}\rightarrow\mathbb{R}$ be a square integrable function on $\mathbb{R}^{2}$, representing a tomographic projection of $\phi$. We say that $P(x)$ as $c$-bandlimited if its two-dimensional Fourier transform, denoted by $\hat{P}(\omega)$, vanishes outside a disk of radius $c$. That is, $P$ is $c$-bandlimited if

[TABLE]

Among all $c$-bandlimited functions, the two-dimensional PSWFs on $\mathbf{D}$ are the most energy concentrated in $\mathbf{D}$, that is, they satisfy (68) when replacing $\mathbf{B}$ with $\mathbf{D}$, while constituting an orthonormal system over ${\mathcal{L}^{2}(\mathbf{D})}$. The two-dimensional PSWFs were derived and analyzed in [57], and were shown to be the solutions to the integral equation

[TABLE]

We denote the PSWFs with bandlimit $c$ as $\psi_{q,t}^{c}(x)$, and their corresponding eigenvalues as $\beta_{q,t}^{c}$, which together form the eigenfunctions and eigenvalues of (82), with $q\in\mathbb{Z}$, and $t\in\mathbb{N}$. Furthermore, the functions $\psi_{q,t}^{c}(x)$ are orthogonal on both $\mathbf{D}$ and $\mathbb{R}^{2}$ using the standard $\mathcal{L}^{2}$ inner products on $\mathbf{D}$ and $\mathbb{R}^{2}$, respectively, and are dense in both the class of $\mathcal{L}^{2}(\mathbf{D)}$ functions and in the class of $c$-bandlimited functions on $\mathbb{R}^{2}$. In polar coordinates, the functions $\psi_{q,t}^{c}(x)$ agree with the form in the right-hand side of (14), and can be expressed as

[TABLE]

where the eigenfunctions $\psi_{q,t}^{c}(x)$ are normalized to have an $\mathcal{L}^{2}(\mathbf{D})$ norm of $1$. Numerical evaluation of the two-dimensional PSWFs was considered in [54].

From the properties of the two-dimensional PSWFs mentioned above, any function $P(x)\in\mathcal{L}^{2}(\mathbb{R}^{2})$ can be expanded in $\mathbf{D}$ as

[TABLE]

Now, considering the truncated expansion

[TABLE]

we are interested in controlling the error

[TABLE]

From (82), the Fourier transform of $\psi_{m,k}$ can be expressed as

[TABLE]

where $\mathbf{1}_{c\mathbf{D}}(\omega)$ is the indicator function on $c\mathbf{D}$, which is analogous to the relation between the three-dimensional PSWFs $\Psi_{\ell,m,s}^{c}$ and their Fourier transforms $\hat{\Psi}_{\ell,m,s}^{c}$ in (78). Continuing, taking the Fourier transform of (85) gives

[TABLE]

where

[TABLE]

We will now relate 2D basis representation error to that of the 3D basis functions. Comparing the 2D expansion (88) with the relation between 2-D and 3-D coefficients (19), while employing (89) and (80) we have

[TABLE]

for $|q|\leq L$, where $\widetilde{a}_{q,t}=0$ for $|q|>L$, and

[TABLE]

where $\gamma_{\ell,s}^{q,t}$ is from (20). Using the Cauchy–Schwarz inequality, we can write

[TABLE]

where we also used the fact that $U^{\ell}(R)$ is a unitary matrix. Finally, taking $Q=L$ and assuming w.l.o.g that $\left\|\phi\right\|_{\mathcal{L}^{2}(\mathbf{B})}=1$, we obtain from (86) and (92) that

[TABLE]

Given a prescribed accuracy $\epsilon$, for every $-L\leq q\leq L$ we choose $T(q)$ to be the smallest integer such that

[TABLE]

which results in

[TABLE]

where $\eta_{\ell,s}^{q,t}$ are computed by evaluating $\gamma_{\ell,s}^{q,t}$ of (20) via numerical integration (using Gauss-Legendre quadratures). Note that the right-hand side of (94) is determined by the decay rate of $\eta_{\ell,s}^{q,t}$ in $t$, which is dominated by the decay rate of the the eigenvalues of the two-dimensional PSWFs $\beta_{q,t}$. Those are known to admit a rapid decay in the form of a super-exponential decay rate following a certain transitional region (see [13, 53]). Hence, if $T(q)$ is sufficiently large then (94) can be satisfied for an arbitrarily small $\epsilon$ with a marginal increase in the number of required terms. Last, we mention that when provided with images sampled on a Cartesian grid, the coefficients $\widetilde{a}_{q,t}$ can be approximated accurately from the images by fast algorithms [36, 37].

Appendix B Linearizing polynomial maps with the Jacobian matrix

In this section, we describe the linearization technique from computational algebraic geometry we used to obtain the uniqueness results in Tables 1, 2, 3 from Section 3. The first paper to apply algebraic geometry techniques to cryo-EM was [5]. Nevertheless, similar Jacobian tests have been used in other applications such as for testing rigidity in sensor network localization, see e.g., [27] and testing whether a matrix can be completed into a low-rank matrix [55].

To state the method, we fix $\mathbb{C}^{N}=\mathbb{C}^{N^{\prime}}\oplus\mathbb{C}^{N^{\prime\prime}}$, let $\pi^{\prime}$ and $\pi^{\prime\prime}$ be projection onto the factors, and consider a polynomial map $F=(F_{1},\ldots,F_{M})\mathrel{\mathop{\mathchar 58\relax}}\mathbb{C}^{N}\rightarrow\mathbb{C}^{M}$ (that is, each coordinate function $F_{i}=F_{i}(x_{1},\ldots,x_{M})$ is a polynomial on $\mathbb{C}^{N}$). While $F$ is generally a nonlinear map, its first derivative at $q\in\mathbb{C}^{N}$ is a linear map represented by the Jacobian matrix

[TABLE]

In addition, we define the fiber in $q\in\mathbb{C}^{N}$ by

[TABLE]

and the projected fiber by

[TABLE]

For $q^{\prime}\in\mathbb{C}^{N^{\prime}}$ and $q^{\prime\prime}\in\mathbb{C}^{N^{\prime\prime}}$, define the specialized fiber by

[TABLE]

Because $F$ is described by polynomials, there is a tight relationship between the dimension of fibers of $F$ (as algebraic varieties) and the dimension of the kernels of $dF$ (as linear spaces). This is summarized by the Jacobian tests below. Somewhat remarkably, the linear algebra tests are done at a single point in the domain of $F$, but imply algebraic geometric statements for almost all points in the domain of $F$.

Theorem 6**.**

Suppose it is known that, generically, the fiber, projected fiber and specialized fiber have dimensions $\geq d_{1},d_{2},d_{3}$, respectively (if we have no such knowledge, then $d_{1}=d_{2}=d_{3}=0$). Choose particular points $q_{0}\in\mathbb{C}^{N}$, $q^{\prime}_{0}\in\mathbb{C}^{N^{\prime}}$ and $q^{\prime\prime}_{0}\in\mathbb{C}^{N^{\prime\prime}}$.

• $-$

Vanilla Jacobian test:* if $\textup{rank}\,dF(q_{0})=N-d_{1}$, then generic fibers have dimension exactly $d_{1}$.*

• $-$

Projected Jacobian test:* if $\dim\,\pi^{\prime}(\textup{ker}\,dF(q_{0}))=d_{2}$, then generic projected fibers have dimension exactly $d_{2}$.*

• $-$

Specialized Jacobian test:* if $\textup{rank }\,d\left(F|_{\mathbb{C}^{N^{\prime}}\oplus\,q_{0}^{\prime\prime}}\right)(q_{0}^{\prime})=N^{\prime}-d_{3}$, then generic specialized fibers have dimension exactly $d_{3}$.*

Several technical remarks are in order. Firstly, in Theorem 6, the fiber, projected fiber and specialized fiber are affine algebraic varieties and hence a dimension is defined for each of their irreducible components according to [20]. The meaning of the theorem is that each component has dimension exactly $d_{1},d_{2},d_{3}$, respectively. Crucially, affine algebraic varieties have finitely many components. Thus the theorem implies “finitely many solutions” up to symmetries, if the symmetries give $d_{1},d_{2},d_{3}$-dimensional ambiguities, respectively. Secondly, “generic” in Theorem 6 is with respect to the Zariski topology. Concretely, there exists some polynomial $G$ on $\mathbb{C}^{N}$ such that for all $q\in\mathbb{C}^{N}$ with $G(q)\neq 0$ the implications in the theorem hold. In particular, any property that holds generically holds on a Lebesgue full measure subset of points. Thirdly, the Jacobian ranks in Theorem 6 take on generic values, as each minor of the relevant matrix is a polynomial in $q,q^{\prime},q^{\prime\prime}$.

Theorem 6 states rigorous conclusions if the Jacobian rank tests are passed. On the other hand, if the tests fail for $q_{0},q_{0}^{\prime},q_{0}^{\prime\prime}$, and $q_{0},q_{0}^{\prime},q_{0}^{\prime\prime}$ were drawn randomly from any continuous distribution on $\mathbb{F}^{N}$, then by genericity of the Jacobian ranks, with probability $1$, the generic fibers, projected fibers, or specialized fibers of $F$ have dimension strictly more than $d_{1}$, $d_{2}$, or $d_{3}$.

We applied the specialized test in subsection 3.2 with $d_{1}=0$, the vanilla and projected tests in subsection 3.3 with $d_{1}=d_{2}=3$ and the vanilla test in subsection 3.4 with $d_{1}=3$. The settings of $3$ reflect the fact, in the latter two subsections, that the fibers are $\textup{SO}(3)$-sets and we are interested in solutions modulo global rotation. The bounds may be seen as instances of the orbit-stabilizer theorem, see [5, Proposition 4.11]. When the Jacobian rank tests were passed, this meant that, generically, there are only finitely many solutions up to global ambiguities.

In practice, we ran the Jacobian tests in floating-point arithmetic and used SVD for robust rank estimation. Namely, we looked at multiplicative gaps between consecutive singular values, and regarded any gap exceeding a predefined threshold ($10^{6}$) as evidence that all lower singular values should be regarded as zero. While these computations fall short of a fully rigorous mathematical proof due to the possibility of rounding errors in floating-point arithmetic, it was typically evident which singular values ought to be counted as zero or non-zero.