Benchmark problems for phase retrieval
aa r X i v : . [ c s . I T ] J un BENCHMARK PROBLEMS FOR PHASE RETRIEVAL
VEIT ELSER ∗∗ , TI-YEN LAN ∗ , AND
TAMIR BENDORY †† Abstract.
In recent years, the mathematical and algorithmic aspects of the phase retrievalproblem have received considerable attention. Many papers in this area mention crystallography asa principal application. In crystallography, the signal to be recovered is periodic and comprised ofatomic distributions arranged homogeneously in the unit cell of the crystal. The crystallographicproblem is both the leading application and one of the hardest forms of phase retrieval. We haveconstructed a graded set of benchmark problems for evaluating algorithms that perform this type ofphase retrieval. The data, publicly available online , is provided in an easily interpretable format.We also propose a simple and unambiguous success/failure criterion based on the actual needs incrystallography. Baseline runtimes were obtained with an iterative algorithm that is similar butmore transparent than those used in crystallography. Empirically,the runtimes grow exponentiallywith respect to a new hardness parameter: the sparsity of the signal autocorrelation. We also reviewthe algorithms used by the leading software packages. This set of benchmark problems, we hope,will encourage the development of new algorithms for the phase retrieval problem in general, andcrystallography in particular. Key words. phase retrieval, crystallography, periodic signals, reconstruction algorithms, bench-mark problems, sparsity
AMS subject classifications.
1. Introduction.
The publication of “Phase retrieval via matrix completion”by Cand`es et al. [10] launched a revival of interest in the phase retrieval problem.Phase retrieval originated in X-ray crystallography, which is still by far the largestapplication. In 2016, about 50,000 crystal structures were deposited in the CambridgeStructural Database [31], each made possible by a phase retrieval algorithm. In brief,phase retrieval seeks to reconstruct a signal from measurements of its Fourier mag-nitudes. This is well-posed, of course, only if additional information is brought tobear on the reconstruction. In the case of X-ray crystallography, where the signal hassupport only at the positions of atoms, the additional information takes the form ofa sparsity constraint.All algorithms currently in use for crystallographic phase retrieval are close de-scendants of algorithms developed in the 1970s, and are viewed as heuristic methodsby today’s standards. By contrast — see [3, 46] for recent surveys — the current waveof phase retrieval research has produced algorithms that can guarantee solutions. Isa migration to such algorithms imminent in crystallography, or does crystallinity —signal periodicity — pose new, fundamental challenges? In this paper, we make acase for the second scenario, and make contributions designed to help resolve whatwe see as an under-appreciated theoretical problem.In Section 2, we describe the crystallographic phase retrieval problem, highlightingthose features that bring it outside the scope of current theoretical research andsignificantly increase its hardness. Our main contribution is presented in Section 3,where we describe a set of benchmark problems. These are synthetic instances ofphase retrieval, designed with sufficient realism that success at solving them wouldtranslate to success with real data sets. Section 4 gives details on the benchmark ∗ Department of Physics, Cornell University, Ithaca, NY, 14853-2501 USA ([email protected],[email protected]). † The Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ,08544-1000 USA ([email protected]). https://github.com/veitelser/phase-retrieval-benchmarks1 VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY constructions and, in particular, defines a new index, peculiar to crystallography, forranking hardness. A representative heuristic method, described in Section 5, exhibitsexponential time behavior with respect to this index. An easy way for a new brandof algorithm to distinguish itself is to demonstrate, at least empirically, behavior withsmaller exponential growth or even sub-exponential behavior.Heuristic methods in phase retrieval are often disparaged for lacking polynomial-time guarantees. This characterization, in the case of crystallographic phase retrieval,should be revisited if the new theoretical ideas fail to produce polynomial-time algo-rithms. As explained in Section 5, the (empirical) exponential behavior of the heuristicmethods is still far superior to na¨ıve exhaustive search. In addition, these methodsappear to produce solutions reliably and their performance on the benchmarks inSection 5 is strong motivation for their continued study.Much of the current theoretical work on phase retrieval, for instance [13, 20, 54],was inspired by the proposal of Cand`es et al. [10] of collecting a greater volume ofmagnitude-only data through the use of sensing vectors (designed or random) imple-mented as phase masks in the X-ray experiment. This proposal overlooks a physicallimit posed by the extreme smallness of the electron/X-ray interaction. For example,in a synchrotron X-ray diffraction experiment to determine the structure of a 30-nmvirus particle, X-rays are elastically scattered into the detector (the source of all infor-mation) only at a rate of 10 − X-rays per particle over the duration of the experiment.Unless data is collected simultaneously from a very large number of identical particles,the number of bits of data arriving at the detector falls far short of what is required todefine the particle’s structure (typically, the positions of many thousands of atoms).Therefore, data is always collected from crystals because it is the only feasible methodof aligning the required large number of particles to enable the simultaneous record-ing from all of them. While periodicity of the particles within the crystal is clearlynecessary for an adequate signal-to-noise ratio, the exact same periodicity is requiredof the proposed phase mask for obtaining additional magnitude-only data. In thissetting, a “phase mask” would be a designed particle that co-crystallizes with theunknown target particle. This is a much used strategy for structure determinationcalled “molecular replacement,” where the designed particle is a previously solvedstructure known to be a constituent part of the unknown particle. However, in molec-ular replacement the known information combines additively with the signal vector,not multiplicatively.Because of physical limits that apply when resolving sub-nanometer detail, thereis much less opportunity for designed measurement than at larger scales, such as ispossible in medical imaging. Signal periodicity is therefore a non-negotiable featureof phase retrieval in structural biology. Section 6 reviews four leading phase retrievalalgorithms for biomolecular crystals. We do not present benchmark results on thesealgorithms but use this section to give an historical account of the evolution of phaseretrieval algorithms, as exemplified by the method in Section 5 for which we do providedetailed results. Section 7 summarizes the paper.There are many non-crystallographic variants of the phase retrieval problem, aris-ing in ptychography [51, 42], laser pulse-shape characterization [52], etc., and wecannot hope to review them here. However, to satisfy the curiosity of the many re-searchers working on the simplified form of phase retrieval [10], we document in theappendix the behavior of our baseline algorithm (
RRR ) adapted for this aperiodic,no-prior-information case.
ENCHMARK PROBLEMS FOR PHASE RETRIEVAL
2. Crystallography.
A crystal is characterized by the periodic arrangement ofa repeating structural unit, also known as the unit cell . In X-ray crystallography, the“signal” is the electron density function of the crystal,(1) ρ c ( x ) = X y ∈ S ˜ ρ ( x − y ) , where ˜ ρ is a compactly supported motif and S is a finite set of translation vectors.The crystal parameters are defined by a lattice Λ ⊂ R D and the set S ⊂ Λ correspondsin practice to a very large and compact subset. If we fix a translation z ∈ Λ whileincreasing the size of S , the relation ρ c ( x + z ) ≈ ρ c ( x ) becomes a better approximation,failing only for x near the surface of the crystal.X-ray experiments measure the Fourier intensities | ˆ ρ c ( q ) | of the crystal, whereˆ ρ c ( q ) = Z R D d x ρ c ( x ) e − i π q · x (2) = X y ∈ S e − i π q · y (cid:18)Z R D d x ′ ˜ ρ ( x ′ ) e − i π q · x ′ (cid:19) (3) = ˆ s ( q ) ˆ ρ ( q ) . (4)As the set S grows (includes more of Λ), the support of the function ˆ s ( q ) is increasinglyconcentrated on the dual lattice Λ ∗ .The q -behavior of the measured Fourier intensities(5) | ˆ ρ c ( q ) | = | ˆ s ( q ) | | ˆ ρ ( q ) | is dominated by the function | ˆ s ( q ) | : the phenomenon of Bragg peaks . In this modelof a perfect crystal, the Bragg peak width is determined by the size of S . In practice,other factors (e.g., small misalignments of domains within the crystal) dominate theBragg peak width. Nevertheless, the counterpart of | ˆ s ( q ) | for imperfect crystals willstill have the periodicity of Λ ∗ and comprise isolated peaks having widths on a scaleso small that the function | ˆ ρ ( q ) | can be approximated as a constant over each peak.These features combine to give the integral of | ˆ ρ c ( q ) | , over just one Bragg peak, a verysimple interpretation. First, since | ˆ s ( q ) | is periodic, its integral over a Bragg peak isindependent of the specific peak and just contributes an overall scale factor. Second,the small width of the peaks ensures that the integrals are just sampling the function | ˆ ρ ( q ) | on the discrete set Λ ∗ . The X-ray measurements thus give the magnitudes ofthe Fourier series coefficients ˆ ρ ( q ), q ∈ Λ ∗ , that define a strictly periodic function on R D / Λ:(6) ρ ( x ) = 1vol(Λ) X q ∈ Λ ∗ ˆ ρ ( q ) e i π q · x , where vol(Λ) denotes the volume of the crystal unit cell, and the phases of ˆ ρ ( q ) stillhave to be retrieved to reconstruct the density ρ ( x ).The phase retrieval problem can also be understood through autocorrelation anal-ysis. In particular, by the convolution theorem, the signal autocorrelation(7) a ( y ) = a ( − y ) = Z d x ρ ( x ) ρ ( x + y ) , VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY is given directly by the inverse Fourier transform of the signal’s Fourier intensities | ˆ ρ ( q ) | . It is instructive to compare the autocorrelation for an aperiodic signal, where ρ and a are defined on R D , to the periodic case where these are functions on R D / Λ.Figure 1(a) shows a signal comprising N = 25 “atoms” that may be interpreted eitheras the repeating motif of a crystal or the aperiodic signal of a single “molecule.” Theautocorrelations a ( y ) in both cases comprise N ( N −
1) peaks associated with all pairsof atom-atom separations (and a strong peak at the origin for the N self-correlations).In the periodic case, Figure 1(b), the same number of autocorrelation peaks is crowdedinto a smaller region than the aperiodic case, Figure 1(c). The resolution of thesignal (atom size) and number of atoms in the illustration was set such that it is justbecoming difficult to resolve individual autocorrelation peaks in the periodic case. Bycontrast, the aperiodic autocorrelation still has well resolved peaks, particularly atthe periphery of the pattern, corresponding to the farthest separated atom pairs.The structure of the autocorrelation a ( y ), for a signal ρ ( x ) comprising atoms,is helpful for understanding how periodicity increases phase retrieval hardness, thatis, recovering ρ ( x ) from a ( y ). First, suppose that all the peaks in a ( y ) are so wellresolved that we know with high precision all the separation vectors, assumed to beunique. The atom positions in ρ ( x ) can then be inferred by hand, starting with a pairof atoms having one of the unique separations, then placing atoms one at a time, eachone constrained by having all separations to the placed atoms included among theunique separations in a ( y ). This is a polynomial-time algorithm and applies equallyin the periodic and aperiodic cases. But now consider the effects of finite resolution,that is, when the O ( N ) number of separations outpaces the number of resolutionelements, which normally scales as the number of atoms, O ( N ). Errors arising frommisidentified separations will occur with greater frequency in the periodic case becausethe autocorrelation peaks occupy a smaller region and have a higher density. Also,the density of peaks in the aperiodic case is not as uniform as it is in the periodiccase, to the degree that sets of widely separated atoms can be reconstructed quiteeasily from the periphery of a ( y ), where the peaks are sparse.The structure of the autocorrelation function, now for the periodic case, also in-forms us on how phase retrieval hardness depends on parameters. Suppose we matchthe size of resolution elements to the size of atoms — a reasonable approximation ofwhat is done in crystallography. If our unit cell has M resolution elements and N atoms, we would characterize the signal sparsity as N/M . However, as a comparisonof Figs. 1(a) and 1(b) makes clear, a signal sparsity of
N/M translates to an auto-correlation sparsity of order N /M . Therefore, we should expect the hardness of theinstances to scale like N /M rather than the standard measure of signal sparsity.An interesting consequence of the scaling difference between signal and autocor-relation sparsity is the possibility of instances that are both hard and have uniquesolutions. Many of the standard hard feasibility problems exhibit the phenomenonthat the onset of solution non-uniqueness, as a function of a parameter in randomlygenerated instances, coincides with instances tuned to the maximum hardness. Insome cases, e.g., logical satisfiability, the only hard instances that can be generatedrandomly are right at the transition point. In crystallographic phase retrieval thesituation is very different. We first note that the information deficit for a generic sig-nal, defined as the ratio of the number of independent measurements to the numberof independent signal values, is 1 /
2. This follows from the fact that for a generalreal-valued signal ρ ( x ) sampled on M resolution elements there are M/ | ˆ ρ ( q ) | . By na¨ıve constraint counting, a unique solutiontherefore requires at least M/ ENCHMARK PROBLEMS FOR PHASE RETRIEVAL apriori , to be zero. And, as we argued above, the retrieval of even such sparse signalspose a challenge when their autocorrelations are not sparse.Phase retrieval has three fundamental, though benign, forms of solution non-uniqueness, frequently referred to as trivial ambiguities. The (real-valued) signalcan be translated, inverted through the origin, or changed in sign, without changingthe Fourier magnitudes. These ambiguities are associated with physical symmetries(arbitrariness of the sign of the electron charge, etc.) and do not compromise in-terpretability. A less trivial form of non-uniqueness, in the case of signals with wellresolved atoms, arises when the geometry of the atom positions has the homometricproperty. Suppose ρ and ρ ′ are two solutions not related by one of the trivial ambigu-ities. Since they have equal Fourier magnitudes, they have matching autocorrelationfunctions and therefore share the same set of atom-atom separation vectors. In theaperiodic case, Rosenblatt and Seymour [43] have shown that this only happens when ρ is the convolution of (atomic) signals ρ ( x ) and ρ ( x ). When this is the case, asymmetry-unrelated signal ρ ′ can be constructed from the convolution of ρ ( x ) and ρ ( − x ) that has the same autocorrelation. However, signals that “factorize” in thissense are highly non-generic because there will be multiple symmetry-unrelated pairsof atoms that have exactly the same separation. In the generic case, when the peaksin the autocorrelation function all correspond to unique atom-atom separations, therewill not be any homometric non-uniqueness. This conclusion also applies to the pe-riodic case. A similar result was recently derived for the one-dimensional discretecase [2, Theorem 2.3].The relevance of periodicity for hardness can be argued on a more abstract levelwith the help of a discrete model, called bit retrieval [23]. In bit retrieval, a signal isrepresented by the integer coefficients of polynomials. The case where the non-zerocoefficients are all equal to 1 comes closest to crystallography (in one dimension),where the polynomial x k corresponds to an atom located at position k . In the ape-riodic case, the signal polynomials b ( x ) are Laurent polynomials, Z [ x, /x ], while inthe periodic case (crystals) the polynomials belong to the ring Z [ x ] / ( x P −
1) for someinteger P that defines the period. In bit retrieval we are given the autocorrelationpolynomial a ( x ) = b ( x ) b (1 /x ) and asked to find a factorization a ( x ) = b ′ ( x ) b ′ (1 /x ),ideally such that b ′ ( x ) also has all nonzero coefficients equal to 1. In the aperiodiccase, the factors can be found in polynomial time by the LLL algorithm [34], whereasno efficient factorization algorithms are known for the ring Z [ x ] / ( x P − G . An element g ∈ G acts on the density func-tion ρ by the composition of a rotation R g (an orthogonal matrix) and a translation T g : g ( ρ ( x )) = ρ ( R g · x + T g ) . Thus in addition to ρ ( x ) = ρ ( x + y ) , y ∈ Λ , the density of a crystal also satisfies ρ ( x ) = g ( ρ ( x )) , g ∈ G. VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY
Fig. 1: (a) A signal ρ comprising N = 25 “atoms.” (b) The autocorrelation functionof ρ in the periodic case (the enclosing square is the unit cell). (c) The autocorrelationin the aperiodic case.The set of rotation matrices R g identify G with a point group (transformations thatfix the origin), while the set of pairs ( R g , T g ), together with the group of latticetranslations Λ, specify the crystal’s space group . The order | G | of the point group isusually denoted Z .Space groups are usually identified prior to phase retrieval, often by the systematicvanishing of Fourier magnitudes | ˆ ρ ( q ) | for special q . If G is non-trivial, it has theeffect of reducing the number of independent density samples in the crystal unit cellby the factor Z . Our benchmarks all have the trivial space group and avoid thiscomplication.
3. Description of the datasets.
In this section we describe our syntheticdatasets and what it means to solve an instance. Details on the construction ofthe benchmark problems are given in the next section and can be skipped by readersjust wishing to solve the benchmarks.Data sets are identified by an integer N , the number of atoms in the signal,and a suffix characterizing the difficulty of the problem: E (easy), M (medium), H(hard). All data have the same format: an M × M table of integers (photon counts)representing the measurements of Fourier intensities | ˆ ρ | of a signal ρ sampled on aperiodic M × M grid. All instances have M = 128. This size was chosen to discouragemethods [10, 33] that represent the signal in terms of a dense, M × M matrix onwhich a rank-1 (or low rank) constraint is imposed. In 3-D protein crystallographythe corresponding dense matrix could not be processed or even stored. Regarding thedimensionality of the signal, we believe this has no effect on phase retrieval complexityin the periodic case; two dimensions was chosen only for ease of rendering the signal.Figure 2 shows a rendering and excerpt of the data file for the easiest instance, data100E . Though the data look random, there is a systematic decrease in the pho-ton counts with increasing spatial frequency. The table of photon counts containsseveral zero entries because, in experiments, intensities are normally measured out to ENCHMARK PROBLEMS FOR PHASE RETRIEVAL data100E . The(0 ,
0) photon count at the upper left corner is not measured and set to zero.frequencies where the Fourier transform is so small in magnitude that few photonsare detected. The (0 ,
0) intensity (photons scattered with zero momentum changeare indistinguishable from photons that did not scatter at all) is never measuredand appears as a 0 in the data file. All other intensities have been symmetrized, | ˆ ρ ( p, q ) | = | ˆ ρ ( − p, − q ) | , because the electron density signal of X-ray crystallographyis real-valued. The data files comprise just the 128 ×
64 half-table of symmetrizedphoton counts. All entries in the 65 th row and 65 th column (not shown in Figure 2)are zero.Solving an instance entails the following. Square roots of the data are taken anddefine the Fourier magnitudes | ˆ ρ ( p, q ) | . The phasing algorithm being demonstratedreconstructs ˆ ρ (0 , > φ ( p, q ) of the periodic signal(8) ρ ( x, y ) = 1 √ M M − X p =0 M − X q =0 e i πM ( px + qy ) | ˆ ρ ( p, q ) | e iφ ( p,q ) . The solution phases must satisfy φ ( p, q ) = − φ ( − p, − q ) in order that ρ ( x, y ) is real.Solutions are required to be consistent with a prior constraint on the support S ofthe signal. The cardinality of S is exactly 8 N , that is, on average each of the N atomsis supported on 8 pixels. In practice, we can take S to be the set of pixels on which ρ ( x, y ) has its 8 N largest values because the signal is not only real but non-negative.Consistency with the support constraint is established by checking a power inequality.From the data, as well as the reconstructed (0 ,
0) intensity, the total Fourier power isgiven by(9) I F = M − X p =0 M − X q =0 | ˆ ρ ( p, q ) | = M − X x =0 M − X y =0 ρ ( x, y ) . In a successful reconstruction, the power in the support(10) I S = X ( x,y ) ∈ S ρ ( x, y ) , matches the Fourier power. Because of Poisson noise (details in Section 4), the powerin the support falls short of the power in the data. However, all the benchmark VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY
Fig. 3: A phase retrieval solution (left) and the as-constructed ground truth (right) forbenchmark data300E . These images are rendered on a grid of the same size, 128 × I S I F > . . An instance is declared to be solved when this criterion is met. Understandingthis criterion, and what makes it difficult to achieve, is the crux of the phase retrievalproblem. Note that the denominator of (11) is mostly known. All that is unknownabout I F is the contribution of the constant term in the Fourier series, | ˆ ρ (0 , | .However, there is very little freedom in the value of this term when ρ is required to besparse. By contrast, the power in the support, I S , depends on all the unknown phases,via equations (10) and (8). Only very special combinations of phases, in combinationwith a special ˆ ρ (0 , ρ that is highly sparse and for which the powerin the small support (nearly) matches I F . Criterion (11) makes no reference to aground-truth solution, nor does it make any assumptions about solution uniqueness.Finally, we note that algorithms for solving the benchmarks are not required towork with signals sampled on M × M grids. Computations may be performed on fineror coarser grids, or without a grid sampling of the signal at all. However, at the endof the computation the algorithm is required to output the phase angles and ˆ ρ (0 ,
4. Construction details.
Figure 3 shows a successfully reconstructed signal forbenchmark instance data300E next to the signal that was constructed, by the methoddescribed in this section, to produce the data (a translation and reflection throughthe origin was applied to the latter to aid comparison). Of the 128 pixels, only8 ×
300 have significant signal (appear gray). These are bandlimited signals and lessnoisy than they appear. Smooth contour plots, formed by zero-padding the Fouriertransforms on a 512 ×
512 grid and back transforming, are shown in Figure 4. We see“atoms” of two types with centers having sub-grid precision.To construct data for an instance with N atoms, first a set of N atom-center pixelson a 512 ×
512 grid was sequentially sampled, uniformly but with the constraint thatthe distance between pixels is at least 12. After downsampling the signal by a factor of4 to avoid atoms artificially centered on grid points, this gives a minimum separationof 3 pixels between atom centers, corresponding to 3 ˚A in physical units (a typicalatomic diameter) when the pixel resolution of the data is 1 ˚A (typical of high quality
ENCHMARK PROBLEMS FOR PHASE RETRIEVAL N/ ×
512 to 128 × p, q ) and ( − p, − q ).The same intensity rescaling factor was used in all the benchmarks. This hasthe effect that individual atoms have the same characteristics (amplitude, width)across all the benchmarks. It also implies the total photon count in the data setsis proportional to N , the number of atoms. This normalization convention can bedefended on information theoretic grounds: to reconstruct the types and positions(in a fixed field) of N atoms, the quantity of information should be proportionalto N . Since each detected photon delivers the same quantity of information, thisproportionality is maintained when the total photon count is also proportional to N .Recall from section 2 that in crystallography the hardness of phase retrieval isexpected to depend on the autocorrelation sparsity rather than the signal sparsity.Since the signal is comprised of N atoms, the autocorrelation (away from the origin)is also “atomic” in character, but with N ( N −
1) peaks corresponding to all the inter-atomic vectors. These peaks are distributed homogeneously in the unit cell becausethe atoms themselves are distributed homogeneously. A small departure from theuniform distribution, see Figure 1(b), is explained by the constraint on the minimumatom-atom distance. For large N , a suitable definition of the autocorrelation sparsity,or the density of autocorrelation peaks, is therefore(12) µ = N V , where V is a measure of the number of pixels that can be resolved by the data.The measure V should correspond to the number of Fourier samples in the data,taking into account the decay in signal with increasing frequency. Because real dataand our synthetic data is well characterized by Gaussian decay of intensity with0 VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY frequency, we chose to define an “effective” number of Fourier samples by the formula(13) V = X q ∈ Λ ∗ e − b | q | , where the sum is over the lattice dual to the crystal lattice Λ and b is a parameter.This definition applies in any number of dimensions. In real 3-D crystals, the intensitydecay is reported as the value of the Wilson B -factor [28], and b = B/
6. For large V , the sum in (13) can be approximated by an integral and we obtain, in threedimensions,(14) V ≈ (cid:18) πB (cid:19) / vol(Λ) , where the last factor is the volume of the crystal unit cell.Numerically performing the sum (13) for the Gaussian filter of the benchmarks,we obtain the formula(15) µ = (cid:18) N . (cid:19) , for the benchmark instances. The numbers of atoms N of the instances was chosento sample µ by roughly equal intervals, from µ = 2 . µ = 39.The benchmark signals all have trivial space group ( P Z , both N and V in (12) should be divided by Z . This has the effect of replacing µ by µ/Z . To the best of our knowledge and withthis generalized definition, µ ≈
13 is the hardest reported case of phase retrieval withreal data for comparable structures (lacking heavy atoms; see Section 6, Table 2).When the atom positions are uniformly and independently sampled, the Fouriercoefficients (before filtering) have (asymptotically) a complex-normal distribution bythe central limit theorem. The corresponding intensities are then exponentially dis-tributed, a phenomenon known as Wilson statistics [28]. This is a reasonable statis-tical model for the benchmarks, since the minimum distance constraints are ratherweak for the densities of atoms considered. In real data it may happen that severalintensities are unusually large by this model, and their existence can be exploited byclever algorithms. Conversely, phase retrieval appears to be more challenging whenthe data lacks such outliers. The extreme case was recently studied for two-valued one-dimensional signals [23], where it is possible to construct signals whose intensities areexactly equal. Since the intensity-distribution characteristics are clearly important,we implemented the following modification in the construction of the atom positions.To quantify the outlier content of the intensity distribution, we used the normal-ized second-moment of the intensities(16) i = h| ˆ ρ | ih| ˆ ρ | i , where h·i denotes a uniform average over the measured, non-zero frequencies. Thisstatistical measure is increased when the high intensity tail of the distribution is en-hanced, and decreases when the distribution is made more uniform. Without anyintervention, when atom positions are uniformly sampled (rejecting positions that vi-olate the minimum distance constraint), we obtain i ≈
4. To make such instances
ENCHMARK PROBLEMS FOR PHASE RETRIEVAL Algorithm 1
Relaxed-Reflect-Reflect (
RRR ) algorithm input | ˆ ρ | , | S | , β Fourier magnitudes, support size,
RRR parameter ρ ← rand() random initial signal i ← repeat ( ρ , S ) ← P ( ρ ; | S | ) support-size projection ρ ← P (2 ρ − ρ ; | ˆ ρ | ) Fourier magnitude projection ρ ← ρ + β ( ρ − ρ ) increment by the projection discrepancy i ← i + 1 increment counter until pow( ρ , S ) > .
95 termination criterion (11) output ρ , i phased input magnitudes (solution), iteration counteasier, we select an atom at random and propose a new random position (still sat-isfying the distance constraint), accepting the proposal whenever the value of i isincreased. These increases are small, and many such moves had to be made to arriveat the value i = 4 . i isdecreased, continuing until i = 3 .
5. Relatively few atom-position re-samplings wereneeded to arrive at the value i = 4 .
5. Baseline results.
To the best of our knowledge, the only known algorithmsthat reliably solve the benchmark problems are heuristic in nature. A common featureof these algorithms is that they act iteratively on the signal. To set a baseline forthe benchmarks, we have selected a simple exemplar called Relaxed-Reflect-Reflect(
RRR ) [23]. This section describes the algorithm, addresses some common misconcep-tions about this type of algorithm, and suggests some standards for reporting results.Almost all crystallographic phase retrieval algorithms repeatedly use a “Fouriersynthesis” step, where a signal is constructed from the measured Fourier magnitudesand some estimate of the phases. The simplest such operation is the Fourier magnitudeprojection, ρ → ρ = P ( ρ ), where ρ inherits its Fourier phases from an arbitrarysignal ρ and combines these with the Fourier magnitudes of the data, when available.When magnitude data is not available, say at frequency q , the Fourier transform at q is simply copied.Most algorithms also repeatedly do “direct space refinement”, where prior infor-mation is imposed on the signal. One of the simplest operations of this kind, whenthe signal is known to be sparse, is the support-size projection ρ → ρ = P ( ρ ), where ρ is unchanged on the | S | highest valued pixels and set to zero on the rest — inthe process establishing the support S of the signal. In RRR , the pixels belonging tothe support S are updated in each iteration. We note that support-size projection issignificantly weaker in constraining the signal than the analogous operation appliedto a known support region. Going further, the case of a known support region S thatis sufficiently compact so it avoids aliasing (in the crystal unit cell) reverts to theeasier, aperiodic phase retrieval problem.Pseudocode for RRR is shown in Algorithm 1. In addition to the operations P and P already described, the algorithm calls on a function to initialize the signaland another function that computes the power ratio (11) for terminating iterations.Empirically, although the evolution of the signal is acutely sensitive to initial condi-tions (see below), the number of iterations required to find solutions, when averagedover different initial conditions, is not. Our implementation used a pseudo-random2 VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY number generator for initialization.
RRR has one parameter, β , that lies between 0and 2. Our baseline results are for β = 0 . RRR iterations,(17) ρ ρ + β ( P (2 P ( ρ ) − ρ ) − P ( ρ )) , we see that the parameter β is like a time-step. However, the “flow” in RRR doesnot correspond to a gradient, when P and P are interpreted as projectors to ageneral pair of constraint sets. With β = 1, the RRR algorithm [23] generalizes, forarbitrary pairs of constraints, Fienup’s “hybrid input-output” (
HIO ) algorithm [26] forphase retrieval with a known support region, and coincides with the Douglas-Rachfordsplitting scheme when applied to partial differential equations [18, 1].Many authors refer indiscriminately to
HIO — iteration (17) with β = 1 — andthe algorithm that simply alternates the projections:(18) ρ P ( P ( ρ )) . In the general convex setting this is “von Neumann’s alternating projections”, whilein microscopy and optics it is referred to as, respectively, “Gerchberg-Saxton” [27]and “Fienup’s error reduction” [26]. Aside from acting on the same space and beingbuilt from the same pair of projections, schemes (17) and (18) have almost nothingin common. Alternating projections is never used for hard (crystallographic) phaseretrieval. In just a few iterations, it converges to one of a multitude of uninterestingfixed points: signals with correct Fourier magnitudes that are proximal to, but notcoincident with, a signal of the correct support size. By contrast, when the
RRR algorithm has a fixed point, it is because the correctly supported signal ρ = P ( ρ ) isin the range of P — it also has the correct Fourier magnitudes. HIO , as originally formulated [26], also has a parameter β , but it does not cor-respond to a time-step and for β = 1, and non-linear P , loses the property thatfixed points are automatically associated with solutions. We note that the support-size choice for P in crystallography, unlike the known-support-region choice for HIO ,is highly non-linear. Marchesini [37] has proposed an optimization interpretation of
HIO that applies for general β . A predecessor of RRR was the “difference map”, adifferently parameterized composition of the two projections with the property thatreversing the sign of β interchanges them [21]. Constraint infeasibility, caused by noisein phase retrieval applications, was addressed by the “relaxed averaged alternatingreflection” ( RAAR ) [36] modification of
HIO .It is also not accurate to characterize
HIO / RRR as “alternating” or “cyclic” in theusual sense, as portrayed in popular accounts [46]. To see this, note that a fixed point ρ of (17) is in general not a solution. Instead it is the signals ρ and ρ generated bythe projections and equal at a fixed point that are solutions (see RRR pseudocode) .The truer sense in which this algorithm alternates is brought out by its equivalence, for β = 1, to the “alternating direction method of multipliers” ( ADMM ) algorithm [22].Contemporary accounts often are dismissive of projection-based algorithms be-cause convergence is not guaranteed, citing reports of “stagnation” in the behaviorof the iterates. While this criticism certainly applies to alternating projections, thedirect opposite is empirically the case for
RRR . The dynamics of
RRR has the char-acteristics of a strongly mixing system in mechanics, where ergodic behavior is therule rather than the exception. It is for this reason that initialization is unimportant. If ρ ∗ = ρ = ρ is a solution, then ρ may be any point in the set P − ( ρ ∗ ) ∩ (cid:16) ρ ∗ − P − ( ρ ∗ ) (cid:17) .ENCHMARK PROBLEMS FOR PHASE RETRIEVAL data100H with RRR . Not only is the final capture by the solution fixed point very brief, so is thetransient from the random initial signal to the family of signals explored in the search.Incorrectly phased signals in the long epoch of search have about 55% power in adynamic support constrained only by size.As in mechanics, the evidence of ergodicity is very strong even while prospects of aproof are dim. Ergodicity of the
RRR dynamics would apply in a strict sense onlyfor infeasible instances, when there are no fixed points. Because of noise, real-worldinstances are infeasible and the algorithm is remarkable in its ability to escape evennear-solutions. Near-solution infeasibility is not an issue for finding phase retrievalsolutions because we terminate the iterations as soon as the power ratio criterion (11)is satisfied. Every one of our runs of
RRR to solve a benchmark problem produced asolution: the success rate was 100%.There is no better illustration of the statistical behavior of the algorithm’s searchthan the time series of the power ratio (11). A typical time series for instance data100H is shown in Figure 5. The solution — marked by the sudden jump — is not constructedincrementally but appears as an isolated event, when the chaotic dynamics arrives bychance at a fixed point’s basin of attraction.
RRR is also used in convex optimization,where the behavior is very different and in fact convergent. However, when thesealgorithms are applied to hard phase retrieval only the final capture-phase of thesolution process — local convergence to the intersection of two affine sets — fallsunder the purview of convex analysis.The power ratio time series also illustrates the limits of using just the supportsize to constrain the signal. Figure 6 shows how the plot in Figure 5 changes as thesupport size increases. The hardest benchmark instances, with N = 400 atoms, arejust short of the point where criterion (11) fails as a valid certificate. It is only for thisreason that the benchmarks do not go beyond N = 400 ( µ = 39). Algorithms thatseek signals with additional prior characteristics — e.g., peaks — could in principlesucceed beyond this limit on the number of atoms. Indeed, algorithm developersare encouraged to exploit any of the prior signal information given in Section 4.Criterion (11) should only be seen as a convenient solution certificate that holds for N ≤ VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY
Fig. 6: Power ratio time series generated by
RRR for instances with progressivelymore atoms showing just the steady state behavior prior to the solution-discoveryjump seen in Figure 5. Beyond 400 atoms the power-in-support fraction is too closeto the value 0.95 (set by noise) to serve as a solution certificate.are best reported in terms that do not depend on implementation details. In thecase of
RRR , the average number of iterations in repeated trials serves this purpose .Success rates, when greater than zero, are more useful when converted to an expectedruntime per solution. Had RRR been run with a bound on the number of iterations,an expected runtime could have been computed from the total number of solutionsfound and the total number of iterations performed.Iteration counts per solution for
RRR , averaged over 20 trials per instance, aregiven in Table 1. Figure 7 shows the behavior with respect to the autocorrelationsparsity parameter µ defined by (15). At each N , almost without exception, the Einstance is easier than the M instance, which in turn is significantly easier than theH instance. The behavior with µ is consistent with a simple exponential. Linear fitsto the logarithms of the iteration counts give the following factors by which the meancount grows when µ is increased by 1: 1.56 (E), 1.72 (M), 1.93 (H).We expect our baseline results to be an easy target. The best outcome of phaseretrieval algorithm development would of course be a fundamentally different algo-rithm, promising subexponential cost in the parameter µ . On the other hand, evenincremental improvements in the case of exponential algorithms will bring substantialdividends. Although we have not performed the experiment, the extrapolation of ourresults indicate that RRR would require roughly three cpu-years to solve data330H ,compared with the single hour needed for data330E .
6. State of the art.
Probably the last time a phase retrieval milestone — onreal data — was hailed was the Shake-and-Bake (
SnB ) solution of triclinic lysozymein 1998 [17]. The
SnB algorithm was the product of a long history of developmentsthat drew inspiration from various disciplines, including signal processing, probabilitytheory and iterative methods for solving non-linear equations. In this section, webriefly review the principles behind
SnB [38] as well as those used by three otherleading crystallographic packages:
SHELXD [48],
SIR2004 [8] and
SUPERFLIP [40].Because the earliest phase retrieval algorithms were developed in the pre-FFT era, The runtime per iteration in our implementation, about 1 msec, is essentially constant acrossall instances.ENCHMARK PROBLEMS FOR PHASE RETRIEVAL
RRR mean iteration counts (log ) N E M H100 1.87 2.15 3.01140 2.37 3.00 3.93175 3.23 3.55 5.20200 3.42 4.57 5.48225 3.47 5.12 6.92245 4.33 5.77 7.03265 5.81 6.02 7.60285 6.06 5.98 7.62300 5.55 6.97 9.15315 6.46 7.29 –330 6.58 8.41 –345 7.83 – –360 6.86 – –375 8.00 – –Fig. 7: Exponential growth of the mean iteration count for
RRR as a function of µ forthe three difficulty grades of instances. Linear fits to the logarithms of the iterationcounts give the following factors by which the mean count grows when µ is increasedby 1: 1.56 (E), 1.72 (M), 1.93 (H).they imposed prior information on the signal not directly in real-space, but indirectlyduring Fourier synthesis. The simplest such strategy, known as David Sayre’s “tangentformula” [44], is based on the observation that in a signal ρ comprised of equal atom-like distributions (e.g., Gaussians), the Fourier phases of ρ and ρ are the same. Thisopens up the possibility that iterating ρ P ( ρ ) might by itself produce as fixedpoints a signal that ( i ) has been synthesized from the known Fourier magnitudesand ( ii ) corresponds to an atomic distribution — at least for crystals of sufficiently6 VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY identical atoms. It was possible to efficiently implement this map working just withthe Fourier coefficients by expressing the transform of the square as the convolutionof the transforms and then approximating the convolution by only the terms whereboth Fourier factors have a large magnitude.
SHELXD and
SIR2004 have the optionto apply the tangent formula operation in alternation with direct space refinementoperations.The tangent formula modification of the Fourier synthesis operation ( P ) is justone way the alternating scheme (18) is made viable again, by eliminating a hostof uninteresting fixed points. The identical-atom model of the signal on which themethod is based is also the premise behind another modification of P . This is the SnB objective function on the phases that is first minimized before phases are combinedwith magnitudes [38]. The simplest form of the objective function is based on theobservation that the distribution of the product ˆ ρ ( q )ˆ ρ ( q )ˆ ρ ( q ) is invariant, for q + q + q = 0, under translation of all the atoms . Therefore, it exhibits anon-trivial dependence on the sum of the corresponding phases, φ + φ + φ . Thedistribution of this “triplet” phase depends on the data via the known magnitude ofthe cubic product. The resulting conditional distribution for the triplet phase can becalculated explicitly for the case of equal atoms uniformly distributed in the crystalunit cell and serves as a model for triplet distributions in a typical crystal [28]. SnB tries to bring the cosines of the triplet phases in line with their expectation in therandom model by minimizing a sum-of-squares objective function.It is important that modified Fourier synthesis, by either the tangent formula orthe
SnB objective function, is combined with a robust direct space refinement opera-tion, such as the support-size projection P . This is because the phase interventions,or modifications of P , are based on approximate models and should only serve tobias the search for phases. When the bias built into P has sufficient strength, oneimproves the probability that the signal P ( P ( ρ )) is also fixed by P and is there-fore a solution. After all, the bias in both methods (tangent formula, SnB objectivefunction) is derived from a direct space model. Phase retrieval often succeeds simplyby alternating a modified P and a direct space refinement P of some type. In SnB , SHELXD and
SIR2004 the P operation can impose prior information on the signalbeyond just the size of its support. This may include knowledge of the minimumatom-atom distance, the presence of a known number of heavy atoms, or even theexpected histogram of the signal values.The phase retrieval algorithm in SUPERFLIP also alternates between Fourier syn-thesis and direct space refinement, but unlike the other packages, it avoids false fixedpoints through a modification of the direct space operation [39]. When written interms of the
RRR projections,
SUPERFLIP iterates a close approximation of the map(19) ρ P (2 P ( ρ ) − ρ ) . The algorithm’s name is derived from the argument of P , wherein the sign of thesignal is reversed wherever it is judged not to be in the support and unchangedotherwise. In a solution, the “charge flipping” step has no effect and the signal isa fixed point of the map. On the other hand, one cannot theoretically rule out thepossibility of exotic non-solution fixed points, where charge flipping only changes theFourier magnitudes — that are then restored and the charge re-flipped by P . The P used by SUPERFLIP [39] differs from the one in
RRR by being parameterized througha small positive lower bound on the signal value in the support rather than a bound This is also the underlying idea of analysis using the bispectrum [53, 4].ENCHMARK PROBLEMS FOR PHASE RETRIEVAL
RRR is not based on an alternation of two operations,it is intriguing and perhaps not a complete coincidence that Fourier synthesis in thisalgorithm is also preceded by charge flipping.Direct comparison of the accomplishments of the crystallographic packages withthe benchmark problems is complicated by a number of factors. First and foremost isthe fact that data sets of sufficient quality, on which sparsity can be imposed, becomeincreasingly rare as the number of atoms in the unit cell grows. In crystals of largeprotein molecules there is too much variability in the electron density from one unitcell to another (solvent disorder) that individual protein atoms cannot be resolved.The signal encoded by the Fourier magnitudes in this case is that of the average ofthe atomic distributions, and the corresponding broadening of the support has madeit too weak to be used as a constraint. Almost all large protein structures are solvedwith the help of additional data, derived from atom-specific inelastic scattering. Alarge share of the credit for structures with N above 1000 that did not rely on suchadditional data and yielded to “direct methods” goes to the crystal growers whomanaged to significantly reduce disorder in their crystals.Another factor that complicates comparisons is the presence of heavy atoms.Large proteins often contain a minority admixture of heavy atoms, and it is well knownthat this makes phase retrieval easier, even when this information is not used. Thispoint is best explained by the structure of the autocorrelation. The presence of heavyatoms in a crystal produces strong autocorrelation peaks at the separations betweenthe heavy atoms, from which the heavy atoms can first be located. Subsequently, thepositions of the light atoms can be determined one at a time by the O ( N ) number ofautocorrelation peaks with intermediate signals, which correspond to the separationsbetween a heavy atom and a light atom. The benchmark problems have equal numbersof atoms of two types and therefore correspond to harder instances in this respect.Table 2 lists what we judge to be the hardest instances of successful real-dataphase retrieval for crystals with ( i ) resolved atoms and ( ii ) the fewest number ofheavy atoms. The highest value of µ in this list is near the middle of the benchmarkproblems.
7. Summary.
Phase retrieval can be decisive in the success of crystal structurediscovery. Available algorithms for periodic signals are heuristic and their success andruntime behavior is poorly documented. It is not normal practice for crystallographersto test algorithms with synthetic data and a known ground truth; failures with realdata are attributed to insufficient resolution or corrupted measurements, and usuallygo unreported. And when data quality is good, Nature does not always cooperate tocreate instances with graded hardness for the study of algorithm behavior.Phase retrieval theory moved into a new era of systematic study when it wastaken up by applied mathematicians about seven years ago. However, the algorithmsgenerated by this development have had no impact on phase retrieval for crystals.Periodicity of the signal in crystallography is not a minor property that can be recov-ered as a special case of the phase retrieval models that have been studied. In fact,the hardness conferred to the phase retrieval problem by periodicity seems to havecompletely escaped the notice of mathematicians.Our benchmark problems were designed to address the circumstances just de-scribed. Crystallographers should evaluate their algorithms on standard benchmarksand applied mathematicians should turn their attention to phase retrieval problemsthat appear in real-world applications. The benchmark problems will serve both ofthese needs and, we hope, open a dialog between the two communities. V E I T E L S E R , T I - Y E N L ANAN D T A M I R B E N D O R Y structure PDB or space group, Z resolution (˚A) N/Z , µ/Z iterations software ref.CSD entry heavy atoms or timehen egg white lysozyme – P
1, 1 0.85 1001, 10S – 3,400
SnB [17]alpha-1 peptide 1BYZ P
1, 1 0.90 408, 1Cl 1.73 4,500
SnB [41]acutohaemolysin 1MC2 C
2, 4 0.85 975, 18S 4.45 –
SnB [35]scorpion toxin II 1AHO P , 4 0.96 500, 10S 4.46 413,000 SnB [49]actinomycin D 1A7Y P
1, 1 0.94 270 1.40 –
SHELXD [45]feglymycin 1W7Q P , 6 1.10 828 6.39 – SHELXD [7]hen egg white lysozyme 4LZT P
1, 1 0.95 1001, 10S 11.06 –
SHELXD [47]human cyclophilin G 2WFI P , 4 0.75 1486, 2Mg 15S 11.17 60 min SHELXD [50]hirustasin 1BX7 P
2, 8 1.20 366, 12S 4.29 546 min
SIR2004 [9]pheromone ER-1 2ERL C
2, 4 1.00 303, 8S 4.72 19 min
SIR2004 [9]Kunitz domain C5 2KNT P , 2 1.20 460, 1P 6S 6.85 22 min SIR2004 [9]bovine ribonuclease 1DY5 P , 2 0.87 1894, 31S 12.53 131 min SIR2004 [9]2C N O PAWVEO P
1, 1 0.80 164 0.58 100
SUPERFLIP [39]2C . N O . GOFMOD P
1, 1 0.80 188 0.66 250
SUPERFLIP [39]apamin – P , 2 0.95 385 4.36 20,000 SUPERFLIP [19]Table 2: Large, nearly equal-atom structures, solved by direct methods. PDB: Protein Data Bank. CSD: Cambridge Structural Database.
ENCHMARK PROBLEMS FOR PHASE RETRIEVAL
PhasePack [14],illustrates the kind of problems that can arise when a healthy dialog is absent. One ofthe algorithms available in
PhasePack is called “Fienup”, and a scientist working in thefield would assume this is Fienup’s
HIO algorithm, the most widely adopted methodand blueprint for the
RRR generalization to general, non-convex, two-constraint prob-lems. However, the “Fienup” algorithm of
PhasePack is just a minor variant ofthe alternating-projection (Gerchberg-Saxton) algorithm and only suitable for convexproblems. It acquired this name, by accident it seems, because Fienup’s article intro-ducing
HIO also featured the alternating-projection (“error-reduction”) algorithm asa poor alternative.Benchmarks are most useful when they uncover trends in behavior. We believethe autocorrelation sparsity parameter µ , introduced here, will serve this purpose.The benchmarks instances sample this hardness parameter over a range that includesthe frontier of real-world phase retrieval and the sampling is fine enough to be usefuleven when runtimes grow exponentially. Acknowledgments.
V.E. and T-Y.L. received support from DOE grant DE-SC0005827. T-Y.L. was also supported by the Taiwan Government Scholarship toStudy Abroad. T.B. would like to thank Amit Singer for his advice and support andYonina Eldar and Mahdi Soltanolkotabi for helpful discussions.
REFERENCES[1]
H. H. Bauschke, P. L. Combettes, and D. R. Luke , Phase retrieval, error reduction algo-rithm, and fienup variants: a view from convex optimization , JOSA A, 19 (2002), pp. 1334–1345.[2]
R. Beinert and G. Plonka , Ambiguities in one-dimensional discrete phase retrieval fromFourier magnitudes , Journal of Fourier Analysis and Applications, 21 (2015), pp. 1169–1198.[3]
T. Bendory, R. Beinert, and Y. C. Eldar , Fourier phase retrieval: Uniqueness and algo-rithms , in Compressed Sensing and its Applications, Springer, 2017, pp. 55–91.[4]
T. Bendory, N. Boumal, C. Ma, Z. Zhao, and A. Singer , Bispectrum inversion with ap-plication to multireference alignment , IEEE Transactions on Signal Processing, 66 (2017),pp. 1037–1050.[5]
T. Bendory, Y. C. Eldar, and N. Boumal , Non-convex phase retrieval from STFT mea-surements , IEEE Transactions on Information Theory, 64 (2018), pp. 467–484.[6]
N. Boumal, B. Mishra, P.-A. Absil, R. Sepulchre, et al. , Manopt, a matlab toolbox foroptimization on manifolds. , Journal of Machine Learning Research, 15 (2014), pp. 1455–1459.[7]
G. Bunk´oczi, L. V´ertesy, and G. M. Sheldrick , The antiviral antibiotic feglymycin: Firstdirect-methods solution of a 1000+ equal-atom structure , Angewandte Chemie Interna-tional Edition, 44 (2005), pp. 1340–1342.[8]
M. C. Burla, R. Caliandro, M. Camalli, B. Carrozzini, G. L. Cascarano, L. De Caro,C. Giacovazzo, G. Polidori, and R. Spagna , SIR2004: an improved tool for crystalstructure determination and refinement , Journal of Applied Crystallography, 38 (2005),pp. 381–388.[9]
M. C. Burla, R. Caliandro, B. Carrozzini, G. L. Cascarano, L. De Caro, C. Giacovazzo,G. Polidori, and D. Siliqi , The revenge of the Patterson methods. i. protein ab initiophasing , Journal of applied crystallography, 39 (2006), pp. 527–535.[10]
E. J. Candes, Y. C. Eldar, T. Strohmer, and V. Voroninski , Phase retrieval via matrixcompletion , SIAM J. Imaging Sci., 6 (2013), pp. 199–225.[11]
E. J. Cand`es and X. Li , Solving quadratic equations via phaselift when there are about asmany equations as unknowns , Foundations of Computational Mathematics, 14 (2014),pp. 1017–1026.[12]
E. J. Candes, X. Li, and M. Soltanolkotabi , Phase retrieval from coded diffraction patterns ,Applied and Computational Harmonic Analysis, 39 (2015), pp. 277–299.[13]
E. J. Candes, T. Strohmer, and V. Voroninski , Phaselift: Exact and stable signal recovery VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY from magnitude measurements via convex programming , Communications on Pure andApplied Mathematics, 66 (2013), pp. 1241–1274.[14]
R. Chandra, C. Studer, and T. Goldstein , Phasepack: A phase retrieval library , arXivpreprint arXiv:1711.10175, (2017).[15]
P. Chen and A. Fannjiang , Fourier phase retrieval with a single mask by douglas–rachfordalgorithms , Applied and Computational Harmonic Analysis, (2016).[16]
Y. Chen and E. J. Cand`es , Solving random quadratic systems of equations is nearly as easyas solving linear systems , Communications on Pure and Applied Mathematics, 5 (2017),pp. 822–883.[17]
A. M. Deacon, C. M. Weeks, R. Miller, and S. E. Ealick , The Shake-and-Bake structuredetermination of triclinic lysozyme , Proceedings of the National Academy of Sciences, 95(1998), pp. 9284–9289.[18]
J. Douglas and H. H. Rachford , On the numerical solution of heat conduction problemsin two and three space variables , Transactions of the American mathematical Society, 82(1956), pp. 421–439.[19]
C. Dumas and A. van der Lee , Macromolecular structure solution by charge flipping , ActaCrystallographica Section D: Biological Crystallography, 64 (2008), pp. 864–873.[20]
Y. C. Eldar and S. Mendelson , Phase retrieval: Stability and recovery guarantees , Appliedand Computational Harmonic Analysis, 36 (2014), pp. 473–494.[21]
V. Elser , Phase retrieval by iterated projections , JOSA A, 20 (2003), pp. 40–55.[22]
V. Elser , Matrix product constraints by projection methods , Journal of Global Optimization,68 (2017), pp. 329–355.[23]
V. Elser , The complexity of bit retrieval , IEEE Transactions on Information Theory, 64 (2018),pp. 412–428.[24]
V. Elser and R. Millane , Reconstruction of an object from its symmetry-averaged diffractionpattern , Acta Crystallographica Section A: Foundations of Crystallography, 64 (2008),pp. 273–279.[25]
A. Fannjiang and W. Liao , Phase retrieval with random phase illumination , JOSA A, 29(2012), pp. 1847–1859.[26]
J. R. Fienup , Phase retrieval algorithms: a comparison , Applied optics, 21 (1982), pp. 2758–2769.[27]
R. W. Gerchberg and W. O. Saxton , A practical algorithm for the determination of thephase from image and diffraction plane pictures , Optik, 35 (1972).[28]
C. Giacovazzo , Direct phasing in crystallography: fundamentals and applications , vol. 8,Courier Corporation, 1998.[29]
M. X. Goemans and D. P. Williamson , Improved approximation algorithms for maximum cutand satisfiability problems using semidefinite programming , Journal of the ACM (JACM),42 (1995), pp. 1115–1145.[30]
M. Grant, S. Boyd, and Y. Ye , CVX: Matlab software for disciplined convex programming ,2008.[31]
C. R. Groom and F. H. Allen , The cambridge structural database in retrospect and prospect ,Angewandte Chemie International Edition, 53 (2014), pp. 662–671.[32]
J. Hoffstein, J. Pipher, and J. H. Silverman , NTRU: A ring-based public key cryptosystem ,in International Algorithmic Number Theory Symposium, Springer, 1998, pp. 267–288.[33]
K. Huang, Y. C. Eldar, and N. D. Sidiropoulos , Phase retrieval from 1D Fourier measure-ments: Convexity, uniqueness, and algorithms , IEEE Transactions on Signal Processing,64 (2016), pp. 6105–6117.[34]
A. K. Lenstra, H. W. Lenstra, and L. Lov´asz , Factoring polynomials with rational coeffi-cients , Mathematische Annalen, 261 (1982), pp. 515–534.[35]
Q. Liu, Q. Huang, M. Teng, C. M. Weeks, C. Jelsch, R. Zhang, and L. Niu , The crystalstructure of a novel, inactive, lysine 49 PLA2 from agkistrodon acutus venom an ultrahighresolution, ab initio structure determination , Journal of Biological Chemistry, 278 (2003),pp. 41400–41408.[36]
D. R. Luke , Relaxed averaged alternating reflections for diffraction imaging , Inverse problems,21 (2004), p. 37.[37]
S. Marchesini , Phase retrieval and saddle-point optimization , JOSA A, 24 (2007), pp. 3289–3296.[38]
R. Miller, G. T. DeTitta, R. Jones, D. A. Langs, C. M. Weeks, and H. A. Hauptman , On the application of the minimal principle to solve unknown structures , Science, (1993),pp. 1430–1433.[39]
G. Oszl´anyi and A. S¨ut˝o , Ab initio structure solution by charge flipping , Acta Crystallo-graphica Section A: Foundations of Crystallography, 60 (2004), pp. 134–141.ENCHMARK PROBLEMS FOR PHASE RETRIEVAL [40] L. Palatinus and G. Chapuis , SUPERFLIP–a computer program for the solution of crystalstructures by charge flipping in arbitrary dimensions , Journal of Applied Crystallography,40 (2007), pp. 786–790.[41]
G. G. Priv´e, D. H. Anderson, L. Wesson, D. Cascio, and D. Eisenberg , Packed proteinbilayers in the 0.90 ˚a resolution structure of a designed alpha helical bundle , Proteinscience, 8 (1999), pp. 1400–1409.[42]
J. M. Rodenburg , Ptychography and related diffractive imaging methods , Advances in Imagingand Electron Physics, 150 (2008), pp. 87–184.[43]
J. Rosenblatt and P. D. Seymour , The structure of homometric sets , SIAM Journal onAlgebraic Discrete Methods, 3 (1982), pp. 343–350.[44]
D. Sayre , The squaring method: a new method for phase determination , Acta Crystallograph-ica, 5 (1952), pp. 60–65.[45]
M. Sch¨afer, G. M. Sheldrick, I. Bahner, and H. Lackner , Crystal structures of acti-nomycin D and actinomycin Z3 , Angewandte Chemie International Edition, 37 (1998),pp. 2381–2384.[46]
Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev , Phaseretrieval with application to optical imaging: a contemporary overview , IEEE signal pro-cessing magazine, 32 (2015), pp. 87–109.[47]
G. Sheldrick, C. Gilmore, H. Hauptman, C. Weeks, R. Miller, and I. Uson , Part 16.direct methods , (2012).[48]
G. M. Sheldrick , A short history of SHELX , Acta Crystallographica Section A: Foundationsof Crystallography, 64 (2008), pp. 112–122.[49]
G. Smith, R. Blessing, S. Ealick, J. Fontecilla-Camps, H. Hauptman, D. Housset,D. Langs, and R. Miller , Ab initio structure determination and refinement of a scorpionprotein toxin , Acta Crystallographica Section D: Biological Crystallography, 53 (1997),pp. 551–557.[50]
C. M. Stegmann, D. Seeliger, G. M. Sheldrick, B. L. de Groot, and M. C. Wahl , The thermodynamic influence of trapped water molecules on a protein–ligand interaction ,Angewandte Chemie International Edition, 48 (2009), pp. 5207–5210.[51]
P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, and F. Pfeiffer , High-resolutionscanning x-ray diffraction microscopy , Science, 321 (2008), pp. 379–382.[52]
R. Trebino , Frequency-resolved optical gating: the measurement of ultrashort laser pulses ,Springer Science & Business Media, 2012.[53]
J. W. Tukey , The spectral representation and transformation properties of the higher momentsof stationary time series , in The Collected Works of John W. Tukey, D. R. Brillinger, ed.,vol. 1, Wadsworth,, 1984, ch. 4, pp. 165–184.[54]
I. Waldspurger, A. dAspremont, and S. Mallat , Phase recovery, maxcut and complexsemidefinite programming , Mathematical Programming, 149 (2015), pp. 47–81.[55]
G. Wang, G. B. Giannakis, and Y. C. Eldar , Solving systems of random quadratic equationsvia truncated amplitude flow , IEEE Transactions on Information Theory, (2017).[56]
H. Zhang, Y. Zhou, Y. Liang, and Y. Chi , A nonconvex approach for phase retrieval: Re-shaped wirtinger flow and incremental algorithms , Journal of Machine Learning Research,(2017).
Appendix A. phase retrieval without prior information.
Beginning withthe publication of [10], “phase retrieval” has become identified in the applied mathe-matics community with the following mathematical problem:Given a sensing matrix A ∈ C m × n and magnitudes | y | ∈ R m ≥ , obtain a signal ρ ∈ C n such that | y | = | Aρ | . (20)Under what conditions on A can a signal that is unique, up to a global phase, berecovered, and how can this be done efficiently? This is conceptually simpler thanthe original, crystallographic, phase retrieval problem in that there are no prior con-straints on ρ , such as non-negativity or sparsity. The loss of information on measuringmagnitudes is entirely offset by having m sufficiently large relative to n for suitable A .Most of the current work on phase retrieval is directed at developing algorithms forsolving (20). In this appendix we compare the performance of three leading algorithms2 VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY with
RRR , which is easily adapted to solve (20) .The most direct real-world application of (20) we are aware of is reconstruct-ing the 2-D projected contrast of an isolated object from its diffraction intensity andknowledge of its 2-D support S . Suppose the measurements extend to spatial frequen-cies that resolve the contrast to pixels of area r . The signal vector ρ will then have n = | S | /r components. Consider the most common case, where these are real-valued.The diffracted intensity will then have centro-symmetry and, when sampled at res-olution r , will be completely described by m = | S − S | /r measurements [24]. It isoften possible to define S by an enclosing mask of known size and shape. The successof algorithms depends critically on the ratio m/n = | S − S | / | S | . For example, atriangular S corresponds to m/n = 3.The sensing matrix A for the application just described would be an m × n subma-trix of a 2-D Fourier matrix. The n columns would correspond to the pixel positionsin S , and the m columns to detector measurement pixels, preferably optimized toimprove the condition number of A . We are not aware of any demonstrations of (20)along these lines, for real or simulated data, by the algorithms recently proposedby the applied math community. This may simply be a reluctance to try the newmethods on sensing matrices with “structure.” The comparisons in this appendix willtherefore be with popular random models for A .In our first comparison, both the entries of A and the ground truth ρ , in eachinstance, were drawn from the normal distribution with zero mean and variance 1 / A and the magnitudes | y | = | Aρ | . Since the signal ρ can only be estimated up to a global phase, we definethe normalized estimation error between ρ and the ground truth ρ by(21) error = min φ ∈ [0 , π ) k ρ − e iφ ρ k / k ρ k . In the two-constraint formulation on which the
RRR algorithm is based, we seekvectors y ∈ C m which ( i ) are in the range of A , and ( ii ) have the known magnitudes | y | . The corresponding constraint projections are P ( y ) = AA † y, (22) P ( y ) = | y | e i arg y , (23)where A † = ( A ∗ A ) − A ∗ is the pseudo-inverse of A . The RRR parameter β was setat 0.5 and ρ was initialized naively, each element drawn from the complex normaldistribution. Iterations were terminated when the normalized norm of the differenceof successive iterates, k y ′ − y k / k y k , dropped below 10 − . The resulting y is ap-proximately fixed by both projections and gives the solution estimate ρ = A † y . Oursolver was written in Matlab and did not require any special packages.We compared
RRR with three recently proposed algorithms designed to solveinstances of (20). The first method, called
PhaseLift , is based on the insight that (20)can be formulated as a set of linear equations with respect to the rank-1 Hermitianmatrix ρρ ∗ . Applying a standard convex relaxation to the rank-1 constraint leads toa semidefinite program (SDP) that can be solved in polynomial time. The signal isestimated as the leading eigenvector of the SDP’s solution. This technique has beenthoroughly analyzed in a series of papers; see for instance [10, 13, 11, 12, 54]. Weused CVX toolbox [30] to solve the SDP. A software interface for a wide range of phase retrieval algorithms is provided in [14].ENCHMARK PROBLEMS FOR PHASE RETRIEVAL ❳❳❳❳❳❳❳❳❳❳ algorithm m/n RRR
NCPC
TWF
PhaseLift — 186.4 19.4Table 3: Average CPU time in seconds per successful reconstruction for different m/n values.The second algorithm, Truncated Wirtinger Flow (
TWF ), was proposed in [16] andminimizes a non-convex least-squares objective by a gradient algorithm. The initialsignal estimate is generated by a spectral algorithm. We used the implementationprovided by the authors with 100 iterations for the initialization and 2000 gradientiterations, with step size 0 .
2, for minimization. Several papers have proposed similartechniques based on different objective functions; see for instance [55, 56]. However,we found these modifications had little effect on performance.The last algorithm, proposed in [5], is called Non-Convex PhaseCut (
NCPC ) and isbased on the PhaseCut algorithm [54]. Whereas the original PhaseCut algorithm usedthe classical Max-Cut SDP relaxation [29], the
NCPC algorithm minimizes insteada non-convex objective on the manifold of phases. In particular,
NCPC employsa Riemannian trust-regions algorithm using the
Manopt toolbox [6]. The signal isinitialized by the spectral algorithm of [16] with 100 iterations. We interrupted thetrust-regions iterations when the norm of the gradient dropped below 10 − .Figure 8 plots the success rates of the four algorithms when recovering complexsignals of length n = 50 as a function of the ratio m/n . For each m/n , 100 trials wereperformed and a trial was declared successful if the error (21) was below 10 − . Thesame ρ and A were given to all the algorithms in each trial. Recall that there are 2 n parameters in a successful reconstruction, the real and imaginary parts of the signal.As can be seen, RRR is the only algorithm that achieves a perfect success rate closeto this information limit.Table 3 gives the total time used by each algorithm over all 100 trials, dividedby the number of successful reconstructions. These values should be interpreted asthe average time needed to achieve a successful reconstruction. We see that
RRR outperforms the other algorithms over all the m/n values considered. These resultsshould be taken with caution since the algorithms used different stopping criteria.For instance, the
TWF code we adopted from the authors’ website halts after a givennumber of iterations (in our case, 2000). By contrast, the stopping criteria of theother algorithms are based on the progress of the algorithm (e.g., when the norm ofthe gradient is smaller than some predefined value.) The parameters for the stoppingcriteria (reported above) were not optimized to minimize the time per successfulrecovery.Cand`es et al. [10] also studied instances of problem (20) where the sensing matrix A corresponds to the following case of the coded diffraction pattern (CDP) model.Given L random phase masks whose entries d ℓ ( s, t ) are drawn independently and http://statweb.stanford.edu/ ∼ candes/TWF/code.html VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY s u cc e ss r a t e RRRNCPCTWFPhaseLift
Fig. 8: The success rates of the
RRR , NCPC [5],
TWF [16], and
PhaseLift [10] algorithmsas a function of m/n for problem (20) with n = 50. A trial was declared successful ifthe error (21) was below 10 − .uniformly from {− , , i, − i } , the resulting measurements are(24) | y ℓ ( p, q ) | = (cid:12)(cid:12)(cid:12)(cid:12) M − X s =0 N − X t =0 ¯ d ℓ ( s, t ) ρ ( s, t ) e − πips/M e − πiqt/N (cid:12)(cid:12)(cid:12)(cid:12) , ℓ = 1 , . . . , L, where ρ is the M × N signal to be reconstructed. The integer L determines theredundancy of the measurements, corresponding to the ratio m/n in the paragraphsabove. By combining indices, { ℓ, p, q } → µ and { s, t } → ν , the entries of the sensingmatrix A are given by(25) A µν = ¯ d ℓ ( s, t ) e − πips/M e − πiqt/N . In order to test the
RRR algorithm on CDP problems, we have to efficientlyimplement the projection P ( y ) = AA † y . This is done by exploiting the fast Fouriertransform (FFT). From equation (25), we can readily obtain the relations( Aρ ) µ = F (cid:8) ¯ d ℓ ( s, t ) ρ ( s, t ) (cid:9) ( p, q )(26) ( A ∗ y ) ν = L X l =1 M − X p =0 N − X q =0 d ℓ ( s, t ) e πips/M e πiqt/N y ℓ ( p, q )= M N L X l =1 d ℓ ( s, t ) F − (cid:8) y ℓ ( p, q ) (cid:9) ( s, t ) , (27)where F and F − denote the 2-D FFT and its inverse . Moreover, the product A ∗ A Here we use the asymmetric normalization convention for the discrete Fourier transform to beconsistent with the FFT codes in
Matlab . The unitary convention was adopted in the main text.ENCHMARK PROBLEMS FOR PHASE RETRIEVAL A ∗ A ) νν = M N L X l =1 | d ℓ ( s, t ) | . The projection P ( y ) = AA † y = A ( A ∗ A ) − A ∗ y can therefore be implemented viapairs of FFTs, without explicit construction of the matrices.Following tradition, we demonstrate CDP phase retrieval on photographic im-ages . Our test of RRR used only L = 3 masks, where previous work [16] used L = 12. For the ground truth signal we used three 320 × RRR image A † P ( y ) after 1 , , RRR iterations is 5 . × − . As can be seen, RRR needs only a few iterations to get a fineestimation of the signal. Figure 10 gives the reconstruction error (21) as a functionof the iterations. CDP reconstructions with as few as a single random phase mask,using
HIO or the Douglas-Rachford algorithm, have been demonstrated when priorconstraints on the signal are also used [25, 15]. This is not the best practice, because the Fourier power in photographs is artificially concen-trated along the lines p = 0 and q = 0 due to image termination (aperiodicity). VEIT ELSER, TI-YEN LAN AND TAMIR BENDORY (a)(b)(c)(d)(e)
Fig. 9: Reconstruction of the Cornell University image using the
RRR algorithm onthe CDP model (24) with L = 3. The top panel (a) is the initial guess (StanfordUniversity) and the following panels (top to bottom) show the convergence of thereconstruction after (b) 1 (c) 3 (d) 5 (e) 250 RRR iterations.
ENCHMARK PROBLEMS FOR PHASE RETRIEVAL − − − − − e rr o r iterationiteration