Using Convex Optimization of Autocorrelation with Constrained Support and Windowing for Improved Phase Retrieval Accuracy

Using Con v ex Optimization of Auto correlation with Constrained Supp ort and Windo wing for Impro v ed Phase Retriev al Accuracy Pietrini, Alb erto alberto.pietrini@icm.uu.se Nettelblad, Carl carl.nettelblad@it.uu.se No vem ber 6, 2021 Abstract In imaging mo dalities recording diffraction data, the original image can be reconstructed assuming kno wn phases. When phases are unkno wn, o versampling and a constrain t on the support region in the original ob ject can b e used to solv e a non-con vex optimization problem. Suc h schemes are ill-suited to find the optim um solution for sparse data, since the recorded image does not corresp ond exactly to the orig- inal w av e function. W e construct a conv ex optimization problem using a relaxed supp ort constrain t and a maximum-lik elihoo d treatment of the recorded data as a sample from the underlying wa ve function. W e also stress the need to use relev ant windo wing techniques to account for the sampled pattern b eing finite. On simulated data, we demonstrate the b enefits of our approach in terms of visual quality and an improv ement in the crystallographic R- factor from .4 to .1 for highly noisy data. 1 In tro duction In imaging applications ranging from astronomy to molecular biology , diffraction- based imaging tec hniques hav e b een demonstrated to be useful. In general, the diffraction from an optically thin ob ject in the far field at low angles can b e appro ximated as the F ourier transform of the scattering p o wer of a pro jection of the original ob ject. Th us, to reconstruct this pro jection, the in verse F ourier transform should b e applied to the collected image. How ev er, imaging detec- tors tend to only record light intensit y (amplitude squared), losing the complex phase of the underlying wa v e. In crystallographic applications, common in imaging of biological particles lik e proteins, this phase r etrieval problem b ecomes under-constrained, based on the physical data alone. The so-called Bragg peaks representing constructive in terference betw een all unit cell rep eats will not, in general, constrain the 1 con tent of the unit cell to a single represen tation, i.e. a single phase assignment. By introducing knowledge of p ossible biological structures in different w ays, the problem still b ecomes tractable [1]. F or isolated particles, phase retriev al based on o v ersampling relative to a supp ort constraint has b ecome the metho d of c hoice. The resulting equation system is non-linear, and solving it in a minim um residual sense results in a non- con vex optimization problem. Th us, most solution metho ds do not mak e claim on pro viding global optima, but try to balance finding optima and exploring a reasonable portion of the domain. P opular solution metho ds use v ariations of alternating pro jection techniques, such as Hybrid-Input Output (HIO, [2]) and Relaxed Averaged Alternating Reflections (RAAR, [3]). A concise comparison of v arious iterativ e sc hemes and their similarities can be found in [4]. Sev eral of these iterative methods were dev elop ed for conditions of relativ ely high signal, with no or only mo derate noise. This is very different from a sparse imaging regime with pixel detectors able to record individual photons, whic h is no w the reality in e.g. FXI (Flash X-Ra y Imaging), where individual biologi- cal particles are imaged in fem tosecond-duration X-ra y pulses delivered using X-ra y free electron lasers, suc h as the Linac Coherent Ligh t Source (LCLS [5]). In addition, some parts of the image migh t b e missing due to detector geom- etry , detector malfunction, or the geometry of the experiment itself. Despite these limitations, successful reconstructions hav e b een rep orted for organic and inorganic samples alike [6, 7, 8, 9]. In this communication, w e in tro duce a relaxed v ersion of the supp ort con- strain t, based on the known relationship b et ween the F ourier transform of the original wa v e, and the F ourier transform of the observed quan tity , the intensi- ties. W e sho w the con vex c haracter of the optimization problem of estimating maxim um-likelihoo d intensities giv en a recorded sparse sampling of the diffrac- tion pattern, a supp ort, and a noise mo del. Con vex optimization problems lac k separate lo cal optima and allow a wealth of existing methodology to b e used to identify p oin ts arbitrarily close to the true global optimum. W e name this method Con vex Optimization of Auto correlation with Constrained Supp ort (CO ACS). CO ACS has similarities to methods suc h as PhaseLift [10, 11] and PhaseCut [12]. How ever, we deliberately disentangle finding optimum in tensities from find- ing the phases. This, together with great care in the n umerical implementation, seems to partially o vercome the n umerical stabilit y concerns that led the au- thors in [11] to conclude that not only would PhaseLift b e unsuitable for phase retriev al of o versampled particles with known supp ort, but that the problem itself is intractable to n umerical solution without adding further constrain ts. W e demonstrate a pro of of concept solution metho d using a pre-existing soft ware pac k age for compressed sensing, TF OCS (T emplates for First-Order Conic Solvers) [13]. W e ev aluate phase retriev al with and without COA CS pre- pro cessing based on TFOCS, realizing mark ed impro v ements. In the remaining sections, w e first describe our optimization problem in greater detail, follow ed b y a detailed accoun t of implementation concerns, ex- p erimen ts, results, and their ev aluation. 2 2 In tensit y healing using CO A CS W e will fo cus on the case of 2D imaging, since it is straightforw ard, yet rich enough to illustrate most of our p oin ts. The ov erall methodology generalizes to 3D, lik e the iterativ e phase retriev al used as the final step in [14], and for that matter to further dimensions. Assume that the 2D pro jection b eing imaged is P . Assuming ideal conditions (a plane incident wa ve, a thin ob ject, a far-field regime, and only considering lo w scattering angles), the scattered wa ve X meeting a square detector satisfies the following relation: X ∝ F ( P ) (1) Since actual photon detectors use discrete pixels cov ering a limited region in space, we will treat a discretized, finite case, where P , X (matrices in b old) are 2D center-cropped discretizations of their contin uous counterparts. F or practical purp oses, we can then consider the following version of eqn. 1: X ∝ F ( P ) (2) Ho wev er, these tw o formulations are not identical. The discrete F ourier transform assumes that the underlying ob ject is perio dic. Figure 1 illustrates the effects of discretizing the infinite ideal diffraction pattern of our sim ulated test particle, with small but non-negligible artifacts pe rmeating to the edge of the image. The physical pro cesses that we mo del are based on eq. 1. T o handle those prop erly in the discretized space, we extend eq. 2 with the Hann window, kno wn from sp ectral tec hniques in signal pro cessing. In 2D, the Hann window in real space amounts to a conv olution with a 3x3 blurring stencil, and in the F ourier space a translated sin usoid curv e dropping off to 0 at the edges. This ensures a representation where the high-frequency conten t smo othly go es to zero. Ha ving no explicit window is equiv alen t to a sharp rectangular windo w, giving rise to the sp ectral leak age effects seen in Figure 1. If we can assume a photon-counting detector with uniform quan tum effi- ciency r , we will then ha ve a sampled diffraction pattern in terms of an in teger matrix B , where: B i,j = Po ( r X i,j X i,j ) = Po ( r | X i,j | 2 ) (3) where Po is a Poisson distribution with the rate parameter λ of appropriate dimension iden tifying the mean (and v ariance). F or high rates, B can be treated as a Gaussian, or ev en as an exact observ ation B i,j = r | X i,j | 2 . Ev en a strong diffraction pattern will tend to hav e a characteristic structure of minima, similar to those from the most trivial slit exp erimen t. Th us, ev en though the pattern as a whole is strong, the Gaussian or identit y approximation might not hold for all pixels within the pattern. 2.1 Supp ort constraints Since a single particle is imaged in isolation, we can in tro duce a supp ort con- strain t, where P S { , or P S { , is 0 (where { signifies the complement). F or phase 3 a) b) c) d) 0.0 1.5 Photons e) 0.00 0.05 A.U. f) 0.0002 0.0000 0.0002 A.U. Figure 1: Sim ulated particle with and without Hann windo wing. a) Original sim ulated diffraction pattern with high-density sphere. b) Center crop of dis- crete F ourier transform of a), sho wing the particle. c) The image in b) with limited range to sho wcase artifacts outside of the ob ject outline. d) The sim u- lated pattern with Hann windo wing applied, resulting in low er in tensity aw ay from the center. e) Discrete F ourier transform of pattern after windo wing. Sligh t blur visible. f ) The image in e) with limited range. Artifacts found in c) mostly absen t. 4 retriev al, one can thus express this as: arg max { X | P S { =0 } P ( X | B ) (4) Alternating-pro jection metho ds are very commonly used for phase retriev al. They w ere pioneered by Fien up [2], and are still popular in practical applica- tions use, as e.g. implemented in [15]. The name comes from the fact that they alternatingly implement some v ariation on the support constrain t, and the in- tensit y (“F ourier”) constrain t. The standard implemen tation is then to assume that the in tensities are exact, rather than sto c hastic as giv en abov e. Both of these constrain ts can b e v aried, to b etter explore the full solution space and accelerate conv ergence [2, 3], or in order to accoun t for measurement error in the in tensities [16]. More recent attempts include using the Alternating Direc- tion Metho d of Multipliers for easily enforcing more general constraints in b oth spaces [17]. These attempts share the prop ert y that they do not claim global conv er- gence. Recent adv ances apply semi-definite and conv ex programming approaches to the full phasing problem, finding low-rank factorizations for matrices in or- der to pro duce the phases, such as PhaseLift [10] and PhaseCut [12]. These approac hes claim to pro duce unique solutions under fav orable conditions with a v anishing probability of not doing so, suc h as observing an ob ject using alternat- ing binary masks, or from unique angles randomly distributed o ver a sphere [10]. While suc h conditions are p ossible to ac hieve in some situations, they would be c hallenging to transfer to the single particle X-ray diffraction experiments. PhaseLift considers the rank-one matrix pp , where p w ould be the flattened v ector from our 2D matrix P . A megapixel image would th us give rise to a trillion-elemen t matrix. It is possible to formulate an alternate op erator based on the F ourier transform that goes from such representations to the diffraction pattern. The problem can then b e solved in terms of satisfying the observ ed data with a suitable matrix in this space, while also minimizing the rank. A con vex represen tation of that problem can be expressed as minimizing the trace. The authors of PhaseLift note that while ov ersampling in theory giv es a unique solution to the phasing problem, uniqueness is not equiv alent to practical n umerical attainability . In fact, they ev en claim, based on simulations including their own metho d as well as alternating pro jections [11] (original italics): The ill-p osedness of the problem is eviden t from the disconnect b et w een small residual error and large reconstruction error; that is to sa y , w e fit the data v ery w ell and y et observe a large reconstruction error. Thus, in stark c ontr ast to what is widely b elieve d, our simula- tions indic ate that oversampling by itself is not a viable str ate gy for phase r etrieval even for nonne gative, r e al-value d images. This argumen t is supp orted b y additional claims about numerical stabilit y . W e will argue that a new analysis, of the supp ort constraint sp ecifically , will allo w a muc h more fav orable n umerical treatment. 5 2.2 A supp ort constrain t for the P atterson function The quan tity that is sampled to form the directly observ able B is Y = X  X (the element wise Hadamard pro duct with the conjugate, forming the absolute v alue squared). In crystallographic applications, the F ourier transform of Y is referred to as the P atterson function. Using the conv olution theorem, we can c haracterize this as: F ( Y ) = F ( X  X ) = P ∗ P (5) Here, ∗ indicates conv olution. The existing supp ort constraint (4) for X can then b e form ulated as a (w eaker) supp ort constraint for Y : arg max Y , F ( P ∗ P ) ( S ∗ S ) { =0 P ( Y | B ) (6) arg max Y , F ( ˆ P ) ˆ S { =0 P ( Y | B ) (7) Based on the structure of the probabilit y distribution function in (3), (7) is a lik eliho o d-optimization problem with a linear transform in to indep enden t v ariables. The previous quadratic term for X has v anished, at the p ossible exp ense of the self-conv oluted v ersion of S , ˆ S , imposing a less strict constrain t on the solution. The self-conv olution effect of treating the intensit y directly has also led some authors to refer to this as the auto correlation function. 3 Implemen tation If the noise structure would ha ve been Gaussian, the resulting problem from (7) w ould hav e amoun ted to an ordinary least squares problem. A solution could b e determined b y solving the o ver-determined homogeneous linear equation sys- tem based on the F ourier op erator for those pixels that are in the complemen t of the self-con voluted supp ort. F or a practical solution, one will how ever also need to add the constrain t of all intensities being non-negative, which necessi- tates sp ecialized conv ex optimization solvers, rather than the most basic least squares routines. Negative in tensities would b e imp ossible to transform bac k to amplitudes in the original phasing problem. F or a long time, so-called interior p oin t metho ds hav e b een a metho d of c hoice for many families of con vex optimization problems [18], with free [19] as w ell as commercial [20] solvers, and more high-lev el mo deling en vironments [21]. These metho ds explicitly form the Hessian of a problem b eing solv ed, with a n umber of elemen ts equiv alen t to the square of the num b er of v ariables. They tend to rely on sparseness in the dependencies betw een v ariables. Unfortunately , whereas the supp ort constrain t is one t yp e of sparsity , the F ourier op erator itself is dense. In some solv er implementations, the op erator w ould also ha ve to b e formed explicitly as a matrix, rather than the far more efficien t F ast F ourier T ransforms. The b enefit of in terior p oin t metho ds is that they tend to conv erge v ery quic kly , in terms of the num b er of iterations. 6 W e hav e instead chosen to rely on first-order methods, where only the gradi- en t is computed. The TF OCS [13] pac k age is a Matlab toolb o x, and is presen ted b y the author as a set of templates from whic h one can build tailored solv ers. There are a couple of different canonical forms for form ulating problems in TF OCS, and w e ha ve used the one that corresp onds to the main function call tfocs : φ ( X ) = f ( A ( X ) + b ) + h ( X ) (8) X opt . = arg min X φ ( X ) (9) The linear op erator A does not need to be implemented as a matrix, but can be expressed in code, assuming it satisfies certain sp ecific criteria. W e can th us rely on an actual efficient F ourier transform implemen tation. The function f needs to b e smo oth and differentiable, with an explicit implemen tation of the gradien t. h on the other hand, only has to b e prox-capable, meaning one can (quic kly) solv e: Φ h ( X , t ) = arg min Z h ( z ) + 1 2 t − 1 || Z − X || 2 2 (10) In our context, f corresponds to the noise mo del, while h represents the supp ort constraint. Both of these are implemen ted as custom functions, for n umerical accuracy reasons outlined in the next section. 4 Numerical concerns As noted b y the authors in [11], n umerical stability and accuracy concerns are of great imp ortance in the phase retriev al application, especially when relying solely on the ov ersampling of a single pattern. Although the p oint of doing a separate step of intensit y healing is to make the phase-retriev al problem more w ell-p osed, these n umerical concerns are of imp ortance for intensit y healing as w ell. The intensit y in a diffraction pattern of, say , a uniform sphere, will deca y in a manner roughly proportional to q 4 , where q is the scattering angle [22]. The dynamic range needed for a problem in intensit y space is higher than in the amplitude space, due to the quadratic relationship b et ween the tw o. Since the F ourier transform adds terms for all frequencies for each resulting pixel, and terms are of v arying sign, loss of precision due to terms of opposite sign canceling out can b e paramount. In an y iterativ e approac h, the numerical accuracy needs not only hold for the final result, but for prop erly ev aluating the c hange undertak en betw een iterates. The steps outlined below try to consistently reduce the chance of: • T runcation error due to v ery small terms being added to large terms. • Loss of precision due to large negative v alues b eing added to p ositiv e v alues of almost equal magnitude, or vice versa. 7 • V ery small steps due to the optimal solution being lo cated at the border of the allow ed domain. The sp ecific places where this has b een taken in to accoun t can b e summa- rized as: • Lipschitz backtrac king structure in TF OCS • T ranslation of the solution v ariables and customized probability function • Relaxed quadratic barrier and accelerated contin uation • Window-a w are supp ort constraints It should b e noted that without these steps, one arriv es at a seemingly math- ematically equiv alen t problem formulation that still results in erratic numerical b eha vior in practice. The sp ecific nature of these steps also relate to design c hoices in the existing TFOCS solver, including the fact that differences b et ween successiv e iterates are computed, and that the ob jective function is handled as a scalar, rather than considering the individual p er-pixel con tributions. On the other hand, these properties are in no wa y unique to TF OCS. 4.1 Lipsc hitz bac ktrac king structure TF OCS can b e described as an alternating-pro jections method tracing mo dified gradien ts, giving far accelerated conv ergence compared to a fixed-step gradient descen t. F or more details on the o verall design of the pac k age, w e refer to [13, 23]. One crucial aspect of the algorithm is a backtrac king approac h to identify the prop er step size. The backtrac king is common to the v arious sp ecific solver sc hemes implemented in TFOCS. The bac ktracking is using a lo cal, adaptiv e estimate of a Lipsc hitz constan t for the gradien t of the ob jectiv e function. As elab orated up on in [13], there are multiple w ays to estimate a b ound on the Lisp c hitz constant based on lo cally ev aluated gradients a nd function v alues. Since the problem is inheren tly related to ev aluating the difference in function v alue b et ween t wo points, there is a great risk for loss of precision. The original TFOCS approach is therefore to use tw o bounds, one of which will giv e less conserv ative estimates, while also b eing more sensitiv e to the details of floating-p oint math, and another that is safer. The TF OCS co de automatically detects whether the relative difference in function v alue b et w een the t wo points is small, and then switc hes to the more conserv ativ e b ound. The conserv ativ e metho d will then b e used un til the next “restart” of the TF OCS algorithm. W e hav e made tw o mo difications to this logic. The first mo dification is to mov e the ev aluation point for c ho osing whic h bound form ulation to use. The original TFOCS co de p erforms one step and then, retrosp ectively , chec ks whether the difference w as small. With this design, a single step in the iteration pro cess might b e taken using an accuracy-deficient Lipschitz estimate. This 8 turned out to be enough to cause drastic divergence for us in some scenarios. W e instead put this test earlier. W e also allow the iteration pro cess to leav e the conserv ativ e regime, without a restart. F urthermore, w e consider not only the difference in function v alues | f ( x ) − f ( y ) | , but also the difference | x − y | . If the step itself is small relative to the magnitude of the v alues, truncation errors as w ell as loss of precision are possible. The total change is less than 20 lines of co de in a single function in the pac k age. It has b een submitted to the authors as a suggested change. This, together with using the “no-regress” restart feature giv es adequate optimization p erformance without div ergence or singularity issues. 4.2 T ranslation of the solution v ariables and customized probabilit y function The F ourier transform is a linear op erator. Thus, the relationship b etw een Y and ˆ P can b e decomp osed: Y = F ( ˆ P ) (11) Y = Y 0 + Y ? = F ( ˆ P 0 ) + F ( ˆ P ? ) (12) In this case, the optimization v ariable for the first-order method in TFOCS can b e Y ? . With Y 0 c hosen prop erly , such as the resulting solution from a previous run of the TF OCS algorithm (the previous “outer iteration”), the norm of Y ? should b e muc h smaller. This c hange reduces the effects of truncation errors in the gradient steps in individual iterations, including the probability of resorting to the conserv ative Lipschitz bound discussed abov e. It also reduces the effects of loss of precision due to large v alues of opp osite sign canceling out within the F ourier transform. The effects do still exist on a global lev el, but they do not hamp er the ev aluation of ob jective v alues, constraint satisfaction, and searc h direction lo c al ly during the iteration pro cess. The equiv alen t translation is not p ossible in the phasing problem, since an y c hanges in the phase estimates affect the full amplitudes, whereas an additive decomp osition is p ossible in the in tensity space. T o fully realize the b enefits of translation, the corresponding shift needs to b e done in the computation of the f term of the ob jectiv e function as well. W e hav e implemented a custom log Poisson implemen tation, which computes the change in ob jectiv e function v alue b et w een Y 0 and Y using Y ? as in- put. T o further accelerate conv ergence, esp ecially for pixels with zero observ ed photons (where the log-Poisson optim um lies exactly at 0, the b order of the non-negativit y constrain t), we replace the log-P oisson with a linear extrap ola- tion b elo w a b order v alue l , and add a quadratic function scaled to still k eep the minimum at 0 for the zero-photon case. Hence, negative v alues are allow ed, but strongly p enalized. In our implemen tation, we ha ve c hosen to implemen t the windowing factor inside f , rather than in A . This makes the effects of the windowing on the re- sulting gradients more visible (the gradients are scaled by the windo w squared). 9 Since the reason to in tro duce l is to ensure a smo oth behavior amenable to opti- mization, it makes sense that the effectiv e gradien t should b e of equal magnitude ev erywhere. Therefore, the actual v alue used for the barrier width is in verted b y the window, resulting in a more lenien t barrier in the high-frequency areas, but also requiring an o verall lo w er v alue of l to guaran tee high accuracy at all angles. T o stabilize the in verse pro cess, the Hann window in amplitude space (whic h is squared in in tensity space) has an extra term of 10 − 3 in order to av oid singularities. 4.3 Relaxed quadratic barrier and accelerated contin ua- tion The extra quadratic barrier l in tro duced in the previous section determines the sharpness of the non-negativity constraint. In our implementation, we initiate l to a high v alue ( l = 4). This v alue is then decreased by a factor of . 5 rep eatedly un til reaching 2 − 46 ≈ 10 − 14 . Due to the Hann windo w scaling of l , this means an effective barrier of 10 − 8 at the edge of the image. This sc heme ensures a v ery quic k conv ergence to the relaxed problem, fol- lo wed b y successiv e sharp enings in a contin uation scheme, starting off from the previous end result. Let X i denote the result from iteration i (starting at 0, with l = 2 2 − i ). Since the ma jor trend influencing the results should b e the reduced fluctuations around 0, the initial guess for iteration i , when i > 2, will b e X i − 1 + 0 . 5( X i − 1 − X i − 2 ), i.e. repro ducing half the c hange that resulted from the previous halving of l . This acceleration scheme is similar to the ones more thoroughly inv estigated for the SCD mo del in [13]. Naturally , the actual X v alues used are translated for the accuracy reasons already discussed. Each con tinuation step with a new l v alue induces a rapid shift in the optim um, even after the acceleration adaption. Therefore, a second set of inner contin uations is applied after a fixed n umber of individual inner iterations of the optimization scheme. In these inner contin uations, a separate acceleration sc heme is used after the second iteration, where the new starting p oin t is displaced with an accumulation factor of . 9. If at an y point the resulting X estimate is found to ha ve a shorter distance to the old end p oin t than the new starting p oint, the acceleration is deemed unsuccessful and that inner iteration is restarted without ac celeration. W e ha ve found th is ad-hoc sc heme to give very reasonable results with an acceptable time usage. The inner con tinuation loop is terminated when the c hange in the ob jectiv e function, scaled b y the v alue of l , falls b elow a pre-defined threshold. 4.4 Windo w-a ware supp ort constraint The supp ort constraint can b e simply implemented as a pro jection to 0 outside of the supp ort. How ev er, this might unnecessarily constrain the mo del in terms of inducing very small steps. W e therefore v entured in to employing a strong, but non-infinite, quadratic penalty , chosen to be 5 e 7 l . 10 The quadratic p enalty , as well as the original zero-pro jection constrain t, has a drawbac k in terms of the sp ectral b eha vior. The sharp cutting at the edge of the supp ort will induce high-frequency signals in the diffraction space, while w e are already employing the Hann window to the original diffraction signal in order to reduce that signal. Therefore, w e are using the same kind of tw o- lev el Hann window in tensit y space for the supp ort penalty matrix in the F ourier space, which amounts to a gradual decrease in support penalty ov er a width of 5 pixels. In practice, similar results were achiev able without this addition, but the time usage w as muc h higher. This can b e understoo d by the fact that in our sim ulations, the support was in fact wider than the ob ject. Therefore, the ideal solution is already zero in the b order region. F or a fully conv erged, ideal, solution, the tw o constrain t formulations are therefore identical. The window ed v ersion will only accelerate con v ergence b y a voiding sp ectral leak age. 5 Data and exp erimen ts Sim ulated diffraction data were pre-treated with our pro of-of-concept imple- men tation of COA CS based on TFOCS, called jackdaw (av ailable at https: //github.com/cnettel/jackdaw ). Sim ulations were done using the Condor [24] pack age. The test particle c hosen in the sim ulations consisted of an icosahedron, with a single v ery high- densit y sphere placed on a p oin t off-center on one of the vertices. The shortest diameter b et ween v ertices w as 20 pixels, imaged on a virtual 256x256 detector. The high-densit y sphere had a diameter of 4 pixels, with a muc h higher scat- tering pow er p er volume than the icosahedron. The presence of a concentrated, mark edly non-symmetric feature makes it p ossible to judge the presence of arti- facts in resulting reconstructions, as well as the ability to resolve the feature at all. At the same time, our simulated sample is as a first approximation roughly spherical, which is appropriate for the types of biological particles one w ould w ant to image. A total of 50 P oisson samplings were created of this particle. A central b eamstop w as sim ulated b y making a 25x25 square stencil in the cen ter missing. This co vered the ma jorit y of the cen tral sp ec kle for the icosahedron, resulting in on a v erage 10,000 photons in the non-mask ed pattern. Intensit y-healed patterns w ere computed using our COA CS implemen tation for this pattern, using a trivial auto-correlation supp ort consisting of a square with a side of 61 pixels. While a Shrinkwrap approac h akin to what is common in phase retriev al [25] would b e p ossible, this c hoice was made to demonstrate the lac k of sp ecific assumptions on the structure of the sample, and ho w it w ould be possible to pro cess a large n umber of shots of v arying particles with little to no man ual tuning. Phasing w as performed using the Python interface to libspimage with the GPU-based HIO phasing method, with the relaxation parameter β = 0 . 9, iden- tically for COA CSed (healed) and un-pro cessed images. The support used w as a square with a side of 31 pixels, i.e. slightly larger than the particle itself. 11 50,000 iterations of HIO w ere p erformed, follow ed by 10,000 iterations of pure ER to find a suitable lo cal minim um. F or eac h pattern, phasing w as repeated 100 times, in order to control for the v ariations due to random initial data. The ev entual reconstruction chosen for eac h pattern w as the a verage of the 10 b est out of these phasing results, chosen on the basis of the real space error, i.e. the norm of the remaining signal outside of the supp ort. 6 Results A qualitativ e comparison of the results p ossible using in tensit y healing are giv en in Figure 2, for the 50 patterns simulated. F or each pair, the CO ACsed pattern is show ed to the right. The high-densit y spherical feature on the edge is clearly visible in all healed reconstructions, and the contours of the icosahedron pro- jection are reasonably traced, while the traditional phasing metho d is only able to pro duce an elongated shap e without any detail, and a muc h more gradual fading into the outer areas of the supp ort. The diffraction pattern with v arious stages of pro cessing is shown for the first random instance in Figure 3. The reconstructed patterns lo ok similar to the original in Figure 1. The phased COA CS patterns and e specially the non- phased CO ACS patterns are clearly sup erior to the patterns achiev ed using traditional phasing. There is no p erceptible ringing effects from the highly regular rectangular supp ort used in the CO ACS patterns. The crystallographic R factor, or equiv alently the L1 norm of the difference b et w een the recov ered wave amplitudes and the true amplitudes, normalized b y the norm of the true amplitudes [26], is presented in T able 1. The R factor with no healing is more than 100% higher than what is obtained when phasing with a healing pre-pro cessing step. How ev er, even this result has an error significan tly larger than what is obtained from the healing step itself. Th us, the phasing pro cess is still a significant con tributor to errors. Since these calculations go all the wa y to the edge of the detector, where the sampled pattern was exceedingly sparse, the v alues are relativ ely high. T able 1: R factor means and standard deviations. Dataset Avg R-factor σ Av erage phasing without CO ACS healing 0.416 0.017 Av erage phasing with CO ACS healing 0.158 0.019 CO ACS healed 0.098 0.005 In Figure 4, the R factors calculated at v arious radii (in pixels) are also sho wn. The traditional phasing metho d is not able to correctly recov er the in tensity in the central missing data region (b eing a square of side 25, this is seen most up to a radius of 12.5). In the o verall minimum for the true signal due to the shape of the high-densit y spherical feature in the pattern at 90 pixels, the direct application of CO ACS gives inferior results due to the more sp ec klish 12 0.00 0.05 A.U. Figure 2: 50 phased reconstructions of sparse patterns based on the simulated particle. P airs of phased healed pattern using CO A CS (left) and phased original patterns (right). Av erage of 10 best individual phasings out of 100 replicates for each of these. 13 a) b) c) d) e) f) 0.0 1.5 Photons Figure 3: Diffraction pattern results for the first random instance (rescaled to undo Hann windo wing). a) Original sim ulated pattern. b) P oisson sampling with simulated 25x25 b eamstop. c) Average diffraction pattern ov er 10 b est phasings. Artifacts at the edges, repro duction of original features dropping off quic kly . d) Average diffraction pattern o v er 10 b est phasings based on COA CS healing. F eatures repro duced to higher angles, but s till problems at the edge, missing the outer ring. e) Diffraction pattern resulting directly from the CO ACS healing. Artifacts present in d) are absent. Outer ring corresp onding to the one found in original pattern. f ) Average based on the separate healing results of all 50 sampled particles. Errors cancel out, res ulting in an image v ery similar to a). 14 nature of that pattern. Although the relativ e error is higher at this p oin t, the absolute error is small nonetheless, since the true signal is close to 0. The healing pro cess required approximately 30 minutes p er pattern on a 12-core serv er with Sandy Bridge I n tel Xeon cores. In the Discussion section, w e outline how this can b e impro v ed. 7 Discussion The time usage for the in tensity healing is muc h longer compared to phasing, 30 min utes vs. less than 10 seconds (for a single reconstruction within the batc h of 100). Nonetheless, b oth perform roughly the same num b er of iterations (on the order of 10 5 ). It should be noted that the curren t CO ACS implementation is a straigh tforward one in Matlab using the TFOCS toolb ox, while the phasing implemen tation libspimage [15] is efficien t GPU-based CUDA co de. When using the TF OCS mo de for coun ting the n umber of op erations, it is clear that the n umber of function and transform ev aluations scales linearly to the num b er of iterations. The adaptive backtrac king of the step-size in TFOCS do es mean that the n umber of ev aluations can b e higher, but in our testing this was still lo wer than a total factor of 5. A preliminary analysis of the computational w orkload indicates that the F ourier transforms should b e the dominating part of the computational load. Thus, an efficient implementation of the TF OCS algorithm should b e quite similar to the existing phasing, in terms of work p er iteration. The time difference is due to the difference b et ween a CPU- based implementation in Matlab for intensit y healing, and an efficient C++ implemen tation on GPUs using CUDA. Another av enue would be to use the existing p ort of TFOCS metho dology into the T ensorFlo w framew ork. W e can also conclude that our curren t acceleration sc heme is ad hoc rather than optimal. In fact, the translation of our problem in to one of general conv ex optimization means that one can muc h more easily identify and b enefit from existing acceleration schemes that hav e shown success in other applications. W e ha v e tried reformulations using the TFOCS SCD sc heme with accompanying relaxations of µ (not shown), but those did not p erform b etter than our ad ho c approac h in terms of ac curacy or performance. W e b eliev e that sligh t additional relaxations lik e the ones already outlined, and an efficient GPU-based implementation, should make in tensity healing achiev- able at time scales similar to phasing. More imp ortan tly , w e hav e demonstrated that a prop er n umerical treatmen t allo ws prop er intensit y healing and phasing without a carefully chosen tight supp ort or Shrinkwrap pro cedure. This dev el- opmen t holds the promise of allo wing routine phasing of recorded exp erimental diffraction patterns, without tedious man ual tuning. With equally tigh t support constrain ts, it is also foreseeable that this method will allow b etter handling of patterns with missing data in some parts of the pattern due to e.g. satura- tion, the presence of a beamstop to protect the detector from the muc h stronger non-diffracted beam, or other asp ects of detector or experiment geometry . F or analyzing real-w orld data, one will probably w ant to augment our curren t Pois- 15 0 20 40 60 80 100 120 Radius [pixels] 0.0 0.2 0.4 0.6 0.8 1.0 R factor Avg phasing w-o. COACS healing Avg phasing w. COACS healing COACS healed Figure 4: R factors (normalized relativ e L1 error) for v arious radius shells in pixels. Curves are av erages ov er the individually computed results for all 50 sim ulated particles. Comparison b et ween results based on av erage phasing of 10 b est reconstructions of original pattern, av erage phasing of 10 b est recon- structions of COA CS-healed pattern, and using the structural factors from the CO ACS pattern directly . COA CS-healing reduces phasing errors, but the phas- ing step is still a significan t contributor to errors. Peak at around 90 pixels is due to the R factor being a relative error measure. This is the lo cation of a minim um due to the shape of the small spherical feature. Hence, absolute errors of the same magnitude are amplified. 16 son probability distribution with a Gaussian comp onent for v alues close to 0, to b etter account for electronic noise in imaging detectors. Another observ ation based on our ev aluation is that the phasing sc hemes are still the w eakest step in the reconstruction procedure. Even the smo oth CO ACS patterns result in reconstruction processes that start to “walk” along the image. This is due to a non-ideal phase ramp b eing induced by the combination of the t wo constrain ts, resulting in mov ement in real space. This, in turn, means that the real space ob ject will repeatedly “bump” in to the supp ort, with new artifacts b eing introduced as the real space constrain t implementation tries to remo ve the signal. While the scop e of this publication is not to improv e phasing p er se, new insigh ts can be gained by clearly separating effects due to sparse data, the difference b et ween the con tinuous and discrete F ourier transforms, and the phasing metho dology itself. W e also note that a few individual reconstructions lo ok far sharp er than the av erages sho wn, although the error metrics do not reliably separate those. It is our opinion that a prop er phasing method should b e able to pro duce R factors similar to those we rep ort for the CO ACS metho d alone in T able 1 and Figure 4. The straightforw ard structure and generality of the con vex optimization for- malism as well as the TFOCS library in particular also mak e it feasible to add additional constraints. Suc h constraints migh t include a total v ariation norm in the auto correlation or F ourier space to regularize the problem. Espe- cially in exceedingly sparse cases, and with c hallenging experiment geometries, suc h additional inspiration from the compressed sensing literature migh t pro ve w orthwhile. 8 Conclusion W e hav e presen ted the COA CS approach to correct the sampled diffraction pat- tern based on a supp ort constrain t. This approac h allo ws for higher-resolution reconstructions, whic h made it p ossible for us to phase simulated data with a wide support and no specific tuning of the parameters with a high lev el of detail. W e hav e also identified several possible future developmen ts, most pressingly to implemen t a GPU-based version of our approach in order to present more com- p etitiv e p erformance in terms of computation time. References [1] P . D. Adams, P . V. Afonine, R. W. Grosse-Kunstlev e, R. J. Read, J. S. Ric hardson, D. C. Ric hardson, and T. C. T erwilliger, “Recen t dev elopmen ts in phasing and structure refinement for macromolecular crystallography ,” Curren t Opinion in Structural Biology 19 , 566–572 (2009). [2] J. R. Fienup, “Phase retriev al algorithms: a comparison,” Applied Optics 21 , 2758–2769 (1982). 17 [3] D. R. Luk e, “Relaxed a veraged alternating reflections for diffraction imag- ing,” Inv erse problems 21 , 37 (2004). [4] S. Marc hesini, “Invited article: A unified ev aluation of iterative pro jection algorithms for phase retriev al,” Review of Scien tific Instrumen ts 78 , 011301 (2007). [5] S. Boutet and G. J. Williams, “The coherent X-ray imaging (CXI) instru- men t at the Linac Coherent Ligh t Source (LCLS),” New Journal of Ph ysics 12 , 035024 (2010). [6] A. V. Martin, J. Andreasson, A. Aquila, S. Ba jt, T. R. Barends, M. Barthelmess, A. Barty , W. H. Benner, C. Bostedt, J. D. Bozek et al. , “Single particle imaging with soft x-rays at the linac coheren t light source,” in “Adv ances in X-ra y F ree-Electron Lasers: Radiation Sc hemes, X-ra y Op- tics, and Instrumentation,” , v ol. 8078 (In ternational Society for Optics and Photonics, 2011), vol. 8078, p. 807809. [7] H. N. Chapman, A. Barty , M. J. Bogan, S. Boutet, M. F rank, S. P . Hau- Riege, S. Marc hesini, B. W. W o o ds, S. Ba jt, W. H. Benner et al. , “F em- tosecond diffractiv e imaging with a soft-X-ray free-electron laser,” Nature Ph ysics 2 , 839 (2006). [8] M. M. Seibert, T. Ekeberg, F. R. Maia, M. Svenda, J. Andreasson, O. J¨ onsson, D. Odi´ c, B. Iw an, A. Ro c ker, D. W estphal et al. , “Single mimivirus particles in tercepted and imaged with an X-ray laser,” Nature 470 , 78 (2011). [9] M. F. Hantk e, D. Hasse, F. R. Maia, T. Ekeberg, K. John, M. Svenda, N. D. Loh, A. V. Martin, N. Timnean u, D. S. Larsson et al. , “High-throughput imaging of heterogeneous cell organelles with an x-ra y laser,” Nature Pho- tonics 8 , 943 (2014). [10] E. J. Candes, T. Strohmer, and V. V oroninski, “Phaselift: Exact and stable signal reco very from magnitude measurements via conv ex programming,” Comm unications on Pure and Applied Mathematics 66 , 1241–1274 (2013). [11] E. J. Candes, Y. C. Eldar, T. Strohmer, and V. V oroninski, “Phase retriev al via matrix completion,” SIAM Review 57 , 225–251 (2015). [12] I. W aldspurger, A. d’Aspremont, and S. Mallat, “Phase recov ery , maxcut and complex semidefinite programming,” Mathematical Programming 149 , 47–81 (2015). [13] S. R. Bec ker, E. J. Cand` es, and M. C. Gran t, “T emplates for con v ex cone problems with applications to sparse signal recov ery ,” Mathematical Pro- gramming Computation 3 , 165 (2011). 18 [14] T. Ekeberg, M. Svenda, C. Ab ergel, F. R. Maia, V. Seltzer, J.-M. Cla verie, M. Hantk e, O. J¨ onsson, C. Nettelblad, G. V an Der Schot et al. , “Three- dimensional reconstruction of the gian t mimivirus particle with an x-ray free-electron lase r,” Physical Review Letters 114 , 098102 (2015). [15] F. R. Maia, T. Ek eb erg, D. V an Der Spo el, and J. Ha jdu, “Ha wk: the image reconstruction pack age for coherent x-ray diffractiv e imaging,” Journal of Applied Crystallography 43 , 1535–1539 (2010). [16] P . Thibault, V. Elser, C. Jacobsen, D. Shapiro, and D. Sa yre, “Reconstruc- tion of a y east cell from x-ray diffraction data,” Acta Crystallographica A 62 , 248–261 (2006). [17] L. Shi, G. W etzstein, and T. J. Lane, “A flexible phase retriev al framework for flux-limited coherent x-ray imaging,” arXiv preprint (2016). [18] Y. Nestero v, Intr o ductory le ctur es on c onvex optimization: A b asic c ourse , v ol. 87 (Springer Science & Business Media, 2013). [19] J. L¨ ofberg, “Y almip: A to olbox for mo deling and optimization in matlab,” in “Computer Aided Control Systems Design, 2004 IEEE In ternational Symp osium on,” (IEEE, 2004), pp. 284–289. [20] M. ApS, The MOSEK optimization to olb ox for MA TLAB manual. V ersion 8.1. (2017). [21] M. Grant, S. Bo yd, and Y. Y e, “Cvx: Matlab soft ware for disciplined con vex programming,” (2008). [22] Absorption and Sc attering by a Spher e (Wiley-Blac kwell, 2007), chap. 4, pp. 82–129. [23] S. Beck er, E. Candes, and M. Grant, T emplates for first-or der c onic solvers user guide 1.3 (2013). [24] M. F. Hantk e, T. Ekeberg, and F. R. Maia, “Condor: a simulation to ol for flash x-ray imaging,” Journal of Applied Crystallography 49 , 1356–1362 (2016). [25] S. Marchesini, H. He, H. N. Chapman, S. P . Hau-Riege, A. Noy , M. R. Ho wells, U. W eierstall, and J. C . Spence, “X-ra y image reconstruction from a diffraction pattern alone,” . [26] W. C. Hamilton, “Significance tests on the crystallographic r factor,” Acta Crystallographica 18 , 502–510 (1965). 19