1

Inexact Alternating Optimization for Phase
Retrieval In the Presence of Outliers

arXiv:1605.00973v2 [cs.IT] 24 Jul 2017

Cheng Qian, Xiao Fu, Member, IEEE, Nicholas D. Sidiropoulos, Fellow, IEEE, Lei Huang, Senior Member, IEEE,
and Junhao Xie, Senior Member, IEEE

Abstract—Phase retrieval has been mainly considered in the
presence of Gaussian noise. However, the performance of the
algorithms proposed under the Gaussian noise model severely
degrades when grossly corrupted data, i.e., outliers, exist. This
paper investigates techniques for phase retrieval in the presence
of heavy-tailed noise – which is considered a better model for
situations where outliers exist. An `p -norm (0 < p < 2) based
estimator is proposed for fending against such noise, and twoblock inexact alternating optimization is proposed as the algorithmic framework to tackle the resulting optimization problem.
Two specific algorithms are devised by exploring different local
approximations within this framework. Interestingly, the core
conditional minimization steps can be interpreted as iteratively
reweighted least squares and gradient descent. Convergence
properties of the algorithms are discussed, and the CramérRao bound (CRB) is derived. Simulations demonstrate that the
proposed algorithms approach the CRB and outperform stateof-the-art algorithms in heavy-tailed noise.
Index Terms—Phase retrieval, iterative reweighted least
squares (IRLS), gradient descent, impulsive noise, Cramér-Rao
bound (CRB).

I. I NTRODUCTION
Phase retrieval aims at recovering a signal x ∈ CN
from only the magnitude of linear measurements. This is
an old problem [2], [3] that has recently attracted renewed
and growing interest. Phase retrieval arises in many fields,
such as X-ray crystallography, coherent diffraction imaging,
and optical imaging and astronomy, where the detectors only
record the intensity information, because phase is very difficult
and expensive to measure.
The mathematical description of the phase retrieval problem
is simple: given the measuring matrix A = [a1 · · · aM ]H ∈
CM ×N and the measurement vector y ∈ RM , where am ∈
CN , and
y = |Ax|,
find x ∈ CN . Early attempts to solve the phase retrieval
problem can be traced back to the 1970s, where techniques
such as Gerchberg-Saxton (GS) [2], Fienup [3] and their variants were proposed. These algorithms have been empirically
Original manuscript submitted to IEEE Trans. on Signal Processing,
October 4, 2016; revised July 25, 2017. Conference version of part of this
work appears in Proc. EUSIPCO 2016 [1].
C. Qian, X. Fu and N. D. Sidiropoulos are with the Department of Electrical
and Computer Engineering, University of Minnesota, Minneapolis, MN 55455
USA (e-mail: alextoqc@gmail.com, xfu@umn.edu, nikos@umn.edu).
L. Huang is with the College of Information Engineering, Shenzhen
University, Shenzhen, 518060 China (e-mail: dr.lei.huang@ieee.org).
J. Xie is with the Department of Electronics and Information Engineering,
Harbin Institute of Technology, Harbin, 150001 China (e-mail: xj@hit.edu.cn).
The work of N. Sidiropoulos was supported by NSF CIF-1525194.

shown to work well under certain conditions, although little
had been known regarding their convergence properties from a
theoretical point of view. Recently, some new algorithms along
this line of work were proposed in [4], where the convergence
issue is better studied.
In recent years, more modern optimization-based approaches for phase retrieval have been proposed. For example,
Candès proposed a semidefinite relaxation approach known as
the PhaseLift algorithm [5] and proved that exact recovery
is possible with high probability in the noiseless case. Hand
recently studied the robustness of PhaseLift and showed that
it can tolerate a constant fraction of arbitrary errors [6].
Wirtinger-Flow (WF) [7] is a more recent approach that combines a good statistical initialization with a computationally
light gradient-type refinement algorithm. The combination
works very well when A is i.i.d. Gaussian. Following [7],
conceptually similar approaches, namely, the truncated WF
(TWF) [8] and truncated amplitude flow (TAF) algorithms [9],
[10] were proposed to handle more challenging scenarios. In
the recent work in [11], [12], the authors proposed a leastsquares feasible point pursuit (LS-FPP) approach that aims to
solve the same optimization problem as LS PhaseLift in the
presence of noise. Simulations in [11] indicate that LS-FPP
approaches the Cramér-Rao bound (CRB) for the Gaussian
measurement model considered in [11].
Some phase retrieval algorithms were originally developed
under an exact measurement model, and subsequently treated
the noisy case by replacing equality constraints with relaxed
inequality constraints [5]-[13]. Most of the existing algorithms
were explicitly or implicitly developed under a Gaussian noise
model. In certain applications, a subset of the measurements
may be corrupted much more significantly than the others, and
heavy-tailed noise may be encountered as well [14]-[16]. One
representative example is high energy coherent X-ray imaging
using a charge-coupled device (CCD), where the impulsive
noise originates from X-ray radiation effects on the CCD,
and the density of impulses increases with the intensity of
X-ray radiation or CCD exposure time [16]-[17]. Under such
circumstance, modeling the noise as Gaussian is no longer
appropriate.
In recent years, robust phase retrieval algorithms, e.g.,
[14], [15], [18], have been proposed to handle outliers. The
framework in [15] considered an undersampled phase retrieval
model corrupted with Laplacian-distributed outliers, but the
approach was designed specifically for sparse x. Variations of
TWF to handle outliers were also considered. For example, the
work in [18] replaces the sample mean that is used in TWF

2

by sample median in the initialization and truncated update,
which exhibits robustness to outliers under certain conditions
[18]. However, as TWF, this approach still heavily relies on
the assumption that Gaussian measurements are employed.
Another important aspect of phase retrieval is how noise
enters the measurement model. Our previous work [11] considered the noise model y = |Ax|2 + n. Another noise model
that is frequently considered in the literature [2]-[4] is
y = |Ax| + n.

(1)

Notice the subtle but important difference between the two
models: whether noise is added to the magnitude or the
squared magnitude. The choice hinges on the experimental
setup, including the measurement apparatus; but (1) is more
widely adopted by experimentalists.
Contributions: We consider the phase retrieval problem under
the model in (1) in the presence of impulsive noise, and
focus on designing robust algorithms to handle it. To fend
against impulsive noise, we adopt the `p -fitting (0 < p < 2)
based estimator that is known to be effective in dealing with
outliers, and devise two optimization algorithms using twoblock inexact alternating optimization. Specifically, the two
algorithms both solve local majorization-based approximate
subproblems for one block, instead of exactly solving the
conditional block minimization problem to optimality. Interestingly, starting from different majorizations, the resulting solutions turn out equivalent to iteratively reweighted least squares
and gradient descent, respectively. Unlike the existing inexact
and exact alternating optimization frameworks that mostly
operate with convex constraints, the proposed algorithms work
with a unit-modulus nonconvex constraint, so convergence
analysis seems difficult. Nevertheless, we prove convergence
of the proposed algorithms to a Karush-Kuhn-Tucker (KKT)
point by exploiting the two-block structure of the problem.
We also derive computationally light implementations using
Nesterov-type and stochastic gradient updates. In order to
assist in performance analysis and experimental design, we
derive the CRB for the model in (1) and under different
parameterizations, in the presence of Laplacian and Gaussian
noise. Curiously, although related CRBs for different noisy
measurement models have been previously derived in [11],
[19]-[22], to the best of our knowledge, there is no available
CRB for the model in (1) – and our work fills this gap. The
proposed algorithms are validated by extensive simulations.
The simulations show that our approaches outperform the
state-of-the-art algorithms in the presence of impulsive noise.
A preliminary conference version of part of this work
has been submitted to EUSIPCO 2016 [1]. The conference
version includes one of the two basic algorithms, the Laplacian
CRB, and limited simulations. The second algorithm, Nesterov
acceleration and stochastic gradient-type updates, additional
CRB results, and, most importantly, proof of convergence of
the iterative algorithms are all provided only in this journal
version, which naturally also includes more comprehensive
simulation results.
Notation: Throughout the paper, we use boldface lowercase
letters for vectors and boldface uppercase letters for matrices. Superscripts (·)T , (·)∗ , (·)H , (·)−1 and (·)† represent

transpose, complex conjugate, conjugate transpose, matrix
inverse and pseudo-inverse, respectively. The Re{·} and Im{·}
denote the real part and imaginary part. E[·] is the expectation
operator, | · | is the absolute value operator, || · ||F is the
Frobenius norm, || · ||p is the vector `p -norm, whose definition
PN
is ||x||p = ( i=1 |xi |p )1/p ,
is the element-wise product,
and diag(·) is a diagonal matrix with its argument on the
diagonal. δij denotes the Kronecker delta function, and ∠(·)
takes the phase of its argument. 0m , 1m , and Im stand for
the m × 1 all-zero vector, m × 1 all-one vector, and m × m
identity matrix, respectively. Furthermore, trace(·) and ∂a/∂x
denote the trace and partial derivative operators, respectively.
II. P ROPOSED A LGORITHMS
A. AltIRLS
Let us first consider the noiseless case where y = |Ax|.
Effectively, it can also be written as
y

u = Ax

(2)

where u = ej∠(Ax) is an auxiliary vector of unknown unitmodulus variables with its mth component being um =
H
ej∠(am x) . In the presence of impulsive noise, `p -(quasi)-norm
has be recognized as an effective tool for promoting sparsity
and fending against outliers [23]-[25]. Therefore, we propose
to employ an `p -fitting based estimator instead of using the
`2 -norm, which has the form of
min

|u|=1,x

M
X

2
|ym um − aH
m x| + 

p/2

(3)

m=1

where 0 < p < 2 is chosen to down-weigh noise impulses
(i.e., outliers), and  > 0 is a small regularization parameter
(e.g.,  ∈ [10−8 , 10−6 ]) that keeps the cost function within its
differentiable domain, which will prove handy in devising an
effective algorithm later.
To handle Problem (3), we follow the rationale of alternating
optimization, i.e., we first update x fixing u, and then we do
the same for u.
Assume that after some iterations, the current solution at
iteration r is (x(r) , u(r) ). At step (r + 1), the subproblem
with respect to (w.r.t.) x is
x(r+1) = arg min
x

M 
p/2
X
H 2
|ym u(r)
m − am x| + 

(4)

m=1

which is still difficult to handle. Particularly, when p < 1,
the subproblem itself is still non-convex; when p ≥ 1, the
subproblem is convex but has no closed-form solution. Under
such circumstances, we propose to employ the following
lemma [26]:
Lemma 2.1: Assume 0 < p < 2,  > 0, and φp (w) :=
p
  p−2
2−p 2
+ w. Then, we have
2
pw
x2 + 

p/2

= min wx2 + φp (w)
w≥0

and the unique minimizer is
 p−2
p 2
wopt =
x + 2 .
2

(5)

(6)

3

By
PM

Lemma
2.1,
an
upper
bound
of

p/2
(r)
H 2
that is tight at the current
m=1 |ym um − am x| + 

solution x(r) can be easily found:
M 
X

|ym u(r)
m

−

2
aH
m x|

+

Algorithm 1 AltIRLS for phase retrieval
1:
2:

(0)

W


p/2

p
2

m=1
M 
X
H
(r)
≤
ym u(r)
wm
m − am x

2


(r)
+ φp wm



(7)

m=1

3:
4:
5:

where
(r)
wm

p
:=
2



ym u(r)
m

−

(r)
aH
mx

6:

 p−2
2

2

+

(8)
7:
(r)

and the equality holds if and only if x = x . Instead
of directly dealing with Problem (4), we solve a surrogate
problem using the right hand side (RHS) of (7) at each iteration
to update x. Notice that the RHS of (7) is convex w.r.t. x and
the corresponding problem can be solved in closed-form:
†

x(r+1) = W(r) A W(r) y u(r)
(9)
where
W(r) = diag

q

(r)

w1

···


q
(r)
wM .

(10)

The conditional problem w.r.t. u is
u(r+1) = arg min

|u|=1

M 
X

(r+1) 2
|ym um − aH
| +
mx

p/2

.

m=1

(11)
Although the problem is non-convex, it can be easily solved to
optimality. Specifically, the first observation is that the partial
minimization w.r.t. u is insensitive to the value of p; i.e., given
a fixed x, for any p > 0, the solutions w.r.t. u are identical.
Second, for all p > 0, the solution is simply to align the angle
(r+1)
of ym um to that of aH
, which is exactly
mx
H

(r+1)

u(r+1)
= ej∠(am x
m

function A LT IRLS (y, A, x(0) )
(0)
Initialize u(0)q= exp(∠(Ax
)) and W(0)

q

),

m = 1, . . . , M.

(12)

We update x and u alternately, until some convergence criterion is met. We see that the way that we construct the upper
bound of the partial problem w.r.t. x is in fact the same as the
procedure in iteratively reweighted least squares (IRLS) [25],
[26]. The difference is that we ‘embed’ this IRLS step into
an alternating optimization algorithm. We therefore call this
algorithm alternating IRLS (AltIRLS), which is summarized
in Algorithm 1.
A relevant question regarding Algorithm 1 is whether or
not this algorithm converges to a meaningful point, e.g., a
stationary or KKT point. Note that the block variable u is
constrained to a non-convex set, and we do not compute the
optimal solution for the block variable x at each iteration
of Algorithm 1. For such a type of algorithm, there is no
existing theoretical framework that establishes convergence.
We therefore need careful custom convergence analysis. We
have the following result:
Proposition 2.1: Assume that 0 < p < 2 and  > 0, and that
A has full column rank. Then, the solution sequence produced

8:
9:
10:

=

diag

(0)

(0)
w1

(0)
ym um − aH
mx

···
 p−2
2
2
+

(0)
wM

=

(0)
wm

with

=

while stopping criterion has not been reached do
x(r) = (W(r−1) A)† W(r−1) (y u(r−1) ).
u(r) = exp(j∠(Ax(r) ))

 p−2
2
2
(r)
(r)
(r)
wm = p2 ym um − aH
, ∀m
x
+

m
q

q
(r)
(r)
W(r) = diag
w1 · · ·
wM
r =r+1
end while
end function

by Algorithm 1 converges to a set K that consists of all the
KKT points of Problem (3).
Proof: See Appendix A.
The result in Proposition 2.1 is interesting: Although the algorithm that we propose to compute the `p fitting-based estimator
solves non-convex subproblems, it ensures convergence to a
KKT point. Such a nice property is proven by exploiting the
two-block structure of the algorithm. Note that the proof itself
does not rely on the specific form of the optimization problem
in (3), and thus can be easily extended to other two-block
alternating optimization cases, which we believe is of much
broader interest.
B. AltGD
Algorithm 1 uses a simple update strategy, but its complexity may still become an issue when the problem size grows.
Specifically, the bottleneck of the AltIRLS algorithm lies in
solving the subproblem w.r.t. x, i.e.,
min
x

W(r) (y

u) − W(r) Ax

2

(13)

2

whose closed-form solution can be computed via (9). It is
observed from (9) that the matrix inversion part requires
O(N 3 ) flops to compute and also O(N 2 ) memory to store,
both of which are not efficient for high-dimensional x – e.g.,
when x is a vectorized image, N is usually very large (more
specifically, for a 100 × 100 image, N is 10,000).
To circumvent this difficulty, one way is to employ a certain
fast least squares solver such as conjugate gradient [28] to
handle Problem (13). This is a viable option, but such solvers
still typically require many iterations to obtain a fairly accurate
solution of (13). Here, we propose to deal with problem (13)
using a simpler approach. Specifically, let us denote
f (r) (x) = W(r) (y

u(r) ) − W(r) Ax

2
2

(14)

and approximate (13) using the following function
g (r) (x) = f (r) (x(r) ) + Re{(∇f (r) (x(r) ))H (x − x(r) )}

4

µ(r)
x − x(r)
2

2

.
2

3

where µ(r) ≥ 0 is a pre-specified parameter. Instead of optimizing f (r) (x), we optimize g(r) (x). By rearranging terms,
the subproblem becomes

 2
1
(16)
x(r+1) = arg min x − x(r) − (r) ∇f (x(r) )
x
µ
2
and the solution is
x(r+1) = x(r) −

1
∇f (x(r) )
µ(r)

(17)

if µ(r) ≥ λmax (AH (W(r) )2 A), and the equality holds if
and only if x = x(r) . Therefore, the algorithmic structure
of Algorithm 2 is the same as that of Algorithm 1, except
that the majorizing functions of the x-block are different –
which means that the proof of convergence of Algorithm 1
also applies here:
Corollary 2.1: If µ(r) ≥ λmax (AT (W(r) )2 A) for all r,
then the whole sequence produced by Algorithm 1 converges
to K.
Algorithm 2 AltGD for phase retrieval

5:
6:
7:
8:
9:

2

10

1

which does not require any matrix inversion operation. Also,
the N × N matrix does not
2 need to be stored, if2 we take
the order Ax(r) → W(r) Ax(r) → AH W(r) Ax(r) to
compute the update of x. Therefore, using this approach, both
memory and complexity requirements are less demanding.
We summarize the algorithm in Algorithm 2. We call this
algorithm alternating gradient descent (AltGD).
Note that when µ is chosen as the Lipschitz constant, (15)
is an upper bound of the cost function in (13). In other words,
we have
f (r) (x) ≤ g (r) (x)

2:
3:
4:

AltGD
Accelerated AltGD

10

i.e., a gradient step with step-size 1/µ(r) , where the gradient
is
2

∇f (r) (x(r) ) = AH W(r) Ax(r) − y u(r)
(18)

1:

10

(15)

|| y−|Ax| ||pp

+

function A LT GD(y, A, x(0) )
Initialize u(0) = exp(∠(Ax(0) ))
while stopping criterion has not been reached do
Choose µ(r−1) as the leading eigenvalue of
H
A (W(r−1) )2 A or the trace of (W(r−1) )2
x(r) = x(r−1) − 1/µ(r−1) ∇f (x(r−1) )
u(r) = exp(j∠(Ax(r) ))
r =r+1
end while
end function

Remark 2.1: Exactly computing λmax (AH (W(r) )2 A) may
be time consuming in practice. Many practically easier ways
can be employed, e.g., the Armijo rule-based line search
[27]. In our simulations, we use a simple heuristic: we
let µ(r) = trace((W(r) )2 ) instead of λmax (AH (W(r) )2 A).
The rationale behind is that we observe that the energy of AH (W(r) )2 A is usually dominated by W(r) ,

0

10

0

20

40

60

80
100
120
Number of iterations

140

160

180

200

Fig. 1. Cost function value versus number of iterations.

and using µ(r) = trace((W(r) )2 ) is a good approximation of trace(AH (W(r) )2 A) that is an upper bound of
λmax (AH (W(r) )2 A). We should mention that this step-size
choice is a heuristic, but works well in practice, as will be
shown in the simulations.
C. Further Reducing Complexity
1) Extrapolation: Algorithm 2 is easier to compute than
Algorithm 1 in terms of per-iteration complexity. However,
first-order methods tend to require more iterations in total.
One way to alleviate this effect is to incorporate Nesterov’s
“optimal first-order” method [29]-[31], i.e., for each update,
we set

t(r−1) − 1  (r)
(r−1)
z(r) = x(r) +
x
−
x
(19)
t(r)
p
1 + 1 + 4(t(r−1) )2
(20)
t(r) =
2
1
x(r+1) = z(r) − (r) ∇f (r) (z(r) )
(21)
µ
In practice, the above ‘extrapolation’ technique greatly expedites the whole process in various applications [23], [32].
Some numerical evidence can be seen in Fig. 1, where a
simple comparison between the plain Algorithm 2 and the
extrapolated version is presented. We choose SNR= 20 dB,
N = 16 and M = 128. Each element in the signal and measurement vectors is independently drawn from the complex
circularly symmetric Gaussian distribution with mean zero and
variance one. The noise is generated from a symmetric α stable
(SαS distribution) which will be described in detail in Section
IV. Fig. 1 shows the convergence of AltGD and accelerated
AltGD using extrapolation when p = 1.3. As we can see, the
accelerated AltGD converges after 40 iterations which is much
faster than AltGD that converges after 200 iterations.
2) Block Incremental / Stochastic Gradient: When N is
very large (thus M > N is larger), even gradient computation
is too much of a burden. In such cases, a pragmatic way
is to “break down” the problem to pieces and do (block)

5

incremental or stochastic gradient. We divide the measurement
matrix into L blocks Γl . Then it is straightforward to get the
gradient for the lth block as
(r) 2
(r) 
∇fΓl (x(r) ) = AH
AΓl x(r) − yΓl uΓl . (22)
Γl WΓl
The estimate of x is updated according to
x(r+1) = x(r) −

1
(r)

µΓl

∇fΓl (x(r) )

(23)

Proof: See Appendix B.
The CRB for real x is a special case of the complex case,
which can be easily derived from Proposition 3.1:
Proposition 3.2: In Laplacian noise, the CRB of real-valued
x is


CRBL,r = trace F−1
(27)
L,r
where
FL,r =

(r)

where µΓl is chosen as the leading eigenvalue of AH
Γl AΓl or
(r) 2 
trace WΓl
. The algorithm can proceed block by block
with revisits, or by randomly picking blocks, resulting in
block incremental gradient and stochastic gradient versions,
respectively.
A subtle point here is that, to maintain robustness, one
should choose a block size larger than one. The reason is
that the robustness of the algorithm is brought by treating
different ym with different weights (more specifically, by
downweighting the outliers). When using only one ym for
updating, this ability vanishes.
III. C RAM ÉR -R AO B OUND A NALYSIS
In this section, the CRB on the accuracy of retrieving x
in (1) is presented. The CRB provides a lower bound on the
MSE of unbiased estimators. In many signal processing applications, one cannot guarantee that an estimator is unbiased,
or even that an unbiased estimator exists, yet the CRB is still
surprisingly predictive of estimator performance [35], [38]. For
phase retrieval under Gaussian noise, over the past few years,
several CRBs have been derived for different models, e.g., 2-D
Fourier-based measurements [19], noise added prior to taking
the magnitude [20] and noise added after taking the magnitude
square [11], [21]. To the best of our knowledge, there is no
available CRB formula for the signal model in (1). Here, we
present the CRBs for two particular types of noise: Laplacian
and Gaussian noise. Although our main interest here lies in
evaluating performance of robust algorithms and the Laplacian
CRB serves this purpose, the Gaussian CRB is of interest in
other application contexts. Note that we use subscripts r, c,
L and G to stand for real, complex, Laplacian and Gaussian,
respectively.
To get started, we should note that our derivations are
based on the assumption that aH
m x is nonzero. That is because
the term |Ax| is non-differentiable at aH
m x = 0, the Fisher
information only exists at aH
x
=
6
0.
With
this caveat, we
m
have the following proposition:
Proposition 3.1: In Laplacian noise, the CRB is


CRBL,c = trace F†L,c
(24)

GL,r = Re{AH diag(Ax)}.

(25)

CRBG = 2CRBL .

GL,c



Re{AH diag(Ax)}
=
.
Im{AH diag(Ax)}

(26)

(30)

Proof: The proof is straightforward by calculating the
Gaussian FIM and comparing it with the Laplacian FIM in
(28). We omit it here.
In certain applications of phase retrieval, we may be more
interested in the performance bound for retrieving the phase
of x. The following proposition provides the CRB on both
the phase and amplitude of x under Laplacian noise. Note
that in the Gaussian noise case, it is straightforward to apply
Proposition 3.3 to compute the lower bound. Similar results
can also be found in [36].
Proposition 3.4: In Laplacian noise, the CRB of the amplitude of x is
CRBL,|x| =

N
X

di

(31)

i=1

and the variance of any biased estimate of the phase of x is
bounded below by
CRBL,∠(x) =

2N
X

di

(32)

i=N +1

where d = [d1 · · · d2N ] contains the main diagonal elements
of F†L which in this case is defined as
2
GL diag(|Ax|−2 ) GTL
σn2

(33)

with
GL =

with

(29)

In Appendix C, we show that when A has full column
rank, FL,c is singular with rank (2N −1) while FL,r is always
nonsingular. Thus, for complex x, we adopt its pseudo-inverse
to compute an optimistic (looser) CRB, which is still a valid
lower bound that can be used to benchmark the efficiency of
any biased estimator [33]-[35]. If x is close to zero, then the
CRB is not tight at low SNRs, and more measurements should
be used to approach the CRB.
Proposition 3.3: In the Gaussian noise case, the CRB is two
times larger than the CRB in Proposition 3.1, i.e.,

FL =
2
GL,c diag(|Ax|−2 ) GTL,c
σn2

(28)

with

where
FL,c =

2
GL,r diag(|Ax|−2 ) GTL,r
σn2


diag(|x|)−1


IN



Re diag (x∗ ) AH diag(Ax)
.
Im{diag(x∗ )AH diag(Ax)}
(34)

Proof: See [11] and Appendix B.

6

In this section, we evaluate the performance of the proposed
methods and compare them with some existing algorithms,
namely, WF [7], TWF [8], TAF [9], AltMinPhase [4], MTWF
[18] and GS [2] in terms of MSE performance, where the
MSE is computed after removing the global phase ambiguity
between the true and estimated x. In the simulations, we
consider an exponential signal x = exp(j0.16πt) comprising
16 samples, i.e., t = 1, · · · , 16. We test the algorithms
using different types of A, i.e., random matrix that is usually
employed in [4], [7]-[9], [18], and the 2D Fourier matrix that
is widely used in real devices [2], [3] (cf. subsection IV-E).
For the random measuring matrix case (used in subsections
A-D), the measurement vectors are generated from a masked
Fourier transformation, which takes the form of

T
A = (DΛ1 )T · · · (DΛK )T
where K = M/N is the number of masks, D is a N × N
discrete Fourier transform (DFT) matrix with DDH = N IN
and Λk is a N × N diagonal masking matrix with its diagonal
entries generated by b1 b2 , where b1 and b2 are independent
and distributed as [5]

1 w. prob. 0.25


(√

 −1 w. prob. 0.25
2/2 w. prob. 0.8
√
b1 =
and b2 =

−j
w.
prob.
0.25
3 w. prob. 0.2.



j w. prob. 0.25
All results are obtained using a computer with 3.1 GHz Intel
Xeon E31225 and 16 GB RAM. For the random measurement
case, the algorithms are all initialized from the same starting
point
is computed by picking the principal eigenvector of
PM that
2
H
y
a
i=1 i i ai , and the stopping criterion is
2
2

− y − |Ax(r−1) |

y − |Ax(r−1) |

2
2

2
2

AltIRLS
AltGD

100

10-5

10-10

IV. S IMULATION R ESULTS

y − |Ax(r) |

105

MSE

Remark 3.1: In deriving the CRB in Proposition 3.4, we
use no additional assumptions on x. Therefore, Proposition
3.4 works for both real and complex x. Specifically, in the
real case, the phase is actually the sign of x. Here, FL is
also singular with rank deficit equal to one, so we adopt its
pseudo-inverse to compute the CRB. We omit the proof of
rank-1 deficiency of FL , since it follows the line of argument
in Appendix C.

≤ 10−7

or the number of iterations reaching 1000. Furthermore, we
use AltGD with extrapolation for the simulations.
A. Selection of p
Before we do the performance comparison, let us study how
p affects the performance of the proposed methods. In this
example, SNR is fixed at 20 dB, where the SNR is computed
via


||Ax||2
SNR = 10 log10
.
||n||2
We assume that 10% of the data are corrupted by outliers that
are generated from the Gaussian distribution with mean zero

10-15

10-20

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

p

Fig. 2. MSE versus p.

and variance 100, and the remaining elements in n are zero.
K = 8 masks are employed to generate the measurements.
Fig. 2 shows the MSE as a function of p. Note that, when
p is smaller than 1, in order to achieve convergence in a
non-convex setting, we initialize our methods in a two-step
way. Specifically, for each of the proposed methods, when
0.6 < p < 1, we use the spectrum output to initialize our
methods with p = 1.3 and 100 iterations, and then choose the
corresponding output to do another 100 iterations with p = 1
and use the output as a final starting point. When p ≤ 0.6,
we take an additional intermediate step with p = 0.7 and 100
iterations to gradually stage the initialization process for our
approaches, since in this p-regime the subproblem is strictly
non-convex. It is observed that generally, p ≤ 1 provides
better performance than p > 1, especially for AltIRLS. When
p < 0.5, AltGD suffers severe performance degradation while
AltIRLS still performs well. Our understanding is that firstorder algorithms are in general more sensitive to problem
conditioning, and a small p can lead to badly-conditioned
W(r) . Working with p < 1 usually requires much more
careful initialization. In our experience, p ≈ 1.3 strikes a
good compromise between numerical stability and estimation
accuracy.
To illustrate the performance provided by our proposed
methods with p < 1 over the state-of-art algorithms including
TWF, WF, TAF, AltMinPhase, MTWF, and GS, we present the
following example, where the parameters are the same as for
Fig. 2, except that 20% of the data are corrupted by outliers.
Note that among these competitors, TAF, AltMinPhase and GS
share the same magnitude measurement model in (1) as ours,
but WF is designed under the energy measurement model (i.e.,
2
ym = |aH
m x| + nm ). In addition, TWF and MTWF adopt a
2
2
Poisson model with λ = |aH
i x| . Hence, y is fed to WF,
TWF and MTWF, where the squaring is element-wise. We set
p = 0.4 for the proposed methods, and perform Monte-Carlo
trials to estimate MSE = 10 log10 (||x̂ − x||22 ). We show the
results of 1,000 random trials in Fig. 3. One can see that for
most trials the benchmark algorithms fail to give reasonable

7

AltIRLS

AltGD
100

600
400

50

200
0
-200

-150

-100

-50

0

0
-200

-150

MSE (dB)
WF
150

150

100

100

50

50

0
-200

-150

-100

-50

-100

-50

0

MSE (dB)
TWF

0

0
-200

MSE (dB)
GS

-150

-100

-50

0

MSE (dB)
AltMinPhase

150

400

100
200

that admit closed-form PDF expression, that is, α = 1 and
α = 2 which respectively yield the Cauchy and Gaussian
distributions. For α < 2, φ(t; α, β, γ, µ) possesses heavy tails,
and thus is considered suitable for simulating impulsive noise.
The parameter α controls the density of impulses - smaller α’s
correspond to more outliers. In the following, we set α = 0.8,
and the other parameters are β = 0, γ = 2 and µ = 0, resulting
in a symmetric α-stable (SαS) distribution with zero-shift.
3) GMM: In our simulations, we also use a two-component
GMM to generate impulsive noise, whose PDF is as follows:


2
X
|n|2
c
p i
exp
−
(35)
p(n) =
2σi2
2πσi2
i=1
where σi2 is the variance of the ith term, and 0 ≤ ci ≤ 1 is the
probability of occurrence of the ith component. Note that we
have c1 +c2 = 1, and we set c2 < c1 and σ22 > σ12 – therefore,
the second component corresponds to outliers. In this example,
we choose c1 = 0.9, c2 = 0.1, σ12 = 0.1 and σ22 = 100, which
corresponds to a situation where strong outliers are present.

50
0
-200

-150

-100

-50

0

0
-200

-150

MSE (dB)
MTWF

-100

-50

0

MSE (dB)
TAF
150

1000
100
500
0
-200

50

-150

-100

-50

0

0
-200

MSE (dB)

-150

-100

-50

0

MSE (dB)

Fig. 3. Signal recovery performance comparison.

results. However, our methods produce much lower MSEs,
which indicates that our proposed methods are efficient in
suppressing outliers.
B. Statistical Performance Comparison
We now compare the MSE performance as a function of
SNR. The results are averaged over 500 Monte-Carlo trials.
To simulate outliers, we use the Laplacian distribution, the
α-stable distribution, and a Gaussian mixture model (GMM)
to generate heavy-tailed noise n, respectively. Specifically, the
outlier generating process is summarized below:
1) Laplacian: The probability density function (PDF) of
the Laplacian distribution is given in (52).
2) α-stable: The PDF of an α-stable distribution is generally not available, but its characteristic function can be written
in closed-form as
φ(t; α, β, c, µ) = exp (jtµ − γ α |t|α (1 − jβsgn(t)Φ(α)))
where Φ(α) = tan (απ/2), 0 < α ≤ 2 is the stability
parameter, −1 ≤ β ≤ 1 is a measure of asymmetry, γ > 0 is
a scale factor which measures the width of the distribution
and µ is a shift parameter. There are two special cases

Figs. 5, 6 and 7 show the MSEs of x̂ under Laplacian, SαS
and GMM noise, respectively. In the simulations, the number
of measurements is set to be M = 8N . Note that for SαS
and GMM noise, we set p = 1.3 for AltIRLS and AltGD,
while for Laplacian noise we let p = 1 which corresponds
to the maximum likelihood estimator. In the Laplacian noise
case, we include the CRB derived in Proposition 3.2 as a
performance benchmark. It is observed from Fig. 5 that our
approaches outperform the GS algorithm and yield the MSEs
which are closest to the CRB. The performance gap between
the proposed methods and their competitors becomes larger
in the SαS and GMM noise cases since the outliers become
stronger; see Fig. 6 and 7. This indicates that the proposed
algorithms work better in more critical situations, i.e., when
more severe outliers exist. Note that the MTWF performs
slightly better than TWF in GMM noise and its performance
is inferior to our methods in the three noise scenarios.
As a reference, we also plot the MSE performance versus
SNR in Gaussian noise in Fig. 8, where we set p = 2 for
having a maximum likelihood estimator. We compare the performance of the algorithms with the CRB in Proposition 3.4.
One can see that the algorithms perform similarly when the
noise follows an i.i.d. Gaussian distribution. This suggests that
the Gaussian noise case may be an easier case - where all the
algorithms under test perform reasonably.
Fig. 4 compares the CPU time as a function of N , where
SNR is 10 dB, N increases from 8 to 128 and the other
parameters are the same as Fig. 7. It is seen that AltGD is
the fastest1 .
C. Success Rate Performance Comparison
Fig. 9 shows the success rate versus the proportion of
outliers in the measurements, where the GMM is employed
1 Note that in the simulation we use the heuristic step size as stated in
Remark 2.1; if the step size µ = λmax (AH (W(r) )2 A) is used, AltGD
may consume more time.

8

103

102

101

AltIRLS
AltGD
MTWF
TAF
TWF
WF
AltMinPhase
GS
CRB:Laplacian

10-1

MSE

CPU time (s)

100

AltIRLS
AltGD
MTWF
TAF
TWF
WF
AltMinPhase
GS

100
10-2
10-1

10-2

0

50

100

150

200

250

10-3

300

N

5

10

15

20

SNR (dB)

Fig. 4. CPU time versus number of variables N in GMM noise.

Fig. 5. MSE versus SNR in Laplacian noise.

to control the outlier fraction. We vary c2 from 0 to 0.5 and
fix σ12 = 0 and σ22 = 100. SNR is 10 dB. For the proposed
methods, we choose p = 0.4 for comparison. We declare
success when the squared Euclidean distance of the estimate
x̂ from the ground-truth x is smaller than or equal to 10−4 ;
i.e., the success rate is computed as

101

100

10-1

(36)

MSE

500

1 X
success rate =
γi
500 i=1

AltIRLS
AltGD
MTWF
TAF
TWF
WF
AltMinPhase
GS

10-2

where
(
γi =

1,
0,

if kx̂i − xk22 ≤ 10−4
otherwise.

10-3

(37)

Fig. 9 shows that our methods outperform the benchmark
methods by a large margin in terms of success rate. When
the outlier fraction is greater than 0.25, our methods can still
estimate x quite accurately, while the other six competitors
all fail. Note that the performance of AltIRLS is better than
AltGD – but the complexity of AltGD is much lower.
Fig. 10 plots the success rate versus sample complexity
(measured in terms of M/N with N = 16), where SNR is 10
dB and 20% of the data are corrupted by outliers. It is obvious
that our two algorithms have the highest success rates when
M/N ≥ 5, while the other algorithms have very low success
rates even when M/N ≥ 8.

10-4

5

10

15

Fig. 6. MSE versus SNR in SαS noise.

101
AltIRLS
AltGD
MTWF
TAF
TWF
WF
AltMinPhase
GS

100

In this subsection, we showcase the performance using real
caffeine molecule data that is often used to test phase retrieval
algorithms [7]. The caffeine molecule data with size 128 ×
128 is the projection of electron density of a 3-D Caffeine
molecule’s density map onto the xy-plane. The objective is to
reconstruct the data from the magnitude of its masked Fourier
transform. Here, we use the same measuring process as in the
previous simulations and K = 8 masks are used to generate
measurements. We compare the performance of AltGD, TAF,

MSE

10-1

D. Performance Comparison on Caffeine Molecule Image
Data

20

SNR (dB)

10-2

10-3

10-4

5

10

15

SNR (dB)

Fig. 7. MSE versus SNR in GMM noise.

20

9

AltIRLS
AltGD
MTWF
TAF
TWF
WF
AltMinPhase
GS
CRB:Gaussian

100

MSE

MTWF, and TWF. We also include the block incremental (BI)
implementation of AltGD (cf. Section II.C2), which is referred
to as BI-AltGD. The outliers are generated from GMM with
c2 = 0.3, σ12 = 0 and σ22 = 100, which indicates that there
are 30% of the data are corrupted by outliers. SNR is 0 dB
and the data is normalized with unit norm. Fig. 11 plots the
retrieved molecule’s density map using AltGD, BI-AltGD with
p = 1.3, TAF, MTWF and TWF. It is seen in Fig. 11 that
our schemes still work well, while the other competitors yield
blurred images.

10-1

E. Performance Comparison with 2D Fourier Measurements

10-2

5

10

15

20

SNR (dB)

Fig. 8. MSE versus SNR in Gaussian noise.

1
AltIRLS
AltGD
MTWF
TAF
TWF
WF
AltMinPhase

0.9
0.8

Success rate

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

Outliers fraction

Fig. 9. Success rate of exact recovery versus outlier fraction in GMM noise.

1
AltIRLS
AltGD
MTWF
TAF
TWF
WF
AltMinPhase

0.9
0.8

Success rate

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
2

3

4

5

6

7

8

9

10

M/N

Fig. 10. Success rate of exact recovery versus M/N under GMM noise.

In the above simulations, the measurement vectors are
generated from coded diffraction pattern (CDP) with random
masks. The success of the existing algorithms such as the
well-known PhaseLift, WF and its variants are all based on
the randomness of measurements. However, to the best of
our knowledge, the CDP model is still considered to be impractical, as there has been no device supporting the masking
technique yet. In practice, most devices record the intensity
of the Fourier transform of the object, where the intensity
is usually proportional to the magnitude of the oversampled
2D Fourier measurements with 2× oversampling. In this case,
the measurement matrix A is a Kronecker product of two
oversampled Fourier matrices.
With 2D Fourier measurements, the most successful method
is Fienup’s hybrid-input-output (HIO) algorithm [3], which
has been observed to give empirically reliable phase retrieval
performance, especially in the noiseless case. In the noisy case,
HIO can also provide good initialization for other algorithms,
e.g., GS. Therefore, HIO and HIO-initialized GS are employed
as baselines in this subsection.
In this simulation, we compare the performance of AltGD,
HIO and GS under 2D Fourier measurements in the presence
of outliers, where GMM-type outliers are generated with
c2 = 0.1, σ12 = 0 and σ22 = 100. Each entry of the signal
matrix X ∈ R16×16 is generated from the normal distribution.
After zero-padding around X, we have X̃ ∈ R32×32 . Then the
measurements are obtained by first calculating 2D fast Fourier
transformation (FFT) on X̃ and then taking the magnitude
of 2D FFT(X̃), where the sizes of left and right Fourier
matrices are both 32 × 32. Note that SNR is 10 dB and is
computed as 10 log10 (kX̃k2F /kNk2F ) where N ∈ R32×32 is
the noise matrix. For AltGD and GS, HIO with 5000 iterations
is employed to provide an initial estimate of X, and then
an additional 5000 iterations are used in AltGD and GS,
respectively. Moreover, HIO with 10000 iterations is included
as a baseline. For AltGD, we choose p = 1.3. The MSEs
obtained from 1000 Monte-Carlo tests are plotted in Fig. 12.
It is seen that HIO does not perform very well, and is inferior
to AltGD and GS. AltGD has better performance than GS,
since most of its MSEs are concentrated around −35 dB, while
those obtained by GS are located around −20 dB.
Finally, Fig. 13 shows the recovered image of 128 × 128
Cameraman from 2D Fourier measurements under GMM
noise, where SNR is 12 dB and the noise is generated in the
same way as for Fig. 12. The number of iterations used for

10

20

20

20

40

40

40

60

60

60

80

80

80

100

100

100

120

120
20

40

60

80

100

120

120

20

40

(a) Original

60

80

100

120

20

(b) BI-AltGD

20

20

40

40

40

60

60

60

80

80

80

100

100

100

120

120
40

60

80

100

60

80

100

120

80

100

120

(c) AltGD

20

20

40

120

120

20

40

(d) TAF

60

80

100

120

(e) TWF

20

40

60

(f) MTWF

Fig. 11. Retrieving molecule image in SαS noise.

50

0
-40

AltGD

-30

-20

-10

0

10

20

30

MSE (dB)
50

0
-40

GS

-30

-20

-10

0

10

20

30

MSE (dB)
40
HIO

20
0
-40

-30

-20

-10

0

10

20

(a) Original

(b) AltGD

(c) GS

(d) HIO

30

MSE (dB)

Fig. 12. Signal recovery performance comparison with 2D Fourier measurements in GMM noise.

initialization and refinement are all 5000. For AltGD, we set
p = 0.6. It is seen that the image recovered from AltGD is
much more clear than those recovered by HIO and GS.
V. C ONCLUSION
In this paper, we considered phase retrieval in the presence
of grossly corrupted data – i.e., outliers. We formulated this
problem as an `p fitting problem, where 0 < `p < 2, and
provided an algorithmic framework that is based on two-block

Fig. 13. Retrieving Cameraman image in GMM noise with 2D Fourier
measurements.

11

inexact alternating optimization. Two algorithms, namely,
AltIRLS and AltGD, were proposed under this framework.
Although the algorithms cannot be analyzed using standard
convergence results for alternating optimization due to a
nonconvex constraint, we managed to show that the algorithms
converge to a KKT point. The tools used for convergence
analysis can also be used for analyzing convergence of other
types of algorithms that involve non-convex constraints and
inexact alternating optimization. Pertinent CRBs were derived
for the noisy measurement models considered. Simulations
showed that the proposed algorithms are promising in dealing
with outliers in the context of phase retrieval using both
random measurements and 2D Fourier measurements, for
simulated and real image data.

≥ g(x̃(r+1) , {ũ(r)
m })
≥
≥

(rj )
j)
g(x̃, {ũm
}) ≥ g(x̃(rj +1) , {ũ(r
m })

≥
≥
≥

Let us denote

=

f (x, u) =

2
|ym um − aH
m x| + 

p/2

m=1

Therefore, our problem can be written as
min

f (x̃, {ũm }) =

kym ũm −

p/2
+

s.t. kũm k22 = 1, ∀m.

(r)
2
(r)
wm
kym ũ(r)
m − Ãm x̃k2 + φ(wm ) .

(39)
g(x̃, {ũ(r)
m })

= arg min
x̃

ũ(r+1)
= arg
m

f (x̃(r+1) , {ũm })

min

kũm k22 =1, ∀m

(40)
(41)

Let us assume the boundedness of x̃(r) and ũ(r) and
which will be shown later. Under this assumption, the
(r)
gradients of g(x̃, {ũm }) and f (x̃, {ũm }) exist. According to
Lemma 2.1, it is easily seen that
(r)
f (x̃, {ũ(r)
m }) ≤ g(x̃, {ũm }) ∀x̃,

f (x̃ , {ũ(r)
m })
(r)
∇x̃ f (x̃ , {ũ(r)
m })

=
=

(r)

g(x̃ , {ũ(r)
m })
(r)
∇x̃ g(x̃ , {ũ(r)
m }).

(42a)
(42b)
(42c)

Now, we have
f (x̃

(44e)

(47)

(48)

f (x̃? , {um }) ≥ f (x̃? , {u?m }).

(49)

ũ?m = 0, ∀m

= g(x̃

(r)

, {ũ(r)
m })

(43a)

(50)

where λm ∈ R2 are dual variables. Combining (47) and (50),
it follows that every limit point of the solution sequence is a
KKT point.
Now, we rigorously show the boundedness of x̃(r) and ũ(r)
and W(r) . Let us first show that W(r) is bounded given
finite x̃(r) and ũ(r) . This is relatively easy, since y and A
are always bounded and thus W(r) is also bounded when
0 < p < 2 and  > 0 by definition (cf. Eq. (8)). Now, we
show that x̃(r) and ũ(r) are bounded. Since the algorithm is
essentially majorization minimization, there is no risk for x̃(r)
being unbounded, if the initialization is bounded, A has full
column rank, and y is bounded. To explain, one can see that
ũ(r) − Ãx̃(0) kpp ≥ ky

ũ(r) − Ãx̃(r) kpp

≥ |kykpp − kÃx̃(r) kpp |
, {ũ(r)
m })

(44d)

Therefore, ũ?m also satisfies the partial KKT condition

ky
(r)

(44c)

(rj )
f (x̃(rj ) , u) ≥ f (x̃(rj ) , {ũm
})

∇ũm f (x̃? , {ũ?m }) + λ?m

(r)
wm ,

(r)

(44b)

and thus

Now our algorithm can be re-expressed as
x̃

(44a)

(45)

∇x̃ f (x̃? , {ũ?m }) = 0.



m=1

(r+1)

j)
f (x̃
, {ũ(r
m })
j +1)
f (x̃(rj +1) , {ũ(r
})
m
(rj+1 )
(rj+1 )
f (x̃
, {ũm
})
(rj+1 )
(rj+1 )
g(x̃
, {ũm
}),

Similarly, by the update rule in (41), we have

(38)
Note that in the new expression (38), all the variables are realvalued. Accordingly, we may define
g(x̃, {ũ(r)
m }) =

(43d)

By (42c), we immediately have

m=1

M 
X

(43c)

The inequality in (45) means that x̃? is blockwise minimizer of
g(x̃, {ũ?m }). Therefore, it satisfies the partial KKT condition
w.r.t. x̃, i.e.,
(46)
∇x̃ g(x̃? , {ũ?m }) = 0.

x̃ = [Re{xT }, Im{xT }]T ,


H
Re{aH
m }, −Im{am }
.
Ãm =
Im{aH
Re{aH
m },
m}

{ũm },x̃

(rj +1)

g(x̃, {ũ?m }) ≥ g(x̃? , u? ).

ũm = [Re{um }, Im{um }]T ,

Ãm x̃k22

(43b)

where (44d) is obtained by the fact that rj+1 ≥ rj + 1
since rj indexes a subsequence. Taking j → ∞, and by the
boundedness of W(r) and continuity of g(·), we see that

.

To simplify the analysis, let us define

M 
X

f (x̃
, {ũ(r)
m })
(r+1)
f (x̃
, {ũ(r+1)
})
m

where (43a) follows (42b), (43b) is obtained because of (40),
(43c) holds since we have the property in (42a), and (43b) is
obtained by the fact that the subproblem w.r.t. u is optimally
solved, i.e., (41).
Assume that {rj } denotes the index set of a convergent sub(r )
sequence, and that {x̃(rj ) , {ũmj }} converges to (x̃? , {ũ?m })
Then, we have

A PPENDIX A
P ROOF OF P ROPOSITION 2.1
M
X

(r+1)

≥

kÃx̃(r) kpp

−

kykpp .

(51)

12

The first inequality holds because block upper bound minimization algorithm always decreases the objective [40]. The
second inequality is due to the triangle inequality for k · kpp
when 0 < p < 1 [39]; for p > 1, one can apply the same
argument to k · kp . From the above, one can see that
kÃx̃(r) kpp

≤ ∞.

Therefore, if A does not have a null space, then x̃(r) has to be
bounded. Consequently, W(r) is always bounded from above.
The next step is to show that the whole sequence converges
to K, the set of all the KKT points. This is relatively
straightforward. Let us assume that the whole sequence does
not converge to K. This means that there exists a subsequence
which does not converge to K. The boundedness of x̃(r) , ũ(r)
implies that every subsequence of {x̃(r) , ũ(r) } has a limit
point. We have just shown that every limit point is a KKT
point, and thus this is a contradiction. Therefore, the whole
sequence has to converge to K.
A PPENDIX B
CRB FOR L APLACIAN N OISE
The likelihood function for Laplacian noise is given by [36][37]
)
( √
M
Y
2
1
H
√
yi − |ai x|
(52)
exp −
p(y; x) =
σn
2σn
i=1
where the noise variance is
is

σn2 ,

and its log-likelihood function

√
ln p(y; x) = −M ln( 2σn ) −

√

M

2X

σn

yi − |aH
i x| . (53)

i=1

The vector of unknown parameters for complex-valued x is
β = [Re{x1 } · · · Re{xN }, Im{x1 } · · · Im{xN }]T . (54)
It is worth noting that σn2 is actually an unknown parameter
which should be considered as a part of β. However, since σn2
is uncorrelated with the real and imaginary parts of xi , their
mutual Fisher information is zero. It will not impact the final
CRB formula for x. For this reason, we do not include σn2 in
β. Thus, the FIM can be partitioned into four parts, i.e.,


FL,rr FL,ri
FL,c =
(55)
FL,ir FL,ii
where

[FL,c ]m,n = E

(
sgn(a) =

(56)

The partial derivative of ln p(y; x) with respective to βm is
√ M
2 X ∂ yi − |aH
∂ ln p(y; x)
i x|
=−
∂βm
σn i=1
∂βm
√ M
H
2 X yi − |aH
i x| ∂|ai x|
=
H
σn i=1 yi − |ai x| ∂βm
√ M
2X
∂|aH
i x|
=
sgn(yi − |aH
(57)
i x|)
σn i=1
∂βm

1,

a>0

−1,

a<0

(58)

and


[Re ai aH

i x ]m


, for βm = Re{xm }

H
H
|a x|
∂|ai x|
 i H
=

∂βm
[Im ai ai x ]m


, for βm = Im{xm }.

|aH
i x|

(59)

Substituting (57) into (56), we have
[FL,c ]m,n =

M M
H
2 X X ∂|aH
i x| ∂|ai x|
σn2 i=1 j=1 ∂βm ∂βn


H
× E sgn(yi − |aH
i x|) · sgn(yj − |aj x|) .
(60)

Next we compute the value of E[sgn(yi − |aH
i x|)sgn(yj −
H
|aH
x|)].
For
notational
simplicity,
let
s
=
sgn(y
i
i − |ai x|).
j
It is obvious that when i = j, we have
E [si sj ] = 1.

(61)

For i 6= j, we first write
E [si sj ] = Pr(si = sj ) × (+1) + Pr(si 6= sj ) × (−1)
= 2Pr(si = sj ) − 1

(62)

where Pr standards for the probability. Then the value of
Pr(si = sj ) is computed as
Pr(si = sj ) = Pr(si = 1|sj = 1)Pr(sj = 1)
+ Pr(si = −1|sj = −1)Pr(sj = −1)
= Pr(si = 1)Pr(sj = 1)
+ Pr(si = −1)Pr(sj = −1)
= 0.5

(63)

where Pr(si = 1) = Pr(si = −1) = 0.5 and the second
equation follows by independence of si and sj when i 6= j.
Substituting (63) into (62) and using (61) yields
(
1, i = j
E [si sj ] =
(64)
0, i 6= j.
Substituting (64) into (60), we obtain
[FL,c ]m,n =



∂ ln p(y; x) ∂ ln p(y; x)
.
∂βm
∂βn

where

M
H
2 X ∂|aH
i x| ∂|ai x|
.
σn2 i=1 ∂βm ∂βn

(65)

Now using (59) and (65), the four sub-FIMs can be easily
derived as
2
FL,rr = 2 Re{AH diag(Ax)} · diag(|Ax|−2 )
σn
× Re{AH diag(Ax)}T
(66)
2
FL,ii = 2 Im{AH diag(Ax)} · diag(|Ax|−2 )
σn
× Im{AH diag(Ax)}T
(67)
2
FL,ri = 2 Re{AH diag(Ax)} · diag(|Ax|−2 )
σn

13

× Im{AH diag(Ax)}T
FL,ir =

FTL,ri .

(68)
(69)

Combining (66)-(69), we obtain
FL,c =

2
GL,c diag(|Ax|−2 ) GTL,c
σn2

(70)



Re{AH diag(Ax)}
=
.
Im{AH diag(Ax)}

(71)

where
GL,c

This completes the proof of Proposition 3.1.
A PPENDIX C
P ROOF OF R ANK P ROPERTY OF FL,c AND FL,r
We note here that FL,c is derived under the assumption
−2
of nonzero aH
) is full rank. As a
m x. Therefore, diag(|Ax|
result, computing the rank of FL,c is equivalent to computing
the rank of GL,c . To this end, define a nonzero vector v =
[ v1T v2T ]T ∈ R2N , which leads to
GTL,c v = Re{AH diag(Ax)}T v1 + Im{AH diag(Ax)}T v2 .
(72)
Now let u = v1 + jv2 , then
n
H o
GTL,c v = Re AH diag(Ax) u
= Re {(Ax)∗

(Au)}

(73)

which equals to zero if and only if
u = jx

(74)

n
o
2
GTL,c v = Re j |Ax| = 0.

(75)

i.e.,

This means that there is only one direction v =
[−Im{x}T Re{x}T ]T , which is non-zero, lies in the null space
of GL,c , thus also in the null space of FL,c .
In the real x case, similar to the proof of FL,c , it suffices to
show that there is no nonzero vector v ∈ RN making FL,r v =
0. It is easy to see that
GTL,r v = Re {(Ax)∗

(Av)} .

(76)

Since v is real-valued, it cannot be set to v = jx to make
(76) zero. Given a nontrivial A, it is impossible to find a v
such that Av = 0 except v = 0. Therefore, FL,r is full rank.
This completes the proof.
R EFERENCES
[1] C. Qian, X. Fu, N. D. Sidiropoulos, and L. Huang, “Inexact alternating optimization for phase retrieval with outliers,” Proceeding of 24th European
Signal Processing Conference (EUSIPCO), pp. 1538-1542, Budapest,
2016.
[2] R. Gerchberg and W. Saxton, “A practical algorithm for the determination
of phase from image and diffraction plane pictures,” Optik, vol. 35, pp.
237-246, 1972.
[3] J. R. Fienup, “Phase retrieval algorithms: A comparison,” Applied Optics,
vol. 21, no. 15, pp. 2758-2769, 1982.
[4] P. Netrapalli, P. Jain and S. Sanghavi, “Phase retrieval using alternating
minimization,” IEEE Trans. Signal Process., vol. 63, no. 18, pp 48144826, 2015.

[5] E. J. Candès, T. Strohmer, and V. Voroninski. “PhaseLift: Exact and stable
signal recovery from magnitude measurements via convex programming,”
Communications on Pure and Applied Mathematics, vol. 66, no. 8, pp.
1241-1274, 2013.
[6] P. Hand, “PhaseLift is robust to a constant fraction of arbitrary errors,”
Applied and Computational Harmonic Analysis, vol. 42, no. 3, pp. 550562, 2017.
[7] E. J. Candès, X. Li and M. Soltanolkotabi, “Phase retrieval via Wirtinger
Flow: Theory and algorithms,” IEEE Trans. Info. Theory, vol. 61, no. 4,
pp. 1985-2007, 2015.
[8] Y. Chen and E. Candès, “Solving random quadratic systems of equations is nearly as easy as solving linear systems,” Advances in Neural
Information Processing System, pp. 739-747, 2015.
[9] G. Wang, G. B. Giannakis, and Y. C. Eldar, “Solving systems of
random quadratic equations via truncated amplitude flow,” arXiv preprint
arXiv:1605.08285, 2016.
[10] G. Wang, G. B. Giannakis and J. Chen, “Scalable solvers of random
quadratic equations via stochastic truncated amplitude flow,” IEEE Trans.
Signal Process., vol. 65, no. 8, pp. 1961-1974, 2017.
[11] C. Qian, N. D. Sidiropoulos, K. Huang, L. Huang and H. C. So, “Phase
retrieval using feasible point pursuit: Algorithms and Cramér-Rao bound,”
IEEE Trans. Signal Process., vol. 64, no. 20, pp. 5282-5296, 2016.
[12] C. Qian, N. D. Sidiropoulos, K. Huang, L. Huang and H. C. So, “Least
squares phase retrieval using feasible point pursuit,” Proceedings of the
International Conference on Acoustics, Speech, and Signal Processing
(ICASSP 2016), Shanghai, China, pp. 4288-4292, 2016.
[13] E. J. Candès, Y. C. Eldar, T. Strohmer and V. Voroninski, “Phase retrieval
via matrix completion,” SIAM Review, vol. 57, no. 2 pp. 225-251, 2015.
[14] J. Sigl, “Nonlinear residual minimization by iteratively reweighted least
squares,” Computational Optimization and Applications, vol. 64, no. 3,
pp. 755-792, 2016.
[15] D. S. Weller, A. Pnueli, G. Divon, O. Radzyner, Y. C. Eldar and J.
A. Fessler, “Undersampled phase retrieval with outliers,” IEEE Trans.
Computational Imaging, vol. 1, no. 4, pp. 247-258, 2015.
[16] S. Chen and H. Lu, “Noise characteristic and its removal in digital radiographic system,” Proceeding of 15th World Conference on Nondestructive
Testing, Roma 2000.
[17] I. Frosio and N. A. Borghese, “Statistical based impulsive noise removal
in digital radiography,” IEEE Trans. Medical Imaging, vol. 28, no. 1, pp.
3-16, 2009.
[18] H. Zhang, Y. Chi and Y. Liang, “Provable non-convex phase retrieval
with outliers: Median truncated wirtinger flow,” Proceeding of International Conference on Machine Learning (ICML), New York, 2016.
[19] J. N. Cederquist and C. C. Wackerman, “Phase-retrieval error: A lower
bound,” Journal of the Optical Society of America A, vol. 4, no. 9 pp.
1788-1792, 1987.
[20] R. Balan, “The Fisher information matrix and the CRLB in a nonAWGN model for the phase retrieval problem,” Proceeding of 2015
Internat. Conf. on Sampl. Theory and Applications (SampTA), pp. 178182, Washington, DC, 2015.
[21] R. Balan, “Reconstruction of signals from magnitudes of redundant
representations: The complex case,” Foundations of Computational Mathematics, pp. 1-45, 2013.
[22] A. S. Bandeira, J. Cahill, D. G. Mixon and A. A. Nelson, “Saving phase:
Injectivity and stability for phase retrieval,” Applied and Computational
Harmonic Analysis, vol. 37, no. 1, pp. 106-125, 2014.
[23] X. Fu, K. Huang, B. Yang, W.-K. Ma and N.D. Sidiropoulos, “Robust
volume minimization-based matrix factorization for remote sensing and
document clustering,” IEEE Trans. Signal Process., vol. 64, no. 23, pp.
6254-6268, 2016.
[24] Y. Liu, Y. Dai and S. Ma, “Joint power and admission control: Nonconvex `q approximation and an effective polynomial time deflation
approach,” IEEE Trans. Signal Process., vol. 63, no. 14, pp. 3641-3656,
2015.
[25] R. Chartrand and W. Yin, “Iteratively reweighted algorithms for compressive sensing,” Proceeding of IEEE Internat. Conf. Acoust., Speech
and Signal Process., Las Vegas, NV, pp. 3869-3872, 2008.
[26] X. Fu, K. Huang, W.-K. Ma, N. D. Sidiropoulos, and R. Bro, “Joint
tensor factorization and outlying slab suppression with applications,”
IEEE Trans. Signal Process., vol. 63, no. 23, pp. 6315-6328, 2015.
[27] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM journal on Imaging Sciences,
vol.2, no.1, pp. 183–202, 2009.
[28] Y. Dai and Y. Yuan, “A nonlinear conjugate gradient method with a
strong global convergence property,” SIAM Journal on Optimization, vol.
10, no. 1, pp. 177-182, 1999.

14

[29] Y. Nesterov, “A method for unconstrained convex minimization problem
with the rate of convergence O (1/k2 ),” Doklady an SSSR, vol. 269, no.3,
pp. 543-547, 1983.
[30] O. Fercoq and P. Richtàrik, “Accelerated, parallel, and proximal coordinate descent,” SIAM Journal on Optimization, vol. 25, no. 4, pp. 19972023, 2015.
[31] Y. Nesterov, “Efficiency of coordinate descent methods on huge-scale
optimization problems,” SIAM Journal on Optimization, vol. 22, no. 2
pp. 341-362, 2012.
[32] Y. Xu, R. Hao, W. Yin, and Z. Su, “Parallel matrix factorization for
low-rank tensor completion,” Inverse Problems and Imaging, vol. 9, no.
2, pp. 601-624, 2015.
[33] P. Stoica and T. L. Marzetta, “Parameter estimation problems with
singular information matrices,” IEEE Trans. Signal Process., vol. 49, no.
1, pp. 87-90, 2001.
[34] A. O. Hero, III, J. A. Fessler and M. Usman, “Exploring estimator
bias-variance tradeoffs using the uniform CR bound,” IEEE Trans. Signal
Process., vol. 44, pp. 2026-2041, 1996.
[35] K. Huang and N. D. Sidiropoulos, “Putting nonnegative matrix factorization to the test: A tutorial derivation of pertinent Cramér-Rao bounds and
performance benchmarking,” IEEE Signal Processing Magazine, Special
Issue on Source Separation and Applications, vol. 31, no. 3, pp. 76-86,
2014.
[36] S. A. Vorobyov, Y. Rong, N. D. Sidiropoulos and A. B. Gershman,
“Robust iterative fitting of multilinear models,” IEEE Trans. Signal
Process., vol. 53, no. 8, pp. 2678-2689, 2005.
[37] S. M. Kay, Fundamentals of statistical signal processing: Detection
theory, Upper Saddle River, NJ: Prentice-Hall, 1998.
[38] S. Basu and Y. Bresler, “The stability of nonlinear least squares problems
and the Cramer-Rao bound,” IEEE Transactions on Signal Processing,
vol. 48, no. 12, pp. 3426-3436, 2000.
[39] R. Chartrand and V. Staneva, “Restricted isometry properties and
nonconvex compressive sensing,” Inverse Problems, vol. 24, no. 3, pp.
035020, 2008.
[40] M. Razaviyayn, M. Hong, and Z.-Q. Luo, “A unified convergence
analysis of block successive minimization methods for nonsmooth optimization,” SIAM Journal on Optimization, vol. 23, no. 2, pp. 1126-1153,
2013.

