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 thatthen the goal is to reconstruct via and the measurement matrixThe 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 aswhere 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 iswhere 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 byFor 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 aswhere . 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.
| ||||||||||||||||||||||||||||||||||||||
Cooling Schedule. An exponential cooling schedule is employed [25]. That is, the temperature in the -th iteration is defined aswhere .
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 iswhere denotes pointwise multiplication. Let the index set be consists of the locations where with positive gradient. Then define neighbours of asNeighbours 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 functionThe 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 setSince , , 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.
| ||||||||||||||||||||||||||||||||||||||||||||||||||
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.
| ||||||||||||||||||||||||||||||||||
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.
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.
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.
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.
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.
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 . Thenwhere 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).