# Binary Sparse Phase Retrieval via Simulated Annealing

Authors: Wei Peng; Hongxia Wang
Journal: Mathematical Problems in Engineering
 (2016)

---

## 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.


---

## Body

## 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 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
}
n from magnitudes of its complex linear measurements without any phase information. That is, given measurement vectors a
i
∈
C
n, suppose that(1)
y
i
=
a
i
,
x
2
,
f
o
r
i
=
1
,
…
,
m
;then the goal is to reconstruct x via y and the measurement matrix(2)
A
=
a
1
,
a
2
,
…
,
a
m
H
.The setting of A is various; here are three common types of measurements: (1)
The most classic type ofA is the discrete Fourier matrix. It is usually impossible to recover 1
-
D 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 matrixA could be complex Gaussian matrix inheriting compressive sensing.
(3)
The measurement matrixA 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(3)
m
i
n
f
x
s
.
t
.
x
i
∈
-
1,1
,where f maps {
-
1,1
}
n 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 Section2. 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 likex and p
(
i
). The ith element of a vector x will be denoted by x
i. The conjugate transpose of a matrix A will be denoted by A
H. Pointwise multiplication will be represented by ⊙. The gradient of a real-valued function F
(
x
) will be denoted by ∇
F
(
x
), the ith element of which is denoted by ∇
i
F
(
x
). The square root of a vector y
∈
R
m is defined pointwisely; that is, y
=
y
1
,
y
2
,
…
,
y
m
T. · with a set S in it denoting the total number of elements in S. {
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(4)
m
i
n
F
x
=
∑
i
=
1
m
a
i
H
x
2
-
y
i
2
s
.
t
.
x
0
≤
k
,
0
<
k
≪
n
x
∈
R
n
,where a
i
∈
C
n is a complex measurement vector. The problem is essentially a quadratic optimization problem.Estimation of Sparsity
k. For a unitary linear transform U, including the Fourier transform, the sparsity is obvious, which could be estimated accurately by(5)
k
=
x
2
2
=
U
x
2
2
=
y
2
2
=
∑
i
=
1
m
y
i
.For complex Gaussian matrix G, the sparsity could be estimated by k
¯
=
∑
i
=
1
m
‍
y
i
/
2
m (see Appendix). In our work, k
=
x
0 is assumed to be given as the prior information.Denote all the possible indexes wherex
i might be one as I, 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 signalx to be binary sparse and the sparsity k is already known, then problem (4) is formulated as(6)
m
i
n
F
x
=
∑
i
=
1
m
a
i
H
x
2
-
y
i
2
s
.
t
.
x
0
=
k
,
0
≤
k
≪
n
x
∈
B
n
,
i
∣
x
i
≠
0
⊂
I
,where y
i
≥
0
(
∀
0
≤
i
≤
m
). It could also be formulated in the matrix form,(7)
m
i
n
F
x
=
A
x
⊙
A
x
-
y
2
2
s
.
t
.
x
0
=
k
,
0
<
k
≪
n
x
∈
B
n
,
i
∣
x
i
≠
0
⊂
I
,where A
=
[
a
1
,
a
2
,
…
,
a
m
]
H and y
=
[
y
1
,
y
2
,
…
,
y
m
]
T.


## 3. Algorithms

As the sparsityx
0
=
k is assumed to be given, the problem is to identify the locations of k ones in x, which is a combinational optimization problem. Size of the feasible set is(8)
x
∈
B
n
∣
x
0
=
k
=
n
k
≈
O
n
k
,
n
≫
k
.The size would be quite huge. For example, letn
=
1024 and k
=
30. Then it will cost about 1
0
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 1 gives a general frame of simulated annealing. The details of the frame will be discussed.Algorithm 1:Simulated annealing.
Require: T
0, T
e
n
d, 0
<
α
<
1, x
0
Initialization: x
n
e
w
←
x
0, T
←
T
0, i
t
e
r
←
0
while
T
≥
T
e
n
d
do
i
t
e
r
⟵
i
t
e
r
+
1
T
⟵
α
T
Select a random neighbour fromN
x
n
e
w according to probability the P defined on N
x
n
e
w
Δ
⟵
F
x
n
e
w
-
F
x
i
t
e
r
if
Δ
≤
0
then
x
i
t
e
r
+
1
⟵
x
n
e
w
else
if
exp
⁡
-
Δ
/
T
>
rand
⁡
(
1
)
then
x
i
t
e
r
+
1
⟵
x
n
e
w
else
x
i
t
e
r
+
1
⟵
x
i
t
e
r
x
n
e
w
⟵
x
i
t
e
r
end if
end if
end while
return
x
o
p
t
=
arg
⁡
m
i
n
0
≤
i
t
e
r
≤
max
⁡
I
t
e
r
F
x
i
t
e
rCooling Schedule. An exponential cooling schedule is employed [25]. That is, the temperature T in the i
t
e
r-th iteration is defined as(9)
T
i
t
e
r
=
T
0
α
i
t
e
r
,where 0
<
α
<
1.Neighbourhood. The neighbourhood is defined to construct the search graph. Suppose we already know x
i
t
e
r and define a set map S
(
x
)
=
{
i
∣
x
i
=
1
}. First of all, ∇
F
(
x
i
t
e
r
) is computed (see computation details in Appendix B). In particular, the gradient for the discrete Fourier transform matrix A is(10)
∇
F
x
=
4
A
H
A
x
2
-
y
⊙
A
x
,where ⊙ denotes pointwise multiplication. Let the index set be(11)
I
1
=
i
∣
∇
i
F
x
i
t
e
r
>
0
∩
S
x
i
t
e
r
.
I
1 consists of the locations where x
i
=
1 with positive gradient. Then define neighbours of x
i
t
e
r as(12)
N
x
i
t
e
r
=
S
x
i
t
e
r
∖
i
1
∪
i
2
∣
i
1
∈
I
1
,
i
2
∈
S
C
x
i
t
e
r
.Neighbours of x
i
t
e
r are the ones which are k sparsity and swap one pair between S
(
x
i
t
e
r
) and S
C
(
x
i
t
e
r
).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
i
t
e
r
) is presented in two stages.Stage 1.
First,i
1
∈
I
1 is selected according to the discrete probability density distribution on I
1. We define the distribution function(13)
p
1
i
=
∇
i
F
x
i
t
e
r
∑
i
∈
I
1
∇
i
F
x
i
t
e
r
,
i
∈
I
1
.The index with larger positive gradient tends to be selected.Stage 2.
Afteri
1 is determined, x
i
1 is set to 0 and a k
-
1 sparsity vector x
m
i
d is obtained. Define the index set(14)
I
2
=
i
∈
I
∣
∇
i
F
x
m
i
d
<
0
∩
S
C
x
i
t
e
r
.Since I
2
⊂
S
C
(
x
i
t
e
r
), ∀
i
2
∈
I
2, setting the i
2th element of x
m
i
d to 1 could derive a k sparsity vector x
n
e
w
∈
N
(
x
i
t
e
r
). The greedy strategy is employed in this step. The index with smallest negative gradient is selected; that is,(15)
i
2
=
a
r
g
min
i
∈
I
2
⁡
∇
i
F
x
m
i
d
.Combined with the specification of the search graph and randomized selection, the Simulated Annealing Sparse PhAse Retrieval (SASPAR) for binary signals is presented as Algorithm2 displays.Algorithm 2:SASPAR.
Require: T
0, T
e
n
d, α, x
0
Initialization: x
n
e
w
←
x
0, T
←
T
0, i
t
e
r
←
0
while
T
≥
T
e
n
d
do
i
t
e
r
⟵
i
t
e
r
+
1
T
⟵
α
T
ObtainI
1
=
i
∈
I
∣
∇
i
F
x
i
t
e
r
>
0
∩
S
x
i
t
e
r
Selecti
1 with the probability:
P
i
1
=
i
=
∇
i
F
x
i
t
e
r
∑
i
∈
I
1
‍
∇
i
F
x
i
t
e
r
x
n
e
w
i
1
⟵
0
Selecti
2 such that
i
2
=
arg
⁡
min
i
2
∈
S
C
x
iter
⁡
∇
i
F
x
n
e
w
x
n
e
w
i
2
⟵
1
Δ
=
F
x
n
e
w
-
F
x
i
t
e
r
if
Δ
≤
0
then
x
i
t
e
r
+
1
⟵
x
n
e
w
else
if
exp
⁡
-
Δ
/
T
>
rand
⁡
(
1
)
then
x
i
t
e
r
+
1
⟵
x
n
e
w
else
x
i
t
e
r
+
1
⟵
x
i
t
e
r
x
n
e
w
⟵
x
i
t
e
r
end if
end if
end while
return
x
o
p
t
=
arg
⁡
m
i
n
0
≤
i
t
e
r
≤
max
⁡
I
t
e
r
F
x
i
t
e
rA 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 I
1 and the index with negative minimal gradient in I
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.Algorithm 3:Greedy 2-opt.
Require: max
⁡
I
t
e
r, x
0
Initialization: x
n
e
w
=
x
0
for
i
t
e
r
=
0
to
max
⁡
I
t
e
r
do
Selecti
1 with positive maximal gradient.
i
1
=
arg
⁡
max
i
∈
S
x
iter
⁡
∇
i
F
x
i
t
e
r
x
n
e
w
i
1
⟵
0
Selecti
2 such that
i
2
=
arg
⁡
min
i
∈
S
C
x
iter
⁡
∇
i
F
x
n
e
w
x
n
e
w
i
2
⟵
1
Δ
=
F
x
n
e
w
-
F
x
i
t
e
r
if
Δ
≤
0
then
x
i
t
e
r
+
1
=
x
n
e
w
else
Terminate the algorithm or restart.
end if
end for
return
x
o
p
t
=
arg
⁡
m
i
n
0
≤
i
t
e
r
≤
max
⁡
I
t
e
r
F
x
i
t
e
r


## 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 toT
e
n
d
=
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 Algorithm3 (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.Figure1 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 k
=
32, 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 to30 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.Figure 2
Max time allowed versus recovery rate.Although the recovery curve in Figure2 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 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.


### 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. Figure3 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 repeated100 times. The Gaussian white noise is added to measurement y. 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 Figure5 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.


## 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.Figure1 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 k
=
32, 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 to30 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.Figure 2
Max time allowed versus recovery rate.Although the recovery curve in Figure2 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 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.


## 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. Figure3 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 repeated100 times. The Gaussian white noise is added to measurement y. 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 Figure5 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.



---

*Source: 8257612-2016-05-31.xml*