Binary Sparse Phase Retrieval via Simulated Annealing

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.


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 x 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 x ∈ {0, 1}  from magnitudes of its complex linear measurements without any phase information.That is, given measurement vectors a  ∈ C  , suppose that then the goal is to reconstruct x via y 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 1 −  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.

Mathematical Problems in Engineering
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 min  (x) where  maps {−1, 1}  to R. 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.

Notation and Setup
The matrices will be represented by upper letters.Vectors and vector-valued functions are denoted by lower bold letters like x and p().The th element of a vector x 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 (x) will be denoted by ∇(x), the th element of which is denoted by ∇  (x).The square root of a vector y ∈ R  is defined pointwisely; that is, √ y = [ √  1 , √  2 , . . ., √   ]  .| ⋅ | with a set  in it denoting the total number of elements in .{0, 1} is also denoted by B. The sparsity ‖x‖ 0 equals the number of nonzeros in x.
A general version of sparse phase retrieval is where a  ∈ C  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  = ∑  =1   /2 (see Appendix).In our work,  = ‖‖ 0 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 x.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 x to be binary sparse and the sparsity  is already known, then problem ( 4) is formulated as where   ≥ 0 (∀0 ≤  ≤ ).It could also be formulated in the matrix form, where  = [a 1 , a 2 , . . ., a  ]  and y = [ 1 ,  2 , . . .,   ]  .

Algorithms
As the sparsity ‖x‖ 0 =  is assumed to be given, the problem is to identify the locations of  ones in x, which is a combinational optimization problem.Size of the feasible set is The size would be quite huge.For example, let  = 1024 and  = 30.Then it will cost about 10 90 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  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 iter-th iteration is defined as where 0 <  < 1.
Neighbourhood.The neighbourhood is defined to construct the search graph.Suppose we already know x iter and define a set map (x) = { |   = 1}.First of all, ∇(x iter ) 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 1 consists of the locations where   = 1 with positive gradient.Then define neighbours of x iter as Neighbours of x iter are the ones which are  sparsity and swap one pair between (x iter ) and   (x iter ).
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 N(x iter ) is presented in two stages.
Stage 1. First,  1 ∈  1 is selected according to the discrete probability density distribution on  1 .We define the distribution function The index with larger positive gradient tends to be selected.

Stage 2.
After  1 is determined,   1 is set to 0 and a  − 1 sparsity vector x mid is obtained.Define the index set Since  2 ⊂   (x iter ), ∀ 2 ∈  2 , setting the  2 th element of x mid to 1 could derive a  sparsity vector x new ∈ N(x iter ).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  1 and the index with negative minimal gradient in  2 .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.

Numeric Simulation
In order to illustrate the performance of SASPAR, several experiments are conducted.To our knowledge, there are few Require:  0 ,  end , , x 0 Initialization: x new ← x 0 ,  ←  0 , iter ← 0 while  ≥  end do iter ← iter + 1  ←  Obtain  1 = { ∈  | ∇  (x iter ) > 0 ∩ (x iter )} Select  1 with the probability: ∇   (x iter ) x new ( 1 ) ← 0 Select  2 such that  2 = arg min Terminate the algorithm or restart.end if end for return x opt = arg min 0≤iter≤max Iter (x iter ) Algorithm 3: Greedy 2-opt.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  end = 20 and  = 0.99 constantly.Time for iterations of all algorithms in this section is limited to 20 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.

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  = 32, which might be explained by the fact that SASPAR requires less time to obtain the minimum while 2-opt algorithm takes more.

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 30 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 256, 512, 1280, and 2048 separately.The max time allowed varies from 2 seconds to 20 seconds.For each length and max time allowed, the computation is repeated 100 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.20 seconds is enough to guarantee a high success     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.

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.

Figure 4 :
Figure 4: An illustration of noisy measurements.

4. 4 .
Recover from Noisy Measurements.SASPAR and GES-PAR 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 100 times.The Gaussian white noise is added to measurement y. Figure 4 displays an illustration of a noisy signal with 4 db SNR.
1 gives Require:  0 ,  end , 0 <  < 1, x 0 Initialization: x new ← x 0 ,  ←  0 , iter ← 0 while  ≥  end do iter ← iter + 1  ←  Select a random neighbour from N(x new ) according to probability the P defined on N(x new ) Δ ←  (x new ) −  (x iter ) opt = arg min 0≤iter≤max Iter (x iter ) opt = arg min 0≤iter≤max Iter (x iter ) Require: max Iter, x 0 Initialization: x new = x 0 for iter = 0 to max Iter do Select  1 with positive maximal gradient. 1 = arg max ∈(x iter ) rate of recovering a 2048-length signal with 30 nonzeros while a 256-length signal with 30 nonzeros might cost more time to guarantee an acceptable success rate.
corresponding zeros to the support of original signals).The sparsity is set to 30, 35, and 40 separately.The experiment is repeated 100 times.Figure3illustrates that SASPAR outperforms GESPAR in the given conditions.SASPAR could also maintain a higher recovery rate with fewer measurements.