## Abstract

In this paper, we study the success rate of the reconstruction of objects of finite extent given the magnitude of its Fourier transform and its geometrical shape. We demonstrate that the commonly used combination of the hybrid input output and error reduction algorithm is significantly outperformed by an extension of this algorithm based on randomized overrelaxation. In most cases, this extension tremendously enhances the success rate of reconstructions for a fixed number of iterations as compared to reconstructions solely based on the traditional algorithm. The good scaling properties in terms of computational time and memory requirements of the original algorithm are not influenced by this extension.

© 2012 Optical Society of America

Inverse problems play an important role in many areas of physics. In an inverse problem we intend to reconstruct physical properties (“scattering potential”) of some object given the result of the interaction (“scattering”) with a well-defined input signal. Typically, propagating wave fields (e.g., acoustic waves, electromagnetic waves, . . .) are used to scan the object. However, for rapidly oscillating electromagnetic wave fields like x-rays, detection devices are only capable of recording the average intensity of the scattered wave field. Its phase information, which also contains important information about the scattering potential, is lost during the measurement process. Despite the lost phase information, x-rays are an important, *non-destructive* tool for the investigation of structural properties of matter because they exhibit high penetration depth and high spatial resolution [1–3]. Note, that 3^{rd} generation synchrotrons [4] provide brilliant x-ray beams with a very high degree of coherence and with very high flux, so that even dot- and wire-like nanostructures can be investigated [2,5–11].

In this paper, we introduce a new algorithm – the (HIO+O_{R})+ER-algorithm – which is capable of reconstructing the missing phase information in a variety of cases for which the traditionally and commonly used HIO+ER-algorithm [12, 13] (reviewed in sec. 2) fails. Before we introduce this new algorithm in sec. 2, we shortly review the phase retrieval problem in sec. 1. In particular, we want to stress the fact that the (HIO+O_{R})+ER-algorithm is based on randomized overrelaxation [14,15] and, thus, does not require additional mathematical constraints [16] compared to the HIO+ER-algorithm. We reformulate our approach using projection polynomials which facilitate a straight forward investigation of randomization of the HIO-algorithm without correlations in the coefficients enforced by overrelaxation. Section 3 discusses numerical aspects how to measure the convergence of phase retrieval. Finally, sec. 4 illustrates the advantages of the (HIO+O_{R})+ER-algorithm with respect to the traditional HIO+ER-algorithms in detail. For this purpose, we focus on three numerical examples each with different specific properties and compare the success rate of reconstructions based on our (HIO+O_{R})+ER-algorithm and on the HIO+ER-algorithm. Moreover, we investigate the sensitivity of our algorithm to the choice of its parameters. In this paper, we focus on the algorithm itself and its improved convergence properties for specific, well-known problems related to phase retrieval based on the HIO+ER-algorithm. A separate publication will illustrate the improved features of the (HIO+O_{R})+ER-algorithm for the reconstruction of the displacement field in inhomogeneously strained nanocrystals.

## 1. The phase retrieval problem

Consider some unknown object *f*_{in}(*x⃗*) in *d*-dimensional direct space (“position space”). We assume that (i) the shape 𝕊 of *f*_{in}(*x⃗*) (i.e., the smallest domain in direct space for which *f*_{in}(*x⃗*) ≠ 0) and (ii) the quantity *A _{k⃗}* ≡ |

**FT**{

*f*

_{in}(

*x⃗*)}| (i.e., the amplitudes of its Fourier transform) are known. It has been proven that this information is sufficient to reconstruct the object

*f*

_{in}(

*x⃗*)

*if*the shape 𝕊 is

*finite*and the dimension

*d*of direct space is at least two [2, 17–19]. Physically, the second condition is fulfilled in scattering experiments which can be described in lowest order Born-approximation (also called kinematical approximation). This is typically true for nanostructures illuminated by coherent x-ray beams.

Suppose we discretize a rectangular region 𝔻 ⊂ ℝ* ^{d}* in direct space, which fully contains our object

*f*

_{in}(

*x⃗*) (i.e., 𝕊 ⊂ 𝔻), with

*N*,

_{i}*i*∈ {1, 2,...,

*d*}, equidistant points. The discrete Fourier transform (DFT) maps this set of points to a mesh 𝕄 in reciprocal space which is related to the Fourier transform of our continuous direct space object

*f*

_{in}(

*x⃗*) [20]. Hence, we have

*x⃗*= (

*x*

_{1},

*x*

_{2},...,

*x*)

_{d}*corresponds to the equidistantly sampled points in direct space defined by 𝔻 and Φ*

^{T}*are the phases associated with the amplitudes*

_{k⃗}*A*in reciprocal space.

_{k⃗}In direct space, every point *x⃗* inside the object domain 𝕊 corresponds to (up to) two real unknowns (magnitude and phase). However, outside the domain 𝕊, *f*_{in}(*x⃗*) is precisely zero for all positions *x⃗*. Hence, each point outside 𝕊 defines two real equations for the determination of the *N* real phases Φ* _{k⃗}*, if we require

*f*

_{in}(

*x⃗*) = 0 ∀

*x⃗*∉ 𝕊. Therefore, we can hope to reconstruct the missing phases Φ

*on the sampled grid 𝕄 once*

_{k⃗}*N*

_{𝕊}is the number of points inside the region 𝕊. The fraction

*σ*=

*N/N*

_{𝕊}is called oversampling ratio. From our elementary discussion, we can conclude that

*σ*= 2 is a lower bound for a successful reconstruction, if no additional a priori knowledge is available. For more details on oversampling, we refer the reader to the investigations of Elser et Milliane [21] and of Miao et al. [22] (and the references therein).

With that knowledge, we restate the phase retrieval problem more precisely: Given the shape 𝕊 of some object *f*_{in}(*x⃗*) and the modulus *A _{k⃗}* of the Fourier transform of

*f*

_{in}(

*x⃗*) on a mesh 𝕄 that satisfies Eq. (2), we try to reconstruct the phases Φ

*and, thus, the object*

_{k⃗}*f*

_{in}(

*x⃗*).

Note, that the reconstruction is only unique up to the unavoidable inherent symmetries of the Fourier transform [2, 17, 18], i.e., a global phase shift (typically irrelevant), a plane wave modulation *e*^{ik⃗·x⃗} in reciprocal space (fixed by the position of the shape 𝕊 inside 𝔻) and the degeneracy between the object *f _{x⃗}* and its complex conjugated, inversion symmetric object
${f}_{\overrightarrow{x}}^{*}$. This last ambiguity constitutes a severe problem if the shape 𝕊 is inversion symmetric, but the object

*f*

_{in}does not possess this symmetry [13,23]. For the remainder of this paper, we restrict to domains 𝕊 which are not inversion symmetric. However, even if we exclude objects with inversion symmetric shape, it is highly non-trivial to formulate a phase retrieval algorithm that works

*automatically*without manual supervision and tuning of parameters during the reconstruction. The (HIO+O

_{R})+ER-algorithm which we propose in the next section is capable of performing this task for a large class of objects

*f*

_{in}.

## 2. Reconstruction algorithms

In this section, we will shortly review the traditional combination of the hybrid input output and error reduction (HIO+ER) algorithm [12, 13] and propose our extension to this algorithm. Throughout the entire reconstruction process, we assume that we know the exact geometrical shape 𝕊 of the object which we reconstruct.

#### 2.1. Hybrid input output and error reduction

A fast and efficient iterative phase retrieval algorithm is based on alternating projections of a trial solution onto the constraints in direct space and in reciprocal space. Hence, a single iterative step (*i*) of this algorithm, which is called error reduction (ER) algorithm [12], is defined by

**P**

_{𝕊}and

**P**

*are projection operators in direct space and reciprocal space respectively, i.e.,*

_{A}**P**

*or*

_{A}**P**

_{𝕊}is transformed to the proper space (i.e., Fourier transform to reciprocal or direct space respectively) before the operator

**P**

*or*

_{A}**P**

_{𝕊}is applied (e.g., ${\mathbf{P}}_{A}{f}_{\overrightarrow{x}}^{(i)}$ means ${\mathbf{P}}_{A}\mathbf{FT}\left\{{f}_{\overrightarrow{x}}^{(i)}\right\}$).

The projection operator **P*** _{A}* is

*non-linear*,

*non-convex*and

*non-unique*[24]. The values $\left|{g}_{\overrightarrow{k}}^{(i)}\right|$ are important for determining the difference to the (experimentally accessible) input data

*A*(see Eq. (18)). The ER-algorithm is a local minimizer of a suitable chosen error metrics [12,13]. However, in practice, phase retrieval problems involve many local minima. Hence, the ER-algorithm will in general not converge to the correct solution

_{k⃗}*f*

_{in}(

*x⃗*), but to some local minimum. In addition, the error metric may stagnate for many iterations before decreasing further. More information on stagnation problems and the ER-algorithm can be found in [12,13,23].

Due to these problems, further algorithms, which try to avoid stagnation and convergence to local minima, have been developed. A very important algorithm is the hybrid input output (HIO) algorithm proposed by Fienup [12, 13, 23]. This algorithm is also an iterative procedure which is defined by the map

*β*is a real parameter typically chosen in the range [0.5; 1.0] [23]. Note, that this definition can be rewritten as

The convergence properties of the HIO-algorithm have been investigated in more detail in [22,25,26].

This algorithm is much stronger in avoiding stagnation and local minima. However, the HIO-algorithm shows its full potential only if it is combined with the ER-algorithm. One step of this combined HIO+ER-algorithm consists of *N*_{HIO} iterations of the HIO-algorithm followed by *N*_{ER} iterations of the ER-algorithm. This combination is more successful than both algorithms on their own as it was observed in practice [13, 23]. Although this algorithm is already very powerful, it is not yet satisfactory for many objects as we will show in sec. 4.

#### 2.2. Hybrid input output with randomized overrelaxation ((HIO+O_{R})+ER)

The HIO+ER-algorithm can be further improved by including randomized overrelaxation [14, 15,27]. Basically, overrelaxation corresponds to replacing a projection operator **P*** _{μ}* by

*λ*is a real constant called

_{μ}*relaxation parameter*[14]. The limiting case

**Q**

*→*

_{μ;λμ}**P**

*is obtained for*

_{μ}*λ*→ 1.

_{μ}Overrelaxation (without any randomization) has been investigated for convex problems [15] and in connection with the ER-algorithm for phase retrieval [14]. Moreover, overrelaxation is also included in the difference map algorithm proposed by Elser in [27]. In the difference map algorithm with overrelaxation, the iterative step for our set of constraints is given by [28]

*λ*=

_{A}*λ*

_{𝕊}=

*β*

^{−1}[27]. A very nice comparison of several phase retrieval algorithms and a benchmark thereof can be found in [28]. Marchesini demonstrated in [29] that for his particular numerical example an additional low-dimensional subspace saddle-point optimization was also able to overcome stagnation of the traditional HIO+ER-algorithm. In our extension, such an additional optimization is not necessary.

We propose to replace the (non-linear and non-convex) projection operator **P*** _{A}* in reciprocal space in the HIO-algorithm itself by its relaxed expression

**P**

_{𝕊}and

**P**

*, we get*

_{A}The deviation *β* (*λ _{A}* − 1) from the identity operator in the first term can neither be represented by the traditional HIO-algorithm for any value of

*β*(see Eq. (6)) nor by the difference map algorithm for any (

*β*,

*λ*,

_{A}*λ*

_{𝕊}) (see Eq. (8)). In both cases, the previous iterative result ${f}_{\overrightarrow{x}}^{(i)}$ is weighted with 1 or projected at least once in either direct (

**P**

_{𝕊}) or reciprocal space (

**P**

*) before being included in the calculation for ${f}_{\overrightarrow{x}}^{(i+1)}$.*

_{A}Finally, we have to choose suitable values for the additional parameter *λ _{A}* in Eq. (9). If we restrict to a fixed set of values

*λ*for all iterations, we risk simply exchanging one trap by another or one local minimum by another. Trying to overcome local minima and stagnation, we propose a new approach: For each iteration, we

_{A}*randomly*select the parameter

*λ*. More precisely, for the remainder of this paper, we investigate a uniform random distribution in the range [1 −

_{A}*ν*, 1 +

*ν*],

*ν*≥ 0, for

*λ*which is reassigned each iteration. Unless stated otherwise, we choose

_{A}*ν*= 0.5. Note, that a fixed value

*λ*≡ 1 for all iterations corresponds to the usual HIO-algorithm. In order to distinguish our new extension from the traditional HIO-algorithm, we call our extension HIO+O

_{A}_{R}-algorithm. The power of this extension is strongly supported by the fact, that the HIO+O

_{R}-algorithm is capable of reconstructing objects without including the ER-algorithm with significant success rates. If the HIO+O

_{R}-algorithm is combined with ER, we call the algorithm (HIO+O

_{R})+ER-algorithm (see Fig. 1).

If we look at Eq. (11) from a bird’s eye view, the most important difference of the operator *Ĥ*_{HIO+OR} (*β*, *λ _{A}*) to the other algorithms which have been proposed [28] is the fact, that the traps and tunnels [14] of the operator in general depend on

*λ*and, therefore, should change from iteration to iteration as a consequence of randomization. Thus, stagnation in a specific trap or tunnel is strongly reduced, as most traps and tunnels are no longer persistent throughout the iterative reconstruction process. However, the true solution

_{A}*f*

_{in}(

*x⃗*) is a fixed point of the iterative procedure defined in Eq. (11) for

*all*values of

*λ*.

_{A}#### 2.3. Projection polynomials

The extension of the HIO+ER-algorithm to the (HIO+O_{R})+ER-algorithm is based on two ingredients: overrelaxation and randomization. In order to gain a deeper understanding of the (HIO+O_{R})+ER-algorithm, it is useful to investigate both ingredients separately. Studying the HIO-algorithm extended by overrelaxation of **P*** _{A}*, but without randomization is straight forward: we simply keep a fixed, predefined overrelaxation parameter

*λ*throughout the entire iterative reconstruction process. In this section, we present a framework for randomization of iterative projection algorithms performed in a way which is not based on overrelaxation of

_{A}**P**

*. Although this framework can be applied for studying randomization of any iterative projection algorithm, we limit our numerical simulations (in section 4) to parameter values whose deterministic contribution coincides the traditional HIO-algorithm.*

_{A}For this purpose, consider the projection polynomial operator

*p*is the maximum number of successive projection operators which should be included in

*Ĥ*

_{Proj}.

This operator exploits the fact, that *any* combination of projection operators
${\mathbf{P}}_{{\xi}_{1}}^{{n}_{1}}{\mathbf{P}}_{{\xi}_{2}}^{{n}_{2}}\dots {\mathbf{P}}_{{\xi}_{m}}^{{n}_{m}}$ of two different kinds *ξ _{j}* ∈ {𝕊,

*A*},

*j*∈ {1,...,

*m*}, reduces to one of the four building blocks

**P**

_{𝕊}(

**P**

_{A}**P**

_{𝕊})

*,*

^{n}**P**

*(*

_{A}**P**

_{𝕊}

**P**

*)*

_{A}*, (*

^{n}**P**

_{𝕊}

**P**

*)*

_{A}*and (*

^{n}**P**

_{A}**P**

_{𝕊})

*for some integer*

^{n}*n*≥ 0. This is a consequence of the defining property ${\mathbf{P}}_{\xi}^{2}={\mathbf{P}}_{\xi}$ (idempotence) of projection operators (which implies ${\mathbf{P}}_{\xi}^{n}={\mathbf{P}}_{\xi}\forall n\ge 1$). Note, that the product of two non-commutative idempotent operators (like

**P**

*and*

_{A}**P**

_{𝕊}) is no longer idempotent. The special case of the identity operator 1 and a single projection operator is included for the case

*n*= 0 in those building blocks. However, the identity operator is included twice, namely for

*n*= 0 for (

**P**

_{𝕊}

**P**

*)*

_{A}*and (*

^{n}**P**

_{A}**P**

_{𝕊})

*. Therefore, one spurious coefficient in the projection polynomial is introduced if both building blocks are included for*

^{n}*n*= 0. This is avoided by restricting those two building blocks to

*n*≥ 1 and including the identity operator separately. Therefore, Eq. (12a) represents the most general polynomial expression which contains a maximum of

*q*successive projections of two different kinds, if we incorporate the idempotence of projection operators.

In order to guarantee, that the true solution *f*_{in}(*x⃗*) is a fixed point of *Ĥ*_{Proj} for any choice of its parameters (*b*, *c⃗*_{𝕊}, *c⃗ _{A}*), it is necessary to enforce one constraint on the coefficients: If the input for the projection operators

**P**

_{𝕊}and

**P**

*coincides with the true solution*

_{A}*f*

_{in}(

*x⃗*), the projection operators

**P**

_{𝕊}and

**P**

*reduce to the identity operator 1. Therefore, given the true solution*

_{A}*f*

_{in}(

*x⃗*) as input, the operator

*Ĥ*

_{Proj}simplifies to the identity operator 1 for any choice of its parameters (

*b*,

*c⃗*

_{𝕊},

*c⃗*)

_{A}*if and only if*the additional constraint

*p*successive projections in the projection polynomial

*Ĥ*

_{Proj}, only 2

*p*free parameters appear in

*Ĥ*

_{Proj}by incorporating the idempotence of projection operators.

Given this fix point property for the true solution *f*_{in}(*x⃗*), we can perform proper randomization of this operator. For this purpose, we split our coefficients *c*_{ξ,n}, *ξ* ∈ {𝕊, *A*}, *n* ≥ 1, in a deterministic part
${c}_{\xi ,n}^{(D)}$ and a random part
${r}_{\xi ,n}{c}_{\xi ,n}^{(R)}$ and assume *r _{ξ}_{,n}* to be uniformly distributed in the range [−1, 1], i.e.,

*r*, we can investigate randomization of the HIO-algorithm without the correlations in the coefficients

_{ξ}_{,n}*c*introduced by overrelaxation.

_{ξ,n}*b*is fixed by Eq. (12b) and, thus, is fully correlated to the values of

*c*.

_{ξ,n}In order to obtain the traditional HIO-algorithm in the limit
${c}_{\xi ,n}^{(R)}\to 0$ for all *ξ* and *n*, we have to choose (at least) *q* = 2 (see Eq. (6)). Therefore, *Ĥ*_{Proj} reduces to

**P**

*by*

_{A}**Q**

_{A;λA}(see Eq. (9)), the coefficients

*c*are parametrized by two parameters

_{ξ}_{,n}*β*and

*λ*as

_{A}*λ*= 1 +

_{A}*γ*.

_{A}*γ*is uniformly distributed in [−

_{A}*ν*,

*ν*]. The deterministic contribution in Eq. (16b) is identical to Eq. (15) and reproduces the traditional HIO-algorithm in the limit

*ν*→ 0. However, the random contribution of each coefficient

*c*is determined solely by the value of

_{ξ}_{,n}*γ*and, thus, not statistically independent, but fully correlated. This is in strong contrast to a statistically independent randomization of each coefficient

_{A}*c*separately.

_{ξ,n}Consequently, we can apply the framework of projection polynomials to investigate whether the benefits of randomization rely on the precise correlations evident in Eq. (16b) or if randomization of the traditional HIO-algorithm implies significant benefits also in absence of these correlations.

## 3. Monitoring convergence of the reconstruction

In order to investigate the success and convergence properties of the reconstruction, we monitor three parameters. First of all, we monitor the change of the reconstructed object *f*^{(i)} from iteration (*i* − 1) to iteration (*i*). Second, we measure the mathematical distance to the (sampled) true solution *f*_{in}(*x⃗*). And last, but not least, we measure the deviation in reciprocal space between the a priori given reciprocal input data *A _{k⃗}* and the magnitude of the Fourier transform of the reconstructed object. All three parameters yield different information:

The distance of the reconstructed object *f*^{(i)} after iteration (*i*) to the true direct space object *f*_{in} is, of course, the best measure for the quality of the reconstructed image. We define this mathematical distance by an angle

*f*

^{(i)}. Moreover,

*φ*

^{(i)}is identical in direct and reciprocal space (invariance of scalar product upon unitary basis transformations like the Fourier transform).

A dimensionless, normalized error measure *ε*^{(i)} which depends only on the a priori known (measured) amplitudes *A _{k⃗}* and not on the true direct space solution

*f*

_{in}(

*x⃗*) is

*is the*

_{p}*p*-norm. Since the construction of the absolute value is a non-linear, non-invertible map, the result can in general be quite different from the angle defined above. In fact, two quite different direct space objects

*f*

_{in}(hence quite different complex Fourier transforms) may posses an extremely similar distribution of the magnitude

*A*in reciprocal space [16,18].

In addition, from a numerical point of view, the change of the reconstructed object *f*^{(i)} from iteration to iteration is very important, too. This must not be confused with the change of the non-invertible map *ε*^{(i)} from iteration to iteration: Moving along a connected line of constant error metric does not change the error metric itself, but the reconstructed object *f*^{(i)} may change arbitrarily. Nevertheless, the algorithm may converge (i.e. no change from iteration to iteration), but not to the true solution *f*_{in} of the problem [14, 16]. For simulated data, we can judge if the algorithm converged to the true solution *f*_{in} by Eq. (17). A suitable choice for observing the change from iteration to iteration is given by the angle

*p*-norm ${\delta}^{(i)}={\Vert {f}^{(i-1)}-{f}^{(i)}\Vert}_{p}$ of the difference is not optimal for measuring the convergence of the reconstruction process, because it does not eliminate global phase shifts from iteration (

*i*− 1) to iteration (

*i*).

## 4. Numerical examples

In this section, we demonstrate the improvements resulting from randomized overrelaxation for three carefully chosen numerical examples: First, we consider two pure phase objects

*ξ*(

*x⃗*) is a

*real*function and Ω

_{𝕊}(

*x⃗*) is the shape function

*f*

_{in}which results in severe problems for the traditional HIO+ER-algorithm. The reconstruction of the second phase object (see Fig. 2(b)) is plagued by many local minima near the perfect solution

*f*

_{in}, in particular, by stripes [23]. Both phase objects have a smooth variation of their phase, i.e., phase variations of 2

*π*extend over several sampling points. As a third test object, we investigate a purely real object

*f*

_{in}(

*x⃗*) with strong variation of its magnitude over very short length scales and a different shape (see Fig. 3(a)). In fact, the success rate of the traditional HIO-algorithm without randomized overrelaxation was lowest for this third example (see Fig. 3(b) and Fig. 4). From the numerical results, we conclude that randomized overrelaxation is a powerful extension for a wide class of objects and not only for objects possessing some particular features.

The first pure phase objects which we used in our numerical simulations is given by

*ξ*

_{1}is defined according to Eq. (20) and the parameters were chosen as

*b*

_{1}= 1.5515 and

*c*

_{1}= −1.835. The shape 𝕊 has been restricted to the triangular region defined by (0, 0), (0,

*a*) and (

*a*, 0) with

*a*= 0.97. The object has been sampled on an equidistant grid with

*N*

_{1}×

*N*

_{2}= 256 × 256 data points in the interval (

*x*,

*y*) = (−1, −1) to (

*x*,

*y*) = (1, 1). The resulting object is depicted in Fig. 2(a).

The second object is defined by

*b*

_{2}= −0.3515 and

*c*

_{2}= 0.535. The phase field in this case is depicted in Fig. 2(b).

A reconstruction is considered successful once the angle *φ*^{(i)} to the input object (see Eq. (17)) falls below some given limit *φ*_{Converged} (=0.05° unless stated otherwise). We repeat the reconstruction process for *N*_{Real} initial trials (random phases at each *k⃗*-point, amplitudes
${A}_{\overrightarrow{k}}^{(0)}$ equal to given amplitudes *A _{k⃗}*). Finally, we calculate the success rate which tells us the percentage of reconstruction processes that have been successful up to iteration (

*i*). A good reconstruction algorithm

- should reach a success rate of almost 100%,
- should not depend on the starting point (as long as no good starting point is available),
- should perform the reconstructions with little computational effort and
- should possess these properties for a wide range of its internal parameters (i.e.,
*β*,*ν*,*N*_{HIO}and*N*_{ER}in case of the (HIO+O_{R})+ER-algorithm).

Each iteration in the HIO+ER-algorithm and its (HIO+O_{R})+ER-extension can be implemented very efficiently exploiting the good scaling properties of the FFT-algorithm. Therefore, as long as the number of iterations is not too high, the third requirement is clearly met by both, the HIO+ER and (HIO+O_{R})+ER-algorithm.

In order to investigate the other points, we plotted the respective success rates for the HIO+ER- and (HIO+O_{R})+ER-algorithm for different internal parameters in Fig. 4. We normalized the computational effort, i.e., one iteration of the (HIO+O_{R})+ER-algorithm with *N*_{HIO} = 130 and *N*_{ER} = 10 corresponds to two successive iterations with *N*_{HIO} = 50 and *N*_{ER} = 20. This prefactor of 2 can be found in the legend of the graph in front of the square brackets for each calculation.

The inclusion of overrelaxation improves the success rate for the first phase object defined by Eq. (22) for each combination of *N*_{HIO} and *N*_{ER}: in contrast to the HIO based calculations, it quickly converged to 100%. In fact, overrelaxation is so powerful that stagnation in the HIO-algorithm is eliminated entirely and intermediate ER-iterations are no longer necessary (see black dash-dotted curve in Fig. 4). The success rate of the traditional HIO+ER-algorithm stagnates on some level far below 100%. This stagnation is significant in the long term behavior, as we illustrate in Fig. 5 by plotting the success rate of the traditional HIO+ER-algorithm for 50 times as many iterations as in Fig. 4(a). Note, that all trials that did not converge after 2500 iterative steps using the traditional HIO+ER-algorithm with *N*_{HIO} = 130 and *N*_{ER} = 10 in Fig. 5 are separated from the true solution by an angle greater than *φ* = 55°.

If we consider the second phase object defined by Eq. (23), the traditional HIO+ER-algorithm performs better here than for the first phase object. For the one specific set of parameters *N*_{HIO} = 130 and *N*_{ER} = 10, the success rate of the traditional HIO+ER-algorithm is even slightly superior to the (HIO+O_{R})+ER-extension, but for this specific set of parameters, the success rate converges to 100% very quickly for both algorithms. However, for the set of parameters (*N*_{HIO} = 50, *N*_{ER} = 20) and (*N*_{HIO} = 40, *N*_{ER} = 30), the success rate quickly converges to 100% only for the (HIO+O_{R})+ER-algorithm, but not for the HIO+ER-algorithm. Hence, the (HIO+O_{R})+ER-algorithm is much more stable with respect to the particular choice of the internal parameters *N*_{HIO} and *N*_{ER} including, in particular, the case *N*_{ER} = 0.

The traditional HIO+ER-algorithm also has severe problems performing a successful reconstruction for the real object depicted in Fig. 3(a): The success rate for the reconstruction process is very poor if we apply our criterion *φ*_{Converged} ≤ 0.05° to determine successful reconstructions. The (HIO+O_{R})+ER-algorithm is much more successful and reaches success rates of 100% with few iterations (see Fig. 3(b)). Furthermore, reaching the success rate of 100% does not sensitively depend on the internal parameters *N*_{HIO} and *N*_{ER}. Again, keeping *N*_{ER} small compared to *N*_{HIO} or eliminating ER-iterations completely yields the best results. Keep in mind, that neither the reality constraint nor the positivity constraint of the object have been exploited during reconstruction.

For the second and third test object, all trials that failed to converge to the true solution *f*_{in} basically evolved towards the correct direction. However, at some point, stagnation in form of stripes close to the true solution *f*_{in} appeared. These stripes are very persistent and prevent the usual HIO+ER-algorithm from further convergence to the true solution *f*_{in}. Typically, the mathematical distance *φ* of the objects containing stripes to the true solution *f*_{in} is in the order of 0.05 to 5.0 degrees. In case of the purely real test object, after 150 iterations all reconstructed objects reached a distance *φ* ≤ 0.5° (= 10 · *φ*_{Converged}) for HIO_{130} − ER_{10} and 31% for HIO_{50} − ER_{20}. The remaining 71% of the objects for HIO_{50} − ER_{20} reached a distance *φ* ≤ 5° (= 100 · *φ*_{Converged}) using traditional HIO+ER-algorithm and 150 iterations. Note, that the change *χ*^{(}^{i}^{)} (see Eq. (19)) from iteration to iteration does not vanish, i.e. numerically, the algorithm did not converge within the given number of iterations.

Figure 6(a) shows the behavior of the success rate of the (HIO+O_{R})+ER-algorithm for different values of the parameter *ν*. In addition, Fig. 6(b) depicts the success rate for different values of the parameter *β*. Within the range *β* ∈ [0.5, 1.0], which is typically used in the traditional HIO+ER-algorithm, the (HIO+O_{R})+ER-extension does not sensitively depend on the particular choice of *β*. Therefore, the (HIO+O_{R})+ER-algorithm does not depend sensitively on any of its internal parameters *N*_{HIO}, *N*_{ER}, *β* and *ν*. Hence, it fulfills all four requirements on a good reconstruction algorithm which we stated above.

Two concepts (overrelaxation and randomization) have been exploited to modify the HIO-algorithm: Fig. 7 depicts the success rate, if overrelaxation is performed with a static parameter *λ _{A}*. We see extremely sensitive behavior of the success rate on the precise value of

*λ*if it is kept fixed. Whereas increasing

_{A}*λ*up to 1.02 improved the success rate for the first object, it completely vanished for

_{A}*λ*= 1.01 for the second phase object. For a deviation from 1.00 greater than ten percent, almost no successful reconstructions have been observed for any of the three test objects. In addition, we do not a priori know the (non-universal) optimum value for

_{A}*λ*which depends on the object

_{A}*f*

_{in}. Consequently, static overrelaxation does not fulfill the requirements for a good reconstruction algorithm stated above.

Based on the discussion in sec. 2.3, we test randomization of the HIO-algorithm for *β* = 0.85 without employing overrelaxation by setting the amplitudes
${c}_{\xi ,n}^{(R)}$ for the random contribution to the coefficients *c _{ξ,n}* in the operator

*Ĥ*

_{HIOR}defined in Eq. (14) equal to ${c}_{\mathbb{S},1}^{(R)}={c}_{A,1}^{(R)}={c}_{\mathbb{S},2}^{(R)}={c}_{A,2}^{(R)}=R=0.2$. The deterministic contribution ${c}_{\xi ,n}^{(D)}$ is set according to Eq. (15). The result for all three test objects is illustrated in Fig. 8 and compared to the result for the HIO+O

_{R}-algorithm. Clearly, pure randomization without overrelaxation performs equally well as the HIO+O

_{R}-algorithm for all three test objects. However, in our opinion, the HIO+O

_{R}-algorithm should be preferred over randomization without overrelaxation, because (i) it requires only two Fourier transformations for each iteration (lower computational effort) and (ii) it has less degrees of freedom (one uniform random distributions instead of four). Note, however, that due to the large number of degrees of freedom of projection polynomials, further optimization of this approach is possible (including more advanced statistical correlations in the coefficients of the projection polynomial). In addition, more elaborate random distributions might improve both approaches further.

In summary, we succeeded in overcoming stagnation of the traditional HIO+ER-algorithm for various objects. For this purpose, the concept of overrelaxation was applied to the reciprocal space projection in the HIO-algorithm which produced additional degrees of freedom. Randomization of these additional degrees of freedom from iteration to iteration was exploited to prevent stagnation in local minima, especially in traps or tunnels. In particular, this was demonstrated for two pure phase objects with smooth phase variation (with different behavior with respect to the traditional HIO+ER-algorithm) and a purely real object with strong amplitude variations over short length scales. The improved algorithm is far less sensitive to the specific choice of the internal parameters *N*_{HIO} and *N*_{ER} than the traditional HIO+ER-algorithm. Moreover, the (HIO+O_{R})+ER-algorithm remains almost insensitive to the particular value of the internal parameter *β* in the range *β* ∈ [0.5, 1.0]. Furthermore, the the good scaling properties in terms of computational time and memory consumption of the HIO+ER-algorithm are fully maintained in the (HIO+O_{R})+ER-algorithm. In conclusion, a reconstruction based on the (HIO+O_{R})+ER-algorithm exhibits significantly enhanced stability and success probability in comparison with the traditional and commonly used HIO+ER-algorithm.

## References and links

**1. **F. v. d. Veen and F. Pfeiffer, “Coherent x-ray scattering,” Phys J..: Condens. Matter **16**, 5003–5030 (2004). [CrossRef]

**2. **R. Millane, “Phase retrieval in crystallography and optics,” J. Opt. Soc. Am. A **7**, 394–411 (1990). [CrossRef]

**3. **U. Pietsch, V. Holy, and T. Baumbach, *High-Resolution X-ray Scattering From Thin Films to Lateral Nanostructures* (Springer, New York, 2004).

**4. **G. N. Afanasiev, *Vavilov-Cherenkov and Synchrotron Radiation: Foundations and Applications* (Springer, Netherlands, 2004).

**5. **I. Robinson and R. Harder, “Coherent x-ray diffraction imaging of strain at the nanoscale,” Nat Mater **8**, 291–298 (2009). [CrossRef] [PubMed]

**6. **J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of x-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature **400**, 342–344 (1999). [CrossRef]

**7. **M. A. Pfeifer, G. J. Williams, I. A. Vartanyants, R. Harder, and I. K. Robinson, “Three-dimensional mapping of a deformation field inside a nanocrystal,” Nature **442**, 63–66 (2006). [CrossRef] [PubMed]

**8. **A. Biermanns, A. Davydok, H. Paetzelt, A. Diaz, V. Gottschalch, T. H. Metzger, and U. Pietsch, “Individual GaAs nanorods imaged by coherent x-ray diffraction,” J. Synchrotron Radiat. **16**, 796–802 (2009). [CrossRef] [PubMed]

**9. **A. A. Minkevich, M. Gailhanou, J.-S. Micha, B. Charlet, V. Chamard, and O. Thomas, “Inversion of the diffraction pattern from an inhomogeneously strained crystal using an iterative algorithm,” Phys. Rev. B **76**, 104106 (2007). [CrossRef]

**10. **A. A. Minkevich, E. Fohtung, T. Slobodskyy, M. Riotte, D. Grigoriev, T. Metzger, A. C. Irvine, V. Novák, V. Holý, and T. Baumbach, “Strain field in (Ga,Mn)As/GaAs periodic wires revealed by coherent x-ray diffraction,” Europhys. Lett. **94**, 66001 (2011). [CrossRef]

**11. **A. A. Minkevich, E. Fohtung, T. Slobodskyy, M. Riotte, D. Grigoriev, M. Schmidbauer, A. C. Irvine, V. Novák, V. Holý, and T. Baumbach, “Selective coherent x-ray diffractive imaging of displacement fields in (Ga,Mn)As/GaAs periodic wires,” Phys. Rev. B **84**, 054113 (2011). [CrossRef]

**12. **J. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. **21**, 2758–2769 (1982). [CrossRef] [PubMed]

**13. **J. Fienup, “Reconstruction of a complex-valued object from the modulus of its Fourier transform using a support constraint,” J. Opt. Soc. Am. A **4**, 118–123 (1986). [CrossRef]

**14. **A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A **1**, 932–943 (1984). [CrossRef]

**15. **D. C. Youla and H. Webb, “Image restoration by the method of convex projections: part 1 – theory,” IEEE Trans. Med. Imaging **1**, 81 –94 (1982). [CrossRef] [PubMed]

**16. **A. A. Minkevich, T. Baumbach, M. Gailhanou, and O. Thomas, “Applicability of an iterative inversion algorithm to the diffraction patterns from inhomogeneously strained crystals,” Phys. Rev. B **78**, 174110 (2008). [CrossRef]

**17. **R. P. Millane, “Multidimensional phase problems,” J. Opt. Soc. Am. A **13**, 725–734 (1996). [CrossRef]

**18. **J. H. Seldin and J. R. Fienup, “Numerical investigation of the uniqueness of phase retrieval,” J. Opt. Soc. Am. A **7**, 412–427 (1990). [CrossRef]

**19. **R. H. T. Bates, “Fourier phase problems are uniquely solvable in more than one dimension,” Optik (Stuttgart) **61**, 247–262 (1982).

**20. **L. Auslander and F. A. Grunbaum, “The Fourier transform and the discrete Fourier transform,” Inverse Prob. **5**, 149 (1989). [CrossRef]

**21. **V. Elser and R. P. Millane, “Reconstruction of an object from its symmetry-averaged diffraction pattern,” Acta Crystallogr., Sect. A: Found. Crystallogr. **64**, 273–279 (2008). [CrossRef]

**22. **J. Miao, D. Sayre, and H. N. Chapman, “Phase retrieval from the magnitude of the Fourier transforms of nonperiodic objects,” J. Opt. Soc. Am. A **15**, 1662–1669 (1998). [CrossRef]

**23. **J. Fienup and C. Wackermann, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A **3** (1986). [CrossRef]

**24. **H. H. Bauschke, P. L. Combettes, and D. R. Luke, “Phase retrieval, error reduction algorithm, and fienup variants: a view from convex optimization,” J. Opt. Soc. Am. A **19**, 1334–1345 (2002). [CrossRef]

**25. **H. Takajo, T. Takahashi, R. Ueda, and M. Taninaka, “Study on the convergence property of the hybrid input–output algorithm used for phase retrieval,” J. Opt. Soc. Am. A **15**, 2849–2861 (1998). [CrossRef]

**26. **H. Takajo, T. Takahashi, and T. Shizuma, “Further study on the convergence property of the hybrid input–output algorithm used for phase retrieval,” J. Opt. Soc. Am. A **16**, 2163–2168 (1999). [CrossRef]

**27. **V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A **20**, 40–55 (2003). [CrossRef]

**28. **S. Marchesini, “Invited article: A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. **78**, 011301 (2007). [CrossRef] [PubMed]

**29. **S. Marchesini, “Phase retrieval and saddle-point optimization,” J. Opt. Soc. Am. A **24**, 3289–3296 (2007). [CrossRef]