Simulated Annealing Applied to HIO Method for Phase Retrieval

Abstract

Phase retrieval utilizing Fourier amplitudes plays a significant role in image recovery. Iterative phase retrieval algorithms have been developed to retrieve phase information that cannot be recorded by detectors directly. However, iterative algorithms face the problem of being trapped in local minima due to the nonconvexity of phase retrieval, and most existing works addressed this by optimizing in multiple runs parallelly to improve the possibility that one of these could reach the global minimum. Alternatively, we propose in this work to increase the probability of reaching the global minimum with one arbitrary initial distribution by adapting simulated annealing in the standard hybrid input-output (HIO) algorithm. Numerical and experimental results demonstrate that the proposed method reconstructs images with mean square errors 50.12% smaller than those reconstructed by HIO. More importantly, the proposed method can be applied to any HIO-based algorithm with multiple runs to further improve the performance.

1. Introduction

The phase retrieval technique plays a crucial part in image recording and reconstruction. Due to the nature of the high oscillation frequency of light, optoelectronic detector arrays cannot follow and record only the magnitude of the light field, i.e., average light intensity over time. In this process, phase, being an essential parameter of the light field and providing significant information of the object, is lost during the measurements. Therefore, efforts have been made for decades to explore different methods for phase retrieval in applications such as holographic imaging [1,2,3], astronomical imaging [4,5,6], and X-ray crystallography [7].

Classic approaches to phase retrieval rely on iterations between Fourier and object domains. The Gerchberg–Saxton (GS) algorithm [8], being one of the earliest successful attempts, performs an iterative procedure to obtain an estimation of the phase. Fienup modified the Gerchberg–Saxton algorithm with a new set of object constraints and an error function, and named the algorithm the hybrid input-output (HIO) algorithm [9]. The HIO method improves converging speed significantly but yields varying results depending on different initials, most of which lead to local minima of the optimization. Therefore, multiple runs need to be performed in HIO parallelly and the operation of getting the average results of those with small errors is essential. To increase the consistency of the results, the geometric average was introduced to the iterations of the HIO algorithm, creating the guided hybrid input-output (GHIO) algorithm [10]. The GHIO method manages to jump out of local minima by running multiple initials parallelly and making good use of a better one between all intermediate results during iterations to guide the direction. However, the success of GHIO depends strongly on multiple runs just like HIO, which come at the cost of more computation. Another multiple initials algorithm utilized particle swarm optimization (PSO) in HIO, which makes random initial particles learn and renew each other, and therefore increases the possibility of reaching the global minimum [11,12]. However, the random nature of stochastic perturbations jeopardized the efficiency of the algorithm.

In the past decade, there have been some new attempts on phase retrieval. Inspired by the idea of compressed sensing [13], Moravec et al. put forward compressive phase retrieval, in which the signal’s sparsity was used as additional prior information to help phase recovery [14]. Different measuring methods, such as tilt illumination [15,16] and optical gratings [17], were introduced to acquire more information for phase retrieval [18].

The emergence of novel methods has significant impacts on the development of phase retrieval, and in the meantime there is still room for improvements of classic algorithms, especially when being combined with modern intelligent optimization algorithms. Here in this work, we propose applying the concept of simulated annealing [19] to HIO to tune the iteration parameter adaptively during the optimization in order to address the problem that many runs need to be performed parallelly in HIO to yield good phase retrieval results. In the process of the algorithm combining HIO and simulated annealing, the present iteration results that own smaller errors in the Fourier domain (ERF) than the previous will be accepted unconditionally; those owning larger errors in the Fourier domain (ERF) than the previous will be accepted according to a simulated annealing probability. Performances of different phase retrieval algorithms are compared quantitatively by numerical simulations and practical experiments. The comparison demonstrates that the proposed method is superior in the converging speed and degree during iterations.

2. Hybrid Input-Output (HIO) Algorithm

In Fienup’s HIO algorithm, phase retrieval is performed by the following steps:

(1)

Generate a random distribution $p_0(x,y)$ in the object domain as an initial distribution. By performing Fourier transform on $p_0(x,y)$, a complex distribution $P_0(u,v)$ is yielded in Fourier domain,

$$P_0(u,v)=F\{p_0(x,y)\}=|P_0(u,v)|\exp [i{\psi }_0(u,v)]$$

(1)

where $|P_0(u,v)|$ is the magnitude, and ${\psi }_0(u,v)$ is the phase;

(2)

Replace the modulus $|P_0(u,v)|$ with the measured magnitudes experimentally known as $A(u,v)$ and keep its original phases to get a combinatorial ${P^{\prime }}_0(u,v)$, i.e.,

$${P^{\prime }}_0(u,v)=A(u,v)\exp \left[i{\psi }_0(u,v)\right],$$

(2)

(3)

Generate an estimation ${p^{\prime }}_0(x,y)$ in the object domain by inverse Fourier transformation as

$${p^{\prime }}_0(x,y)=F^{-1}\{{P^{\prime }}_0(u,v)\}=F^{-1}\{A(u,v)\exp [i{\psi }_0(u,v)]\},$$

(3)

(4)

For the $(k+1)th$ iteration, i.e., the next iteration, the estimation of the object is determined as

$$p_{k+1}(x,y)=\{\begin{array}{l}{p^{\prime }}_k(x,y), (x,y)\notin S\\ p_k(x,y)-\beta {p^{\prime }}_k(x,y), (x,y)\in S\end{array},$$

(4)

where $\beta$ is a constant and $S$ are the regions at which ${p^{\prime }}_k(x,y)$ breaches the object-domain constraints, i.e., wherever ${p^{\prime }}_k(x,y)$ is negative or it goes beyond the known area of the object [20]. A brief procedure of the HIO method is shown in Figure 1a.

In this way, the image is reconstructed iteratively and its quality is judged by the mean-square errors (MSE) [11], i.e.,

$$MSE=\frac{1}{M\times N}{\displaystyle \sum _{x=1}^M{\displaystyle \sum _{y=1}^N{\left[p_{rec}(x,y)-p_{ori}(x,y)\right]}^2}}$$

(5)

where $p_{rec}(x,y)$ and $p_{ori}(x,y)$ mean the reconstructed and original images, respectively; $M$ and $N$ represent the two-dimensional size of the image.

The parameter $\beta$ in Equation (4) has a crucial impact on the results of the HIO algorithm. The impact could be embodied as the fact that the same initial distribution with different choices of $\beta$ leads to dramatically different reconstructions, and there is no guarantee that the best $\beta$ can be found to run the HIO algorithm due to the nonconvexity of phase retrieval. Therefore, most of the time, a better result is obtained by performing many runs of HIO with different initial distributions or different $\beta$ and averaging their results, which requires large computational time.

However, HIO does not have to run alone because it can be looked on as a suboptimal version of a more general approach [20]. This implies that the HIO method could be combined with other methods and be regarded as a part of a more universal method. Therefore, we expect to propose a more general approach, in which the rule of utilizing fixed $\beta$ in a single run of HIO is changed by introducing other algorithms. In this way, a better result could be brought about with fewer runs and less computational time. In other words, a better reconstruction could be realized by the modified method in a single run.

3. Simulated Annealing Hybrid Input-Output

Inspired by some existing algorithms that apply simulated annealing for phase retrieval [21,22,23], in this work we use the concept of simulated annealing in HIO to address the problem of fixed $\beta$ during the iterations.

As an effective probabilistic technique, simulated annealing plays an important role in approximating the global optimum of a given function. It derives from the physical process of annealing in metallurgy, to make perfect crystals by heating and controlled cooling. In this way, particles in crystals can be arranged in a state with minimal possible energy, provided that the beginning temperature is sufficiently high and the cooling is sufficiently slow [24].

There are usually four steps in simulated annealing.

(1)

Start with an appropriate annealing temperature $T$;

(2)

Give a perturbation to get a random state close to the current one and calculate the cost function $F$ of the new one corresponding to the energy physically;

(3)

Accept the new state if $\Delta F<0$ where $\Delta F=F^{new}-F^{old}$, and otherwise accept it on a probabilistic basis under the Metropolis criterion [25]. The annealing probability is given by

$$p=\exp (-\Delta F/k_{B}T),$$

(6)

where $k_B$ is the Boltzmann constant;

(4)

Decrease the temperature and start the next iteration from (2).

Equation (6) demonstrates that the probability $p$ increases with $T$ and decreases with $\Delta F$, which means that particles still tend to accept the new state if it is a bit worse than before in a high temperature.

The proposed simulated annealing hybrid input-output algorithm hereafter referred to as SAHIO, is divided into five steps:

(1)

Generate a random initial distribution $p_0(x,y)$ in the object domain and give an initial temperature $T$ and a value of $\beta$. In our experiment, the initial temperature $T$ could be set as

$$T=\left|\right|A(u,v)|{|}^2$$

(7)

(2)

Get the Fourier transformation $P_0(u,v)$ just like Equation (1). Then substitute the modulus of $P_0(u,v)$, namely $|P_0(u,v)|$, with the measured data experimentally referred to as $A(u,v)$ and keep its original phase angles to get integrated ${P^{\prime }}_0(u,v)$. Next, perform inverse Fourier transformation to acquire an estimate ${p^{\prime }}_0(x,y)$ in the object domain, which is described in Equation (3);

(3)

For the $(k+1)th$ iteration, i.e., the next iteration, calculate $p_k^{″}(x,y)$ according to Equation (4) and then evaluate the errors of $p_k(x,y)$ and $p_k^{″}(x,y)$ respectively using an error function, i.e.,

$$ERF(p)={\displaystyle \sum _{u,v}\left|\right|A(u,v)-P(u,v){\Vert }^2}$$

(8)

Furthermore, $\Delta ER{F}_k$, regarded as an analogy with $\Delta F$, is obtained by

$$\Delta ER{F}_k=ERF\left(p_k^{″}\right)-ERF\left(p_k\right)$$

(9)

(4)

Accept $p_k^{″}(x,y)$ as $p_{k+1}(x,y)$ at the beginning of the next iteration if $\Delta ER{F}_k<0$ or $\Delta ER{F}_k⩾0$ but it satisfies

$$p=\exp \left[-\frac{\Delta ER{F}_k}{T}\right]\ge rand(1)$$

(10)

where $rand(1)$ is a random number generated by the computer, which satisfies uniform distribution between 0 and 1. In this way, we could simulate the annealing event according to probability $p$ because of $Pr\left[rand\left(1\right)\le p\right]=p$. In other words, $p\ge rand\left(1\right)$ means accepting $p_k^{″}(x,y)$ as $p_{k+1}(x,y)$ according to the simulated annealing probability p. Otherwise, reject $p_k^{″}(x,y)$ and make $p_k(x,y)$ as $p_{k+1}(x,y)$ to start the next iteration;

(5)

Change the value of β randomly between 0 and 1 for the next iteration. Moreover, decrease the temperature according to some criteria.

The whole process of the SAHIO method is displayed in Figure 1b. As for the criteria of cooling, there are two points to be noted:

(6)

Lowering temperature too often and fast is not advised sometimes; that is to say, it is better to drop $T$ during every $N$ iteration, where $N$ is often greater than 1 and it could be set empirically;

(7)

The strategy of reducing $T$ is various and in this work, it decreases exponentially, i.e.,

$$T_{k+1} = \gamma T_k$$

(11)

where γ is a constant between 0 and 1 and in this work it is set as 0.8.

To show the operability of the algorithm, a piece of pseudo-code of SAHIO is given in Algorithm 1: The pseudo-code of SAHIO.

Algorithm 1: SAHIO (Simulated Annealing Hybrid Input-Output)

1:  Generate a random initial distribution $p_0$ and give the constraint regions $S$;

2:  Set the initial temperature $T=\left|\right|A|{|}^2$ and set $\beta =0.5$;

3:  Give the number of internal and external iterations, e.g., $Gen$ and $Iter$;

4:  Caculate $ERF(p_0)$ according to Equation (8);

5:  for $gen=0,1,2,\cdots Gen$ do

6:   for $iter=0,1,2,\cdots Iter$ do

7:    set $k=Gen\times Iter+iter$;

8:    set $P_k=F\left\{p_k\right\}=\left|P_k\right|\exp [i{\psi }_k]$, ${P^{\prime }}_k(u,v)=A\exp \left[i{\psi }_k\right]$ and ${p^{\prime }}_k=F^{-1}\left\{{P^{\prime }}_k\right\}=F^{-1}\{A\exp [i{\psi }_k]\}$;

9:    caculate $p_{k+1}$ according to Equation (4);

10:    caculate $ERF(p_{k+1})$ according to Equation (8) and $\Delta ER{F}_k=ERF\left(p_{k+1}\right)-ERF\left(p_k\right)$;

11:    if $\Delta ER{F}_k<0$ then

12:     set $p_{k+1}=p_{k+1}$;

13:    else then

14:     generate a random number $rand(1)$;

15:     if $p=\exp \left[-\frac{\Delta ER{F}_k}{T}\right]\ge rand(1)$ then

16:      $p_{k+1}=p_{k+1}$;

17:     else then

18:      set $p_{k+1}=p_k$ and $ERF(p_{k+1})=ERF(p_k)$;

19:     end if

20:    end if

21:    change the value of $\beta$ randomly;

22:   end for

23:   set $T=\gamma T$;

24:  end for

In the above algorithm proposed, temperature $T$ and simulated annealing probability $p$ are both crucial parameters. Due to the fact of the different magnitude of $\Delta ER{F}_k$ for different objects, an appropriate initial temperature $T=\left|\right|A(u,v)|{|}^2$ is chosen at the beginning of the algorithm to constrain the ratio of $\Delta ER{F}_k$ to $T$ in Equation (10) to guarantee the probability $p$ of not more than 1 during iterations. Changing $\beta$ randomly between 0 and 1 could give a different step or direction for the HIO method to improve the performance while fixed beta is easier to lead to local minima.

Although it is an effective indicator for $\Delta ER{F}_k$ to judge the quality of the HIO’s current iterative result, greed and myopia are hidden in it meanwhile. A “temporarily worse” intermediate result $\Delta ER{F}_k\ge 0$ could lead to a “finally better” result. Therefore, we decrease the temperature according to some criteria such as Equation (11) to increase the probability of accepting the “temporarily worse” results in the later stage of iterations, which helps jump out the local minima and reach the global minimum.

4. Results

To compare the performance of different algorithms, a series of numerical simulations were implemented as follows.

First, a comparison between HIO and SAHIO methods was given. We performed the phase retrieval simulation with the same two original images and an identical initial distribution for different methods to guarantee uniqueness. It is worth mentioning that there exists a best initial temperature for the SAHIO, but how to find it is not in the scope of this work and we used Equation (7) to generate our initial temperature. Figure 2 shows the two original pictures; Figure 3 and Figure 4 demonstrate the reconstructed images and their MSEs after 200 iterations by HIO and SAHIO methods, respectively; Figure 5 displays the average MSEs during iterations.

From Figure 3 and Figure 4, it is obvious that the images reconstructed by the SAHIO method, with smaller MSEs, were clearer than the others. On average, the MSE yielded by the SAHIO methods was 64.42% lower than that of the HIO method for 2 images.

However, it is hard to reflect the advantages of the SAHIO algorithm, generally due to the small sample size (only two original images and one initial distribution). Therefore, to verify the universality of the SAHIO method, a large numerical simulation was carried out with 35 general images, in which 10 generated initial distributions are used to start up 2 algorithms, respectively (i.e., 10 runs are performed parallelly for each image of 35, and their results are averaged at last). The average MSEs during iterations are shown in Figure 5 as well.

Some observations and deductions could be made from Figure 5:

(1)

Overall, the MSEs yielded by two algorithms decreased with increasing numbers of iterations with different speeds and results, which means that both methods have dissimilar degrees of the ability to reconstruct images.

(2)

Generally speaking, the SAHIO method declined faster than HIO. However, it is worth noting that in the early iterations, SAHIO, at a slight disadvantage for decreasing the MSE, showed a bit greater fluctuation owing to randomly changing. It was the random process of the exploration that delayed the time in the former steps but gave more chances to avoid the local minima in the latter stage.

(3)

At the end of 200 iterations, for 35 images the average MSE yielded by the SAHIO method converged at the lowest value, i.e., 0.00745, while the HIO method ended up with 0.01493; that is to say, the MSE generated by the SAHIO method was 50.12% smaller than the HIO method, demonstrating its strength of jumping out of the local minima and searching for the global optimum.

We concluded that the SAHIO algorithm is better than the HIO algorithm both on speed and degree of convergence. Furthermore, theoretically, all algorithms based on HIO but replacing HIO with SAHIO are expected to get better results. Therefore, a second simulation was carried out to investigate the effect of the method combining the SAHIO and GHIO algorithms, termed as SAGHIO. The same 35 images with 6 initial distributions were used in this experiment, and GHIO was regarded as a baseline to be compared with SAGHIO. The average MSEs during 200 iterations are shown in Figure 6, and some inferences can be drawn:

(1)

The MSEs curve of SAGHIO converged faster than that of GHIO. In the beginning, SAGHIO performed a little worse; SAGHIO caught up quickly, but the speed advantage seemed less obvious with iterations going on;

(2)

At the end of 200 iterations, the MSEs of SAGHIO converged to a lower number, i.e., 0.00152, while GHIO was 0.00266. Due to the restriction of the axis scale, the difference of final MSEs between SAGHIO and GHIO was inapparent, but SAGHIO’s result was 42.857% better than GHIO’s indeed, which means SAGHIO could reconstruct images clearer in general.

Moreover, the physical experiment for testing the effect of the SAHIO method was designed. A simple imaging system shown in Figure 7 was set up. The light emitted by the laser ($\lambda =532 \mathrm{nm}$) passed through the beam expander and got collimated into parallel light by the convex lens group ($f_1=40 \mathrm{mm}$,$f_2=50.8 \mathrm{mm}$) and then illuminated the object, i.e., a $1 \mathrm{mm}$ triangle aperture. The Fourier magnitude of the triangle aperture was recorded on the focal plane of the third convex lens ($f_3=100 \mathrm{mm}$) by a CCD (Daheng Optics, Beijing, China, MER-030-120UM). Phase retrieval was performed by HIO ($\beta =0.9$) and SAHIO respectively in three runs, and the reconstructed objects, Fourier amplitudes and phases distributions are shown in Figure 8. Figure 9 shows their ERFs during 200 iterations, which are used to judge the quality of reconstructed objects.

Similar observations and deductions could be made from Figure 9. The ERFs of iterative results from SAHIO could get a lower level than those from HIO, which means that the quality of reconstructed images by SAHIO is better than that by HIO. Therefore, it can be concluded that SAHIO is more capable of reaching the global minima than HIO.

5. Discussion

Although the SAHIO method could address the plight of getting trapped in local minima to some extent, there still exists some room for improvement.

First, methods of generating and updating temperature could be explored much more. We introduced temperature $T$ to regulate the process of annealing in the hybrid modified algorithm, and the choice of initial temperature and the criteria of decreasing temperature are both of great importance. An appropriate initial temperature could not only ensure the probability $p$ is not more than 1 during iterations but also stimulate the potential of searching for the global minimum for the SAHIO algorithm. Proper criteria of cooling could control the process of the whole algorithm to speed it up or deepen it. In the experiments of this article, we used Equations (7) and (11) to generate and update the temperature, but there may be other better methods to do the same job. Therefore, it deserves more effort to explore better methods for all kinds of applications.

Second, not only MSE and ERF but also other ways to judge the quality of reconstructed images should be considered. The features of a reconstructed image could be reflected from many aspects and multiple methods could be integrated to evaluate the quality of the phase retrieval results comprehensively.

Last, changing $\beta$ randomly in the SAHIO algorithm of this article was a bit inefficient sometimes; a more appropriate and directional change is expected to replace the random one. All in all, more efforts on phase retrieval should be made for improvement in the future.

6. Conclusions

In this work, we incorporated simulated annealing into the HIO method to propose a combined phase retrieval scheme, namely the SAHIO algorithm. By the benefit of the SA method, the SAHIO algorithm lessens the possibility of locking into the local minima compared with the HIO method. A series of ideal universal pictures were used as the original images for different algorithms to compare their results and test the feasibility and superiority of the SAHIO scheme. It was concluded in the numerical simulations that the SAHIO method has more potential to deal with the nonconvex optimization, with the MSE 50.12% smaller than the HIO method when 35 original images were used in 10 runs. Although the SAHIO scheme lagged in speed at the beginning, it surpassed HIO gradually. Furthermore, SAHIO was applied in the interaction process like SAGHIO and it still led to better results compared with basic GHIO. Moreover, a physical experiment of retrieving phases of a triangle aperture was carried out to compare SAHIO with HIO, and the similar superiority of SAHIO was exhibited.

Both numerical and physical experiments indicate that SAHIO can be utilized in the Fourier spectrum based on image recovery techniques to reconstruct images better. The problem of introducing additional parameters was analyzed qualitatively, and a direction of improvement was given meanwhile. More efforts on phase retrieval algorithms need to be made in the future.