# A Survey of the Use of Iterative Reconstruction Algorithms in Electron Microscopy

Authors: C. O. S. Sorzano; J. Vargas; J. Otón; J. M. de la Rosa-Trevín; J. L. Vilas; M. Kazemi; R. Melero; L. del Caño; J. Cuenca; P. Conesa; J. Gómez-Blanco; R. Marabini; J. M. Carazo
Journal: BioMed Research International
 (2017)

---

## Abstract

One of the key steps in Electron Microscopy is the tomographic reconstruction of a three-dimensional (3D) map of the specimen being studied from a set of two-dimensional (2D) projections acquired at the microscope. This tomographic reconstruction may be performed with different reconstruction algorithms that can be grouped into several large families: direct Fourier inversion methods, back-projection methods, Radon methods, or iterative algorithms. In this review, we focus on the latter family of algorithms, explaining the mathematical rationale behind the different algorithms in this family as they have been introduced in the field of Electron Microscopy. We cover their use in Single Particle Analysis (SPA) as well as in Electron Tomography (ET).


---

## Body

## 1. Introduction

Electron Microscopy has been established as one of the key players in Structural Biology with the goal of elucidating the three-dimensional structure of macromolecular complexes in order to better understand their function and molecular mechanisms [1–3]. One of the most important steps in the image processing pipeline is the 3D reconstruction of a map compatible with the projections acquired at the microscope [4]. In practice, projections of the macromolecule are contaminated by a huge amount of noise (typical Signal-to-Noise Ratios in the order of 0.01; i.e., there is 100 times more noise power than signal power, Hosseinizadeh et al. [5]), and the 3D reconstruction emerges, in a simplified way, as the 3D “average” of thousands of projections, each one looking at the molecule from a different point of view [4]. In the 3D space, the signal coming from the macromolecules is reinforced by the averaging process. However, random noise tends to be canceled by this averaging. Currently, the 3D reconstruction step is no longer seen as a limiting step (except for its execution time) in Single Particle Analysis due to the large number of particles involved in the reconstruction (between tens and hundreds of thousands), and direct Fourier inversion methods are currently the standard de facto [6–8]. These latter methods are especially well suited to handle a large number of projections thanks to their computational speed and their accuracy when the angular coverage of the set of projections fully fills the 3D Fourier space, which is currently normally the case in SPA (a word of caution should be expressed in those cases in which subsequent rounds of 3D classification significantly reduce the number of images per 3D class). However, in the past, there was an intense work of research in selecting the best reconstruction algorithm between four different families of reconstruction algorithms: (1) direct Fourier inversion [6–8], (2) back-projection algorithms [9–18], (3) iterative algorithms [19–43], and (4) Radon inversion methods [44–46].The situation in Electron Tomography (ET) is, in general, specially involved, due to the combined effect of a smaller number of projection images, the existence of a missing wedge, and the fact that images are very large, posing an important challenge to the computational resources (especially in terms of memory; for a comparison between Single Particle Analysis and Electron Tomography, see Jonic et al. [47]). The research of 3D reconstruction algorithms trying to reduce reconstruction artifacts and the execution time has been rather intense, especially in the last decade [48–61]. Similarly, pure reconstruction algorithms in ET may be combined with a priori information incorporating the fact that the reconstruction is sparse in some domain or exploiting the discrete nature of the object being imaged [62–76].In this review, we focus on the family of iterative reconstruction algorithms, also known as series expansion methods. The classical algorithms (ART, Block ART, and SIRT) have been followed by a number of more modern algorithms like Conjugate Gradient, Subgradient Descent, Projected Subgradient, Superiorization, ADMM, and so forth, allowing all kinds of regularizations with a special emphasis on sparsity promoting regularizers. We start by showing how the 3D reconstruction problem can be posed as a problem of solving a linear system of equations. Then, we put the more classical methods in a common algebraic framework. Finally, we introduce the rationale behind the more modern methods.


## 2. 3D Reconstruction as an Equation System Problem

Images collected by the electron microscope can be understood as the parallel projection of the Coulomb potential of the macromolecule being imaged. The relationship between the 3D model,V(r~):R3×1→R, and a projection image, I(s~):R2×1→R, is given by [77](1)IA~s~=∫-∞∞VA~-1H~Ts~dt,where r~=(x,y,z,1) and s~=(x,y,1) are the homogeneous coordinates of a 3D point and a 2D point, respectively; H~T is(2)H~T=10001000t001,with t being the integration variable, and(3)A~=r11r12r13txr21r22r23tyr31r32r3300001=Rt0T1is the matrix that specifies the point of view from where the projection has been taken (R is a rotation matrix normally specified by 3 Euler angles and t is an in-plane displacement). Every projection has its own matrix A~ reflecting the different projection directions and individual shifts of each image.In practice, this ideal image is never observed but it is corrupted by random noise (whose nature is related to the structure of the ice surrounding the macromolecule, the random arrival of electrons, etc.) and the projected image is not known at any positions~ but at a discrete set of positions (s~i, normally the centers of the pixels of the acquired image). In this way, we may rewrite the above equation as(4)IA~,i=IA~s~i=∫-∞∞VA~-1H~Ts~idt+ni.Note that this image formation model uses additive noise, since the main source of noise is not the low number of electron counts at each pixel but the structure of the amorphous ice embedding the molecule. This noise has a Gaussian distribution as was shown in [78].Let us assume now that we express the volume as a linear combination of functions,b(r~), shifted to a set of known positions r~j (e.g., usually these positions are regularly distributed on a grid [30]):(5)Vr~=∑jxjbr~-r~j.The basis functions may be voxels or any other function with some interesting property from the tomographic point of view (for instance, Kaiser-Bessel modified functions, also known as blobs, [32, 36] are known to reduce the noise in the 3D reconstruction).The goal of the tomographic problem is to determine the basis coefficients,xj, such that the experimental projection is actually the line integral of the volume V. Substituting (5) into (4) and disregarding the noise term (it will be recovered later), we have(6)IA~,i=∫-∞∞VA~-1H~Ts~idt=∫-∞∞∑jxjbA~-1H~Ts~i-r~jdt=∑jxj∫-∞∞bA~-1H~Ts~i-r~jdt=∑jxjaA~,ij,where aA~,ij is the line integral of the basis function located at r~j at the location s~i along the direction determined by the geometrical transformation A~, that is, the contribution of the j-th basis function of the volume to the i-th pixel in the image. We may write the image pixels in a single vector using lexicographical order and we may rewrite the above equation in a matrix form:(7)IA~=AA~x. If we have images of size P×P pixels and assuming that the volume is reconstructed in a cube of size B×B×B, then IA~∈RP2 is a vector in a space of dimension P2, x∈RB3 is a vector in a space of dimension B3, and AA~ is a matrix of size P2×B3. If we look again at the pixel level, we may write(8)IA~,i=aA~,iTx=aA~,i,x,where aA~,iT is the i-th row of the matrix AA~. This is the equation of a hyperplane in the RB3 space, meaning that every single pixel of the projection image provides a hyperplane constraint for locating the volume coefficients x. The true x, in the ideal noiseless case, must be at the intersection of all the hyperplanes defined by all pixels (the noise observed in the measurements simply randomly shifts each hyperplane along the direction perpendicular to the hyperplane normal). In this way, a single projection provides P2 equations, but we have B3 unknowns so that there is an infinite number of structures compatible with a single projection. In Electron Microscopy (EM), we collect thousands of projection images (let us say M) and we may stack all the measurements in a big vector. The corresponding equation system would then be(9)IA~1IA~2⋮IA~M=AA~1AA~2⋮AA~Mxwhich in general can be written as the linear equation system(10)b=Ax, where b∈RMP2 is a vector collecting all the experimental images and A is a matrix of size MP2×B3 collecting all the projection matrices. Note that this matrix is very sparse since, for common basis functions like voxels, blobs, and so forth, just a few coefficients in the volume contribute to each pixel (only those coefficients along the integration line passing through that pixel). Just to give some typical numbers in EM, Let us say we have M=50,000 images of size P×P=256×256. We would have MP2=3.3 billion equations and B3=P3=16.8 million unknowns (if voxels are used as basis function). The equation system is inconsistent because of the measurement noise and some kind of least squares solution must be sought [79]. We will do so in the next section.We have presented the linear equation system only in connection to the data collection geometry. However, the Contrast Transfer Function (CTF) of the microscope (i.e., how the microscope blurs the ideal images) can easily be incorporated. The convolution can be represented by a matrix multiplication using a Toeplitz matrix, such that instead of theAA~ matrices that contain purely geometrical information we may use the matrices CAA~, where C is the CTF Toeplitz matrix (each experimental image may have its own CTF matrix).Interestingly, posing the tomographic problem as a linear equation system can also be done in Fourier space. Thanks to the Central Slice Theorem, the relationship between the 3D Fourier transform of the macromolecule model and the 2D Fourier transform of the projection can be expressed as [77](11)I^A~S~=e-iS~,H0RtV^A~′TH~0TS~,where A~′=R00T1, H~0 is the matrix defined at (2) with t=0, and H0 is the same matrix except the last row (the one adding its homogeneous nature). S~ is the 2D frequency coordinate and I^A~ and V^ are the 2D and 3D Fourier transforms of the projection image and the macromolecular model, respectively. Equation (11) means that, to evaluate the Fourier coefficient at the 2D frequency coordinate S~, we need to evaluate the Fourier transform of the volume at the coordinate A~′TH~0TS~ and then multiply by the corresponding phase term, e-iS~,H0Rt, to account for the shift in the images. In practice, we do not have the Fourier transform of the volume at any possible location. We will have to interpolate its value from the 3D Fourier coefficients in its surroundings. Let us concentrate at the frequency S~j and let us call the corresponding 2D Fourier coefficient I^j. Let us refer to the j-th 3D Fourier coefficient by V^j. Then, we may interpolate the 2D Fourier coefficient as(12)I^i=∑jV^ja^A~,ij,where the terms a^A~,ij comprise the phase shift as well as the interpolation weights. This equation is formally identical to (6), meaning that the 3D reconstruction problem in Fourier space can also be expressed as a linear system of equations, and we have to be careful to perform the complex interpolation. Actually, this is the position taken in Scheres [80] and Chen and Förster [59]. Including the CTF in this framework is even easier than in the real space case, since it suffices to multiply the a^A~,ij terms by the ci CTF coefficients.


## 3. (Weighted) Least Squares Solution

Given the inconsistent equation systemAx=b, we may try to look for the x that minimizes its distance to all the hyperplanes given by the experimentally measured pixel values. This is achieved by minimizing the norm of the residuals vectors Ax-b. In tomographic terms, Ax may be interpreted as the expected experimental projections given a macromolecular model x, while b are the experimentally measured projections. The idea is to find the model x such that its reprojections are as similar as possible to the experimentally acquired images. The residual of the comparison between the experimental images and the reprojections should be just the noise present in the experimental images:(13)x=arg⁡minx⁡Ax-b2, where the norm is calculated using the standard inner product:(14)u2=u,u=uTu. This is a linear least squares problem. Note that the solution is not unique (any vector x0 in the kernel of the matrix A will yield the same error norm, because Ax+x0=Ax+0=Ax). In practice, x0 refers to any structure whose Fourier coefficients are in areas not measured by the experimental images. Developing the norm in the optimization problem, we have(15)x=arg⁡minx⁡Ax-bTAx-b=arg⁡minx⁡b2+xTATAx-2xTATb. If we now differentiate and equate to 0, we have(16)∂Ax-b2∂x=2ATAx-2ATb=0 or equivalently(17)ATAx=ATb.This latter equation is known as the normal equation of the least squares problem. It can be shown that any solution of the normal equations is also a solution of the least squares problem, and, conversely, any solution of the least squares problem must also be a solution of the normal equations [81, Chapter 6]. Except for degenerate cases, ATA is a positive definite matrix, implying that the matrix is invertible, and consequently the normal equations have a unique solution. For complex numbers, the norm is defined as(18)u2=uHu with H being the conjugate transpose, so that the normal equations in Fourier space would be(19)AHAx=AHb.We may now introduce weights in the minimization if there are some measurements that are more reliable than others (in real space, this is more difficult, but in Fourier space, low frequency components are overrepresented with respect to high frequency ones and they are consequently downweighted, Chen and Förster [59]). Given a diagonal, positive definite matrix W of size MP2×MP2, we may define the W inner product as(20)u,vW=uHWv. The normal equations become in this case(21)AHWAx=AHWb.The most straightforward approach to solving this equation system is the use of the Moore-Penrose pseudoinverse:(22)x=AHWA-1AHWb=A†b. The matrix AHWA is of size B3×B3 and its direct inversion is normally out of the reach of any current computer. Numerical algorithms that avoid the need to invert this matrix are normally used and they are explained in the next section.


## 4. Iterative Solutions of a Linear Equation System

Note that all equation systems posed so far (see (6), (12), (17), (19), and (21)) are of the generic form:(23)Ax=b.In this section, A and b are the generic system matrix and data terms, independently of their actual expressions in terms of the microscopic projection images, CTFs, Fourier or real space versions, and so forth, since these details have already been carefully presented in the sections above.A very simple, but general, approach to produce an iterative algorithm solving the equation system in (23) is to decompose matrix A as the difference of two other matrices, M and N, as in(24)A=M-N. Then,(25)M-Nx=bMxk+1=Nxk+bxk+1=M-1Nxk+M-1b=xk+M-1b-Axk.This iterative scheme has the interesting property that if the residuals are b-Axk=0 (i.e., we have successfully found a solution of the equation system), then there is no update of the current solution.(i) Jacobi. A is decomposed in its diagonal matrix D and two lower L and upper U triangular matrices:(26)A=D-L+U⟹xk+1=xk+D-1b-Axk.(ii) Gauss-Seidel. A is decomposed in a similar way to the Jacobi decomposition above, but the L matrix is used differently:(27)A=L+D-U⟹xk+1=xk+L+D-1b-Axk.(iii) Richardson. A is decomposed using the identity matrix as reference:(28)A=I-I-A⟹xk+1=xk+b-Axk.The reasons why these numerical schemes succeed in solving the equation systemAx=b can be nicely illustrated in Richardson’s scheme. We may rewrite the iterative algorithm as(29)xk+1=b+I-Axk. We start with x0=0. Then, the estimates of x would be(30)x1=bx2=b+I-Abx3=b+I-Ab+I-A2b⋮xk=∑m=0k-1I-Amb.If A is an invertible matrix, then ∑m=0∞I-Am converges to A-1 as long as for all eigenvalues of the matrix A, λi, it is fulfilled that 0<1-λi<1, which for real eigenvalues translates into 0<λi<2 [82, Section 3.5.1]. Actually, this family of iterative algorithms is normally modified to introduce a relaxation factor μk (which may be different for each iteration):(31)xk+1=xk+μkM-1b-Axk. This relaxation factor helps to increase the radius of convergence of the algorithm by modifying the eigenvalues of the matrix involved in the convergence analysis (we have seen the analysis for Richardson’s iteration but different matrices are needed for the rest of the schemes). There are works analyzing the convergence of the reconstruction algorithm in terms of the properties of the sequence of numbers μk [83].SIRT is one of the most popular reconstruction algorithms used in Electron Tomography. As we show below, SIRT is the result of the Jacobi algorithm applied to the normal equations with a particular weighting scheme. Let us consider the normal equations of the Weighted Least Squares problem in real space:(32)ATWAx=ATWb. The Jacobi update with relaxation parameter is(33)xk+1=xk+μkD-1b-Axk. Let us decompose matrix A in its different rows (remember that the i-th row indicates how each basis function in the volume contributes to the i-th pixel),(34)A=a1Ta2T⋮aMP2T, and columns (the j-th column indicates how the j-th basis function affects all the pixels in the measurements),(35)A=α1α2⋯αB3. Then,(36)A=ATWA=∑i=1MP2wiaiaiTb=ATWb=∑i=1MP2wibiaiAxk=∑i=1MP2wiaiaiTxk=∑i=1MP2wi∑m=1B3aimxmkai.Let us now concentrate on a given basis function j:(37)djj=∑i=1MP2wiaij2=αjW2bj=∑i=1MP2wibiaij=αj,bWAxkj=∑i=1MP2wi∑m=1B3aimxmkaij=αj,AxkWxjk+1=xjk+μkdjj-1bj-Axkj=xjk+μkαj,b-AxkWαjW2.Interestingly, if the current image residual b-Axk is W-orthogonal to αj, that is, changing the coefficient of the j-th basis function does not affect the residual, then, as expected, the coefficient of the j-th basis function is not changed.In EM, we are used to formulate SIRT as(38)xk+1=xk+μk∑i=1MP2wibi-ai,xkai2ai.In the following, let us show that both formulations (see (37) and (38)) are not equivalent. Actually SIRT is not a single algorithm but a full family of reconstruction algorithms [84]; each one provides a different insight into the reconstruction process. If we look at a particular basis function j in (38), we have(39)xjk+1=xjk+μk∑i=1MP2wibi-ai,xkai2aij=xjk+μk∑i=1MP2wiai2bi-ai,xkaij.Let us now work on (37):(40)xjk+1=xjk+μkdjj-1bj-Axkj=xjk+μk1αjW2∑i=1MP2wibiaij-∑i=1MP2wi∑m=1B3aimxmkaij=xjk+μk1αjW2∑i=1MP2wibi-ai,xkaij.The algorithms in (39) and (40) are clearly different and, however, both are called SIRT and both belong to the SIRT family of algorithms. Any algorithm of the form(41)xjk+1=xjk+μk1γj∑i=1MP2wiρibi-ai,xkaijfor suitable γj and ρi numbers is also considered a SIRT algorithm [84]. A particular class for which convergence has been proven [84] is for(42)γj=∑i=1MP2wiaijαρi=∑j=1B3aij2-α0<α≤2,0<μk<2.The Jacobi iteration resulting in (40) corresponds to α=2; the typical EM SIRT (see (39)) corresponds to α=0 and is the one implemented in Xmipp [85, 86] and TomoJ [87]. However, ASTRA toolbox [88], also used in EM, has a SIRT algorithm with α=1.The case whereα=0 (see (38)) makes an interesting connection to the theory of feasibility problems [89–91]. Consider the hyperplane defined by all volumes compatible with the i-th experimental measurement:(43)Hi=x∈RB3∣ai,x=bi. Solving the tomographic problem amounts to finding a volume x such that it is compatible with all measurements (note that, in this formulation, we are disregarding the effect of noise which makes such an intersection be empty):(44)x∈H1∩H2∩⋯∩HMP2.Given a volumexk, we may orthogonally project it onto the hyperplane given by the i-th measurement by(45)ProjHixk=xk+bi-ai,xkai2ai.In this way, we may rewrite the SIRT iteration withα=0 as(46)xk+1=xk+μk∑i=1MP2wiProjHixk-xk.That is, at every iteration, we update the volume with a weighted sum of the orthogonal projections of the current solution onto the set of hyperplanes defined by the experimental measurements. The relaxation factorμk (normally chosen between 0 and 1, although there are convergence theorems for values between 0 and 2) could be understood as how much we go from our current position x(k) to the desired position ProjHixk. In a noiseless case, we are rather certain about the update and may set μk=1. In a noisy environment, we may be more conservative and use a smaller relaxation factor reflecting our distrust in the experimental measurements.Actually, we may update our estimate of the current solution after a single hyperplane projection; that is, we do not have to wait to “see” all measurements at the same time but we may update the volume just after seeing each pixel value:(47)xk+1=xk+μkwikProjHikxk-xk.The indexi(k) will go over all the experimental measurements. This scheme is known as ART (algebraic reconstruction technique) and the difference between SIRT and ART is the same as the difference between a Jacobi update (SIRT) in the solution of a linear equation system and a Gauss-Seidel update (ART). In the ART case, we update the volume as soon as we have new information available, while in the SIRT case, we update the volume with a consensus of all the information available. This property of ART can be exploited to use ART in a streaming mode as data is acquired and stop data acquisition as soon as the reconstruction achieves a specific criterion (this could minimize beam damage, for instance) [92]. ART converges much faster than SIRT, although it tends to produce little bit noisier reconstructions. This problem can be alleviated through the choice of a relaxation factor that decreases over time. A trade-off between updating after seeing every pixel (ART) and updating after seeing all pixels (SIRT) is given by Block-ART or Simultaneous ART (SART). The volume is updated after seeing a small set of pixels Sk (normally all those in the same experimental image):(48)xk+1=xk+μk∑i∈SkwiProjHixk-xk.Kunz and Frangakis [73] studied how the order in which these blocks were chosen affected the reconstruction’s quality.Additionally, we do not need to be restricted to orthogonal projections, and oblique projections can be undertaken. Given a symmetric, positive definite matrixG, the oblique projection onto the hyperplane Hi is defined as [93](49)ProjHiGxk=xk+bi-ai,xkaiG-12G-1ai.The study of these alternatives has given rise to a whole family of iterative algorithms (Cimmino, Component Averaging (CAV), Block-iterative CAV (BiCAV), Block-Simplified SART, iterative algorithms with Bregman projections, Block-iterative Underrelaxed Entropy Projections, Averaging Strings, etc.) [93–98]. A different variant of ART is a multiplicative ART (MART), where the iterative step is multiplicative instead of additive. Although some of these algorithms have been tested in Electron Microscopy [40, 99], none of these more advanced variants have made their way into a massive adoption. For extensive reviews of the algebraic reconstruction techniques, the reader is referred to Gordon and Herman [100], Gordon [101], and Herman [102].All SIRT algorithms can be written in a compact matrix form:(50)xk+1=xk+μkCATRb-Axk,with A being the system matrix (A projects the current estimate of the solution onto the image projection space, while AT back-projects the residual into the volume space) and R and C being suitable diagonal matrices (C acting on the basis function space cjj=1/γj; and R acting on the experimental measurements space rii=wi/ρi).Stated in this matrix form, we may easily find another SIRT algorithm that is very popular in image processing: the Landweber iteration. A possible solution to the Weighted Least Squares problem,(51)x=arg⁡minx⁡Ex=12Ax-bW2=arg⁡minx⁡12Ax-bTWAx-b=arg⁡minx⁡12bW2+12xTATWAx-xTATWb, is to use gradient descent iterations:(52)xk+1=xk-μk∇Exk=xk-μkATWAx-ATWb=xk+μkATWb-Ax. This latter iteration is called Landweber iteration and it can be easily recognized that it fits into the generic SIRT form of (50).Another interesting variant of this family of traditional algorithms is the possibility to use unmatched projectors, for example, a relatively complicated forward operator,A, including the projection geometry and multiple defoci effects and a simple backward operator, B, considering only the projection geometry [103]:(53)xk+1=xk+μkCBTRb-Ax. The convergence of this kind of algorithms depends on the eigenvalues of the matrix I-CBTRA and the specific details can be found in [103]. In EM, this strategy has been implemented in Xmipp (the forward projection considers geometry and CTF, but the backward projection only considers geometry). Although it was shown to successfully recover the underlying 3D structure, this method has not been further pursued because of the general problem of most iterative algorithms: their computational speed compared to direct Fourier inversion algorithms.In recent years, one of the most popular algorithms for solving a linear equation system is Conjugate Gradient [59, 104–106]. The success of this algorithm is demonstrated in its ability to converge to the local solution in B3 steps (as many steps as the dimension of the solution being sought). Note that this number is finite and much smaller than the number of steps of ART or SIRT, whose convergence theorems are on the limit as the number of iterations goes to infinity. Given the equation system(54)Ax=b, the trick is to use a set of directions pk which is orthogonal with respect to the inner product induced by A (k≠k′⇒pk,pk′A=0; these directions are said to be conjugate with respect to A). Then the solution sought can be written in the form(55)xk=∑k′=0kpk′,bpk′A2pk′. The set of conjugate directions is also iteratively constructed as(56)pk+1=b-Axk-∑k′=0kpk′,b-Axkpk′A2pk′. Actually, these summations can be nicely reorganized so that in real implementation we only need to keep a vector for the current solution, a vector with the residuals, and a vector with the current conjugate direction [104].In the context of EM, the Conjugate Gradient was used by Chen and Förster [59] in Fourier space. However, there are a number of variants to the basic Conjugate Gradient algorithm (Conjugate Residuals, Biconjugate Gradient, Stabilized Biconjugate Gradient, Lanczos method, Generalized Minimal Residuals (GMINRES), Bi-Lanczos, Conjugate Gradient Squared, Quasi Minimal Residuals, etc.; most of them belong to a family of algorithms called Krylov subspace algorithms [107, 108]), none of which have been tested in EM.In practice, none of these iterative algorithms are run to convergence. Instead, the algorithms are typically run for a fixed number of iterations (typicallyM or 2M in the case of Block ART, 20 for CG, and 100–150 in the case of SIRT). How deep these algorithms have gone in the objective function landscape depends on the conditioning number of the equation system:(57)κA=σmaxσmin, where σmax and σmin are the maximum and minimum singular values of matrix A. If this ratio is close to 1, iterative algorithms will quickly converge to the minimum of the error function. If the ratio is much larger than 1, then the problem is said to be ill-posed (small perturbations in the b vector may translate into large variations in the solution vector x) and the convergence speed of iterative algorithms is slow. In case of solving normal equations, we have(58)κA=κATA=κ2A. That is, the ill-posed character of matrix A is worsened. For this reason, in the theory of linear equation solving, it is customary to use a preconditioning matrix P so that we do not solve the problem(59)Ax=b but(60)P-1Ax=P-1b.The preconditioningP matrix is chosen in such a way that(61)κP-1A<κA. Although several preconditioners have been tested in different kinds of tomographic setups (computerized tomography [109], optical diffusion tomography [110], positron emission tomography [111], electrical capacitance tomography [112], acoustic waveform tomography [113], etc.), in EM, there has not been any attempt to use any preconditioning, although this is an issue the community is aware of [114].


## 5. Constrained 3D Reconstruction

The feasibility problem of (44) can be complemented with some additional constraints representing a priori knowledge about the reconstruction. We may do so by imposing the fact that the volume also belongs to some convex set (a set is convex if, for any two volumes in this set, x1 and x2, the linear combination 1-λx1+λx2, with 0≤λ≤1, also belongs to that set; the set of all nonnegative volumes is convex as well as the set of all volumes defined within a mask, the set of all volumes with a given symmetry, the set of all volumes bandlimited to a given frequency, etc. [64]). Given a collection of convex sets representing our a priori knowledge about the reconstructed particle, C1,C2,…,CK, we may try to find a feasible volume at the intersection(62)x∈H1∩H2∩⋯∩HMP2∩C1∩C2∩⋯∩CK. This problem has been extensively explored in the EM community under the name Projection Onto Convex Sets (POCS) and the interested reader is referred to Carazo and Carrascosa [115, 116], Carazo [117], García et al. [31], Sorzano et al. [64], and Deng et al. [118]. The idea is to alternate between projections onto the subspaces defined by the experimental measurements (Hi) and projections onto the convex sets (the projector onto the set of nonnegative volumes simply sets all negative values to 0; the projector onto the set of volumes defined within a mask applies that mask to the current solution; the projector onto the set of symmetric volumes symmetrizes the current solution; the projector onto the set of bandlimited volumes applies a low-pass filter to the current solution; etc.). In the signal processing community, POCS algorithms have been generalized by algorithms using the so-called Proximity operators. Despite being out of the scope of this review, because these algorithms have not been introduced in EM, the interested reader may follow the review of Combettes and Pesquet [119] and the references therein.Landweber iterations can also be set in a constrained setup. Let us assume that we have the a priori knowledge thatx is in a convex set C. Then, the solution of the constrained problem(63)x=arg⁡minx∈C⁡12Ax-bW2 can be found with the iterative algorithm [119]:(64)xk+1=ProjCxk+μkATWb-Ax.As seen in the two examples above, one of the most interesting ideas of this constrained optimization is the possibility to alternate between the standard tomographic update (x(k)+μkATWb-Ax) and the projection onto convex sets (POCS, ProjC).These ideas of constrained optimization can be further extended to nondifferentiable, convex functions (thel1 norm of the residual is a function of this kind). Let us assume that we are minimizing a nondifferentiable, convex function ϕx:(65)x=arg⁡min⁡ϕx.At a differentiable point ofϕ, a gradient descent iteration is perfectly suitable to find the minimum of the objective function:(66)xk+1=xk-μk∇ϕxk.The problem at a nondifferentiable point (which occurs normally at the frontiers of the intersection of convex sets) is that the gradient is not well defined. We may define instead the subgradient. A vectorgk∈RB3 is a subgradient of ϕ at xk if for any x we have(67)ϕx-ϕxk≥gk,x-xk. Intuitively, the gradient defines at differentiable points a unique hyperplane that is tangent to ϕ and, since ϕ is convex, this tangent plane is always below ϕ. At nondifferentiable points, there are an infinite number of hyperplanes touching at xk,ϕxk and below ϕ. The normal vectors to these hyperplanes constitute the set of subgradients. The subgradient method iteration is then, at these nondifferentiable points [120, 121],(68)xk+1=xk-μkgk. We may easily add convex constraints to this method obtaining constrained subgradient minimization (known as Projected Subgradient minimization) [122, 123]:(69)xk+1=ProjCxk-μkgk.Superiorization has been proposed as an alternative to Projected Subgradients [43, 124–126]. The idea of Superiorization is to steer the reconstruction iterations towards the constraints. One of the main differences between Projected Subgradient and Superiorization is that if there are several convex constraints, Projected Subgradient requires the projection onto the intersection of all of them, the feasible set, while Superiorization requires the sequential projection onto each of the constraints. Note that the intersection may be empty. For instance, the set of all volumes defined within a finitely supported mask and bandlimited is empty. Also, the projector onto the intersection set may be more complicated than the projector onto each one of the individual constraints. This limits the applicability of Projected Subgradient, while this is not a problem for Superiorization. Although these algorithms are available in the tomography community and are well characterized in other domains, none of them have made their way into EM.A different approach to constrained optimization which has been actively explored in EM is by defining new equations that must be simultaneously solved along with the equation system coming from the measurements [64]. For instance, a mask can be easily expressed in terms of basis functions. Let us assume that there is no density coming from the molecule at a given location r0. Then, we may particularize the series expansion at r0 (see (5)) and add the linear equation(70)Vr~0=∑jxjbr~0-r~j=0. In this way, we handled a priori information coming from masks, symmetry, total molecular mass, and nonnegativity. Adding this a priori information was certainly relevant when the number of projections was small. However, as the number of projections increased, the improvement by the a priori information was less noticeable, showing that most of the information was already present in the experimental dataset. Now, with increasing pursuit of high-resolution results, there is again room for exploring these extra sources of information available, although for the moment this line of research has not been resumed.Iteratively steering algorithms tend to promote reconstructions with certain characteristics. For instance, [65, 71, 72, 127, 128] focused on the problem of reconstructing discrete valued objects (the reconstructed volume could only have a few, normally two (background and foreground), values) or objects with very few active voxels. They use projectors similar to those used in convex sets; for instance, bgART [71] defined(71)Projxjk=mkxjk<mk+Kskxjkxjk≥mk+Ksk, where mk is the estimated mean of the background gray values at iteration k, sk is their estimated standard deviation, and K is a user-selected multiple (normally a number between 3 and 6). The main difference between these iterative steering algorithms and those projecting onto convex sets is that the former project onto a nonconvex set; additionally, the set onto which the reconstruction is projected changes from iteration to iteration. The convergence of this kind of algorithms was studied in [129]. Reference [75] presented an algorithm for Electron Tomography in which the steering is driven by a nonlinear diffusion denoising algorithm. After each reconstruction step, the reconstructed volume is denoised by applying a step of a nonlinear diffusion algorithm (the projector). For these steering algorithms to converge, the reconstruction algorithm must be “perturbation resilient” [124] and the perturbation must be “small enough” as not to destroy the work of the reconstruction algorithm.


## 6. Sparse 3D Reconstructions

Sparse representations has been one of the most active research fields in the area of image and signal processing in the last 10–15 years [130–132]. The idea is that natural images and objects have a representation in some appropriate space in which very few nonnull coefficients are needed. This space may be fixed (e.g., wavelet transform or discrete cosine transform and DCT) or computed ad hoc for a particular problem (dictionary based algorithms). Knowing that our object has a sparse representation helps the algorithm to concentrate the energy in a few coefficients, preventing the energy dispersion normally caused by noise. In this problem setup, vector norms different from the Euclidean are normally employed. In general, the lp norm is used:(72)xp=∑jxjp1/p. Euclidean norm is obtained for p=2, Manhattan norm is obtained for p=1, and for p=0 the norm of a vector is simply the count of the number of nonzero coefficients in the vector (technically, l0 is not a norm because it does not fulfill the condition λx0=λx0). The goal of sparsity is to find a representation of x such that(73)x=argminx0s.t.Dx=x, where D is the dictionary of elements available to represent x (the wavelet, DCT, or ad hoc dictionary) and x is the representation of x in that dictionary. The aim is to have as few coefficients different from 0 as possible.As stated above, sparse representations are mostly interesting inl0-norms. However, having this norm in the objective function requires combinatorial optimization techniques, known to be NP-hard in computational complexity. Interestingly, lp norms with 0<p<2 tend to promote sparse representations (x has relatively few nonzero coefficients) and efficient algorithms have been developed in the recent years for these minimizations.The first class of algorithms we will review involve anl0-norm and are those solving any of these two related problems:(74)x=arg⁡minx0s.t.b-ADx22≤ϵ(75)orx=arg⁡minb-ADx22s.t.x0≤N.Note that our “tomographic dictionary” (AD=ATAD) includes the standard dictionary representation, D, as well as information about the projection structure of the tomographic setup A=ATA. In the following, let us refer to this “tomographic dictionary” by D and to its j-th column by dj (in the dictionary jargon, this would be the j-th atom and we will assume that there are J≪B3 atoms in the dictionary). Many different algorithms exist to solve this problem [130]; none of them have been tested in EM. Among them, Matching Pursuit [133] and Orthogonal Matching Pursuit [134] are two of the most popular algorithms (actually the latter was tested in a tomography setting by researchers working on EM, although it was not applied to EM data [135]). For its simplicity, let us show the iterations in Matching Pursuit to illustrate the flavor of these algorithms and how they differ from the iterative algorithms presented so far.Letxk represent our current reconstruction. In the first iteration, it will be x0=0. We now choose the atom that is maximally aligned with the residual of our equation system, r(k) (note that r(0)=b), calculate its coefficient, and update the residual. This process is iterated until the desired number of atoms is reached. The following algorithm is executed from k=0 to k=N-1.(i) Step 1. Seek the atom maximally aligned with this residual(76)j=arg⁡maxj′⁡rk,dj′.(ii) Step 2. Update the current solution with the projection of the residual onto this atom(77)xj′k+1=xj′kj′≠jrk,djdj2j′=j.(iii) Step 3. Update the residual(78)rk+1=rk-xjk+1dj.This algorithm is greedy and finds a suboptimal solution of the reconstruction problem (see (75)) [133]. However, it is very easy to implement. A significant improvement is provided by Orthogonal Matching Pursuit in which Step 2 is modified to update all coefficients (not just one) by orthogonally projecting b onto the subspace spanned by all the atoms employed so far.The second class of algorithms substitute thel0-norm by an l1-norm: (79)x=arg⁡minx1s.t.b-ADx22≤ϵ(80)orx=arg⁡minb-ADx22s.t.x1≤δ.This allows more efficient optimization algorithms to be employed. Classical algorithms are Least Absolute Shrinkage and Selection Operator (Lasso also known as Basis Pursuit), Iterative Reweighted Least Squares (IRLS), Iterative Shrinkage algorithms, Least Angle Regression (LARS), and any of their variants [130]. For their simplicity, Iterative Shrinkage algorithms have found their way into Electron Microscopy [136, 137] (for a theoretical background, you may also see [138, 139]). However, the EM versions are aimed at solving a regularized problem,(81)x=arg⁡min⁡b-ADx22+λx1, rather than the constrained problems above. In their most simplified formulation, these algorithms alternate between a standard reconstruction update (let us, for instance, take a SIRT step) and a soft-thresholding step in some suitable, sparse space (e.g., the wavelet space). The rationale is relatively simple and resembles the line of thought of the iterative steering algorithms. After applying a standard reconstruction step, the solution is steered towards a sparse solution by setting to 0 all small coefficients in some space known to promote sparsity (as the wavelet space). A basic iterative step would be(82)xk+1/2=xk+μkCATRb-Axkxk+1=∑j=1JSλxk+1/2,djdj,where Sλ(x) is the soft-thresholding function:(83)Sλx=x+λ2x<-λ20x≤λ2x-λ2x>λ2.In the above iteration, we have assumed an orthonormal dictionary or transformation, but equivalent formulas can be found for nonorthonormal dictionaries.The use of Lagrangian augmented objective functions of the form(84)x=arg⁡minx⁡Fb-Ax+λGx is a widespread technique in data analysis. The first term, F, is called the data fidelity term, while the second term, G, is a regularizer. Many different schemes respond to this structure. For instance, a typical regularized Weighted Least Squares problem would be given by Tikhonov quadratic regularization (also known as ridge regression) involving some “preprocessing” matrix L,(85)x=arg⁡minx⁡b-AxW2+λLx22, whose normal equations would be(86)ATWA+λLTLx=ATWb. Assuming Gaussian errors for the measurements and a Gaussian prior for the volume coefficients, we would have L=Σx-1/2.(87)ATWA+λΣx-1x=ATWb. If we think that x and b are in Fourier space, these would be the normal equations associated with the same Bayesian problem that Relion is solving [7]. W contains the CTF information, A the data collection geometry (including the probability of each projection having a given projection direction and shift), and Σx our prior about the energy of the coefficients in Fourier space. A difference between Relion and these normal equations is that Relion reestimates the prior after each iteration.We might have gone one step further to the generalized Tikhonov regularization,(88)x=arg⁡minx⁡b-AxW2+λLx-x0Q2, whose normal equations are(89)ATWA+λLTQLx=ATWb+QLx0.If we assume independence between the preprocessed coefficients, most maximum a posteriori (MAP) algorithms could be written using(90)GLx=∑jϕLxj, where ϕx is the a priori, negative log likelihood of observing a value of x (ϕx=-log⁡pxx). Assume that a Gaussian distribution of the x coefficients results in L=Σx-1/2 and ϕx=x2+C (C is a constant that does not affect the MAP optimization). Assume that a centered Laplace distribution of standard deviation σ (pxx=1/2σe-2|x|/2σ) results in L=2/σI and ϕx=x+C with G being an l1-norm (instead of an l2-norm as in the case of the Gaussian prior). In an EM setup, Moriya et al. [74] assumed a Median Root Prior which favors locally monotonic reconstructions.If we takeL to be a volume derivative operator and take G to be the l1-norm of this derivative (that is, a Laplace prior on the derivative coefficients),(91)Gx=Lx1, then we have a total variation regularization (also named as TVL1). The problem with the l1-norm is that it cannot be differentiated and sometimes it is substituted by(92)Gx=∫R3∂Vr∂x2+∂Vr∂y2+∂Vr∂z2dr, where V(r) is the volume reconstructed using the x coefficients (see (5)). This was explored for EM by Zhu et al. [62], Aganj et al. [63], Li et al. [66], Goris et al. [57], and Zhuge et al. [76].The approach of Albarqouni et al. [140] also falls into this regularization category. The function G is, in their case, a function borrowed from robust statistics, the Huber function:(93)ϕτx=0.5x2x≤ττx-0.5τ2x>τ. This function is half-way between a quadratic and a linear penalization (τ controls the switch between these two behaviors). For low values of the derivatives, the function behaves as a quadratic term and for high values it behaves as a linear term. The idea is not to let large derivatives dominate the optimization process.We can see that the above regularized problems include either some a priori knowledge on the volume (x) or a property of that volume (Lx, normally its derivative). However, we could include both through a new algorithm (alternating-direction method of multipliers, ADMM) recently introduced for EM [141]:(94)x=arg⁡minb-Ax22+λ1G1x+λ2G2us.t.u=Lx.The actual problem being solved is an augmented Lagrangian:(95)x=arg⁡minx,u,α⁡Lμx,u,α=b-Ax22+λ1G1x+λ2G2u+αTLx-u+μ2Lx-u2, where α is the set of Lagrangian multipliers. The ADMM proceeds iteratively as(96)xk+1=arg⁡minx⁡Lμx,uk,αkuk+1=arg⁡minu⁡Lμxk+1,u,αkαk+1=αk+μLxk+1-uk+1.Related to these sparse reconstruction problems is the one of compressed sensing. The idea is to perform a 2D or 3D reconstruction problem starting, not from a full image, but from a “few” incoherent points from the projection images [142–144]. The trick is that the a priori knowledge that the solution is sparse allows reconstructing the full image or volume from fewer points than Nyquist theorem requires (there is a lower limit on the number of measurements required). This theory has not been explored in biological EM applications, but it has been studied in material science EM [145–150], especially when Scanning TEM (STEM) is used [151]. For biological samples, Energy Filtered TEM (EFTEM) is a clear candidate to benefit from a compressed sensing acquisition [149].


## 7. Conclusions

The field of iterative reconstruction algorithms has been very much studied, particularly in its application to Electron Microscopy data, as we have shown in this review. The 3D reconstruction problem is no longer seen as a bottleneck in Single Particle Analysis (this is a technique in which many single particles, assumed to come from an homogeneous population but at different angular orientations, are combined into a single 3D map; Jonic et al. [47]), and the few attempts to complement the data with a priori knowledge are not in widespread use, probably due to their higher computational cost. However, as we approach resolutions of 2-3 Å, it may be worthy to retake this line of research as a way to increase the reconstruction resolution. In Electron Tomography, the situation is different due to the reconstruction artifacts induced by the maximum tilt angle limitation, the presence of gold beads, and the low number of projections. Additionally, the particular data collection geometry (normally, single tilt axis) favors the adoption of a pure 2D reconstruction approach that reduces the reconstruction problem in one dimension, implying a great reduction in the computational cost. In this field, iterative algorithms capable of incorporating a priori information are still a very active field of research. Traditional iterative algorithms as ART or SIRT still dominate the “market.” However, very powerful reconstruction algorithms with more modern approaches incorporating convex, nonconvex, and sparsity constraints are continuously appearing and most likely, in the near future, one of these algorithms will eventually become the standard de facto.



---

*Source: 6482567-2017-09-17.xml*