Mathematical Problems in Engineering
Volume 2016 (2016), Article ID 8257612, 7 pages
http://dx.doi.org/10.1155/2016/8257612

Research Article

Binary Sparse Phase Retrieval via Simulated Annealing

Wei Peng and Hongxia Wang

College of Science, National University of Defense Technology, Changsha, Hunan 410073, China

Received 11 March 2016; Revised 23 April 2016; Accepted 4 May 2016

Academic Editor: Pasquale Memmolo

Copyright © 2016 Wei Peng and Hongxia Wang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

This paper presents the Simulated Annealing Sparse PhAse Recovery (SASPAR) algorithm for reconstructing sparse binary signals from their phaseless magnitudes of the Fourier transform. The greedy strategy version is also proposed for a comparison, which is a parameter-free algorithm. Sufficient numeric simulations indicate that our method is quite effective and suggest the binary model is robust. The SASPAR algorithm seems competitive to the existing methods for its efficiency and high recovery rate even with fewer Fourier measurements.

1. Introduction

Phase retrieval, well known as the problem of recovering a signal from magnitudes of its Fourier transform, has numerous applications including crystallography [1], optical imaging [2], speech processing [3], and astronomical imaging [4].

Reconstruction of a signal from its magnitudes of a linear transform is a difficult task, which is usually formulated as a nonconvex optimization problem. Algorithms for the problem could be divided into the projection type and the lifting type. Projection algorithms like error reduction [5], Fienup [6], and HPR [7] project the signal onto two feasible sets alternately, at least one of which is nonconvex. These algorithms have been reported to be variants of their convex counterparts essentially [8]. Though they might work well in engineering applications, they usually lack convergence to the global minima. Phase lifting technique based on matrix completion [9] is developed to transform phase retrieval problems to be convex. However, the technique would square the number of variables for executing semidefinite programming (SDP), which restricts the lifting technique to small scale problems. More algorithms are proposed in recent; see Wirtinger Flow [10] and alternating minimization [11].

With the popularity of compressing sensing, the sparsity of signals also attracts attentions in phase retrieval. For the sparse cases, algorithms on different strategies are presented; see Bayesian Network based message passing [12], greedy strategy based GESPAR [13], projection based sparse-Fienup [14], lifting based SDP algorithms [15], and more. Analysis of uniqueness for sparse phase retrieval could be found [16, 17].

In the field of compressing sensing, sparse binary signals are also studied [18]. In this paper, we are interested in recovering a binary or boolean vector from magnitudes of its complex linear measurements without any phase information. That is, given measurement vectors , suppose that

then the goal is to reconstruct via and the measurement matrix

The setting of is various; here are three common types of measurements:

(1)The most classic type of is the discrete Fourier matrix. It is usually impossible to recover signal uniquely regardless of trivial transforms (shifting, mirroring, and phase shifting) without extra information provided. Fortunately, constrained by sparsity, the original signal might be reconstructed uniquely.(2)The measurement matrix could be complex Gaussian matrix inheriting compressive sensing.(3)The measurement matrix could be several masked discrete Fourier matrices, which increases the number of measurements.

The binary signal or boolean signal exists widely, such as the computer system and boolean objects in microscopic imaging systems. However, insufficient studies focus on the binary phase retrieval. In general, binary phase retrieval belongs to pseudo-boolean optimization, which is written as

where maps to .

Simulated annealing (SA) [19], as a probabilistic method to approximate the global optimum of the objective function, is often considered when the feasible set is discrete. Moreover, simulated annealing has been employed in the previous classic phase retrieval studies. Reference [20] performs simulated annealing on the level-quantized phase variables for front-wave imaging while [21] also applies SA to the phase variables. Reference [22] introduces SA for finding periodic phase relief structures via phase perturbations. In [23, 24], the SA methods are also applied to the phase retrieval problems but only require 0-order information without gradients. However, for the recently prevailing sparse phase retrieval, it seems that simulated annealing is still waiting to be exploited. Therefore, simulated annealing for Sparse PhAse Retrieval (SASPAR) is presented in this paper. We will identify the nonzeros in signals via simulated annealing based on the local gradient. The numerical experiments indicate the outstanding performance of SASPAR to recover signals from magnitudes of their Fourier transform. The recovery rate even outperforms GESPAR, which is the state-of-the-art algorithm for sparse phase retrieval. Especially in the cases with fewer measurements, our algorithm could also recover the signal with a high probability.

Notations and formulations for the binary sparse phase retrieval are given in Section 2. Algorithms are presented in Section 3. Several numerical simulations are conducted in Section 4 and Section 5 concludes the paper.

2. Notation and Setup

The matrices will be represented by upper letters. Vectors and vector-valued functions are denoted by lower bold letters like and . The th element of a vector will be denoted by . The conjugate transpose of a matrix will be denoted by . Pointwise multiplication will be represented by . The gradient of a real-valued function will be denoted by , the th element of which is denoted by . The square root of a vector is defined pointwisely; that is, . with a set in it denoting the total number of elements in . is also denoted by . The sparsity equals the number of nonzeros in .

A general version of sparse phase retrieval is

where is a complex measurement vector. The problem is essentially a quadratic optimization problem.

Estimation of Sparsity . For a unitary linear transform , including the Fourier transform, the sparsity is obvious, which could be estimated accurately by

For complex Gaussian matrix , the sparsity could be estimated by (see Appendix). In our work, is assumed to be given as the prior information.

Denote all the possible indexes where might be one as , which is also called the support of . In Fourier cases, zeros are added to the end of signals for obtaining high frequency Fourier measurements (see details in [13]), which is called oversampling.

In this paper, as we restrict the signal to be binary sparse and the sparsity is already known, then problem (4) is formulated as

where . It could also be formulated in the matrix form,

where and .

3. Algorithms

As the sparsity is assumed to be given, the problem is to identify the locations of ones in , which is a combinational optimization problem. Size of the feasible set is

The size would be quite huge. For example, let and . Then it will cost about evaluations to search the whole feasible set, which seems impossible for an ordinary computer at present. Simulated annealing is an effective tool for this type of problems. Algorithm 1 gives a general frame of simulated annealing. The details of the frame will be discussed.

Require: , , ,

    Initialization: , ,

    while do

        Select a random neighbour from according to probability the defined on

        if then

        else

            if then

            else

            end if

        end if

    end while

    return

Algorithm 1: Simulated annealing.

Cooling Schedule. An exponential cooling schedule is employed [25]. That is, the temperature in the -th iteration is defined as

where .

Neighbourhood. The neighbourhood is defined to construct the search graph. Suppose we already know and define a set map . First of all, is computed (see computation details in Appendix B). In particular, the gradient for the discrete Fourier transform matrix is

where denotes pointwise multiplication. Let the index set be

consists of the locations where with positive gradient. Then define neighbours of as

Neighbours of are the ones which are sparsity and swap one pair between and .

Randomized Selection. In order to avoid trapping by a local minima in a single loop, we select a neighbour randomly. The discrete distribution defined on is presented in two stages.

Stage 1. First, is selected according to the discrete probability density distribution on . We define the distribution function

The index with larger positive gradient tends to be selected.

Stage 2. After is determined, is set to and a sparsity vector is obtained. Define the index set

Since , , setting the th element of to could derive a sparsity vector . The greedy strategy is employed in this step. The index with smallest negative gradient is selected; that is,

Combined with the specification of the search graph and randomized selection, the Simulated Annealing Sparse PhAse Retrieval (SASPAR) for binary signals is presented as Algorithm 2 displays.

Require: , , ,

    Initialization: , ,

    while do

        Obtain

        Select with the probability:

        Select such that

        if then

        else

            if then

            else

            end if

        end if

    end while

    return

Algorithm 2: SASPAR.

A more naive and straight algorithm of searching the minima is also presented. Without probabilistic selection of the swapped indexes, simply greedy 2-opt [26] strategy is employed in Algorithm 3, that is, exchanging the index with positive maximal gradient in and the index with negative minimal gradient in . Greedy 2-opt binary phase retrieval is an algorithm which is easy to implement and parameter-free. The comparison is conducted in the simulation part.

Require: ,

    Initialization:

    for to do

        Select with positive maximal gradient.

        Select such that

        if then

        else

            Terminate the algorithm or restart.

        end if

    end for

    return

Algorithm 3: Greedy 2-opt.

4. Numeric Simulation

In order to illustrate the performance of SASPAR, several experiments are conducted. To our knowledge, there are few algorithms specified for binary phase retrieval. Since we lack specified existing algorithms, the state-of-the-art algorithm for sparse phase retrieval GESPAR is employed. SASPAR is compared with GESPAR by embedding binary signals to the field of real numbers; thus combinatorial optimization is converted to continuous optimization.

Fourier measurements are employed in the section. The initial temperature of SASPAR is usually set to be half of the signal length. The stopping temperature and cooling parameter are set to and constantly. Time for iterations of all algorithms in this section is limited to seconds by default. After the cooling procedure is terminated, the SASPAR will restart until it exceeds the time limit.

The evaluation includes its comparison to Algorithm 3 (greedy 2-opt), the sensitivity to max time allowed, recovery from fewer measurements (compared with GESPAR), and robustness to noise (compared with GESPAR). The following experiments run on desktops with i7 4790 CPU and 4 GB RAM. All the scripts implemented by MATLAB for the experiments are available online https://github.com/dynames0098/binary-sparse-phase-retrieval-SA.git.

4.1. Deterministic versus Heuristic

The deterministic greedy algorithm (greedy 2-opt) and the heuristic simulated annealing algorithm (SASPAR) are compared. The max time allowed is fixed to 20 seconds for both algorithms.

Figure 1 illustrates that the greedy 2-opt algorithm is still a competitive algorithm to SASPAR. Their performances are quite close to each other. The success rate is high with sparsity from 24 to 30. However, SASPAR outperforms the 2-opt algorithm from sparsity , which might be explained by the fact that SASPAR requires less time to obtain the minimum while 2-opt algorithm takes more.

Figure 1: SASPAR versus greedy 2-opt.

4.2. Max Time Allowed for Computation

As is known, more time allowed for searching the optimum results in a higher success rate. In this experiment, signals are recovered from magnitudes of Fourier measurements. The sparsity is fixed to and we double the amount of sampling by adding zeros to signals, which also double the length of signals. The size of signals is set to , , , and separately. The max time allowed varies from seconds to seconds. For each length and max time allowed, the computation is repeated times. The result is displayed in Figure 2.

Figure 2: Max time allowed versus recovery rate.

Although the recovery curve in Figure 2 is not smooth, it is clear that a higher recovery rate is obtained as the time allowed increases. Sparser signals tend to be recovered more successfully. seconds is enough to guarantee a high success rate of recovering a -length signal with nonzeros while a -length signal with nonzeros might cost more time to guarantee an acceptable success rate.

4.3. Fewer Fourier Measurements

The high frequency part of Fourier transform might not be available due to the limit of devices or other restrictions. It becomes more difficult to reconstruct signals from fewer Fourier measurements. A comparison between SASPAR and GESPAR is conducted. The length of support of signals is fixed to 512 while the measurements vary from 600 to 1000 (by concatenate corresponding zeros to the support of original signals). The sparsity is set to 30, 35, and 40 separately. The experiment is repeated 100 times. Figure 3 illustrates that SASPAR outperforms GESPAR in the given conditions. SASPAR could also maintain a higher recovery rate with fewer measurements.

Figure 3: Recover from fewer Fourier measurements: SASPAR versus GESPAR.

4.4. Recover from Noisy Measurements

SASPAR and GESPAR are compared to recover noisy signals with 4 db to 14 db SNR. The length of signals is fixed to 256 and 256 zeros are concatenated for Fourier oversampling. The experiment is repeated times. The Gaussian white noise is added to measurement . Figure 4 displays an illustration of a noisy signal with 4 db SNR.

Figure 4: An illustration of noisy measurements.

To our surprise, the binary phase retrieval is quite robust to noises. As Figure 5 shows, we could recover signals; even the SNR is as low as 4 db, where the signal seems totally overwhelmed by the noise. With proper algorithms selected, sparse binary signals are able to be recovered from seriously noisy signals. However, we still lack the theoretical explanations for robustness of binary sparse phase retrieval.

Figure 5: SNR versus recovery rate.

5. Conclusion

In this paper, we introduce the binary sparse phase retrieval model and present the efficient heuristic algorithm SASPAR. We specified the settings of the simulated annealing procedure and conducted a series of numerical experiments in various cases to testify its efficiency, robustness, and scalability.

Though SASPAR is employed to solve sparse binary phase retrieval problems, it should not be restricted to the binary case, which is quite special in phase retrieval. It might provide an efficient support selection strategy for all sparse recovery problems. As is accepted widely, GESPAR works significantly well in small scale problems. However, we notice that GESPAR requires more measurements for recovery than SASPAR. It implies that the greedy 2-opt support selection strategy might hinder its performance as experiment 1 suggests. Therefore, SASPAR will be adapted to recover the sparse real or complex cases in the future work. In addition, the sparse binary model is also attractive for its simplicity and robustness. As little theoretical evidence has been provided, our future work will also focus on the analysis of binary phase retrieval.

Appendix

A. Estimation of the Sparsity

Suppose that complex matrix and , where . are composed of i.i.d. standard real Gaussian variables . Let . Denote . Then

where random variable , since

B. Calculation of the Gradient x

Define the objective function

Suppose that , where and are real vectors, . Then the gradient is computed via the chain rule,

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

The authors are grateful for the support from National Natural Science Foundation of China (nos. 61571008 and 61201327).

References

R. W. Harrison, “Phase problem in crystallography,” Journal of the Optical Society of America A, vol. 10, no. 5, pp. 1046–1055, 1993.

A. Walther, “The question of phase retrieval in optics,” Journal of Modern Optics, vol. 10, no. 1, pp. 41–49, 1963.

R. Lawrence and B.-H. Juang, Fundamentals of Speech Recognition, 1993.

C. Fienup and J. Dainty, “Phase retrieval and image reconstruction for astronomy,” in Image Recovery: Theory and Application, pp. 231–275, Academic Press, New York, NY, USA, 1987.

R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of the phase from image and diffraction plane pictures,” Optik, vol. 35, no. 2, pp. 237–250, 1972.

J. R. Fienup, “Phase retrieval algorithms: a comparison,” Applied Optics, vol. 21, no. 15, pp. 2758–2769, 1982.

H. H. Bauschke, P. L. Combettes, and D. R. Luke, “Hybrid projection-reflection method for phase retrieval,” Journal of the Optical Society of America A: Optics and Image Science, and Vision, vol. 20, no. 6, pp. 1025–1034, 2003.

H. H. Bauschke, P. L. Combettes, and D. R. Luke, “Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization,” Journal of the Optical Society of America A, vol. 19, no. 7, pp. 1334–1345, 2002.

E. J. Candes, Y. Eldar, T. Strohmer, and V. Voroninski, “Phase retrieval via matrix completion,” SIAM Review, vol. 57, no. 2, pp. 225–251, 2015.

E. J. Candès, X. Li, and M. Soltanolkotabi, “Phase retrieval via wirtinger flow: theory and algorithms,” IEEE Transactions on Information Theory, vol. 61, no. 4, pp. 1985–2007, 2015.

P. Netrapalli, P. Jain, and S. Sanghavi, “Phase retrieval using alternating minimization,” in Proceedings of the 27th Annual Conference on Neural Information Processing Systems (NIPS '13), December 2013.

P. Schniter and S. Rangan, “Compressive phase retrieval via generalized approximate message passing,” IEEE Transactions on Signal Processing, vol. 63, no. 4, pp. 1043–1055, 2015.

Y. Shechtman, A. Beck, and Y. C. Eldar, “GESPAR: efficient phase retrieval of sparse signals,” IEEE Transactions on Signal Processing, vol. 62, no. 4, pp. 928–938, 2014.

S. Mukherjee and C. S. Seelamantula, “An iterative algorithm for phase retrieval with sparsity constraints: application to frequency domain optical coherence tomography,” in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP '12), pp. 553–556, Kyoto, Japan, March 2012.

Y. Shechtman, Y. C. Eldar, A. Szameit, and M. Segev, “Sparsity based sub-wavelength imaging with partially incoherent light via quadratic compressed sensing,” Optics Express, vol. 19, no. 16, pp. 14807–14822, 2011.

H. Ohlsson and Y. C. Eldar, “On conditions for uniqueness in sparse phase retrieval,” in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP '14), pp. 1841–1845, IEEE, Florence, Italy, May 2014.

Y. C. Eldar, P. Sidorenko, D. G. Mixon, S. Barel, and O. Cohen, “Sparse phase retrieval from short-time fourier measurements,” IEEE Signal Processing Letters, vol. 22, no. 5, pp. 638–642, 2015.

U. Nakarmi and N. Rahnavard, “BCS: compressive sensing for binary sparse signals,” in Proceedings of the IEEE Military Communications Conference (MILCOM '12), pp. 1–5, Orlando, Fla, USA, November 2012.

D. Pham and D. Karaboga, Intelligent Optimisation Techniques: Genetic Algorithms, Tabu Search, Simulated Annealing and Neural Networks, Springer Science & Business Media, 2012.

N. Yoshikawa and T. Yatagai, “Phase optimization of a kinoform by simulated annealing,” Applied Optics, vol. 33, no. 5, pp. 863–868, 1994.

D. Morris, “Simulated annealing applied to the Misell algorithm for phase retrieval,” IEE Proceedings—Microwaves, Antennas and Propagation, vol. 143, no. 4, pp. 298–303, 1996.

J. Turunen, A. Vasara, and J. Westerholm, “Kinoform phase relief synthesis: a stochastic method,” Optical Engineering, vol. 28, no. 11, Article ID 281162, 1989.

M. Nieto-Vesperinas, F. J. Fuentes, and R. Navarro, “Performance of a simulated-annealing algorithm for phase retrieval,” JOSA A, vol. 5, no. 1, pp. 30–38, 1988.

Y. Lin, T. J. Kessler, and G. N. Lawrence, “Design of continuous surface-relief phase plates by surface-based simulated annealing to achieve control of focal-plane irradiance,” Optics Letters, vol. 21, no. 20, pp. 1703–1705, 1996.

Y. Nourani and B. Andresen, “A comparison of simulated annealing cooling strategies,” Journal of Physics A: Mathematical and General, vol. 31, no. 41, pp. 8373–8385, 1998.

C. H. Papadimitriou and K. Steiglitz, Combinatorial Optimization: Algorithms and Complexity, Courier Corporation, 1982.