Coordinate Descent Algorithms for Phase Retrieval

I Introduction

Phase retrieval refers to the recovery of a complex-valued signal from only intensity or squared-magnitude measurements of its linear transformation . It has been a very active field of research because of its wide applicability in science and engineering, which include areas of optical imaging , crystallography , electron microscopy , neutron radiography , digital communications , astronomy and computational biology . The first model for phase retrieval investigates the problem of recovering a signal from the squared-magnitude of its Fourier transform. To address various applications, the power spectrum measurement model has been extended to different formulations, including the short-time Fourier transform , coded diffraction patterns , and random measurements , –. Nevertheless, in all these models, observations of the signal-of-interest (SOI) are obtained via a linear mapping, and we can only measure the intensity.

Early approach to phase retrieval is based on error reduction, which includes the most representative Gerchberg-Saxton (GS) algorithm and its modified version proposed by Fienup , as well as other variants –. In essence, the error reduction techniques apply the concept of alternating projection. That is, at each iteration, the current SOI estimate is projected onto one constraint set such that the magnitudes of its linear mapping match the observations, and then the signal is projected onto another constraint set to conform to the a priori knowledge about its structure . This methodology works well in practice but its convergence is unclear because projection onto nonconvex sets is involved. Recently, the guarantee of convergence to global solution for the GS algorithm is proved under the condition of resampling . The number of measurements required in the resampled GS scheme is on the order of $N\log^{3}N$ with $N$ being the signal length. Nevertheless, it will be clear later that this sampling complexity is not optimal compared with other advanced methods.

In fact, phase recovery corresponds to a nonconvex optimization problem. To be specific, it requires solving a system of quadratic equations, or equivalently, minimizing a multivariate fourth-order polynomial, which is generally known to be NP-hard . The convex relaxation based methods, including PhaseLift and PhaseCut , relax the original nonconvex problem into a convex program. The PhaseLift converts the quadratic equations into linear ones by lifting the $N$-dimensional signal vector to an $N\times N$ rank-one matrix. Then it approximates the minimum rank problem using trace norm minimization, which is convex and can be solved by semidefinite programming (SDP). The sampling complexity of PhaseLift is $\mathcal{O}(N\log N)$, which is lower than that of the resampled alternating projection method and is nearly optimal . On the other hand, the PhaseCut recasts phase retrieval as a quadratically constrained quadratic program (QCQP) which is then approximately solved via semidefinite relaxation . It has similar sampling complexity to PhaseLift and both exhibit good retrieval performance. However, the computational load of the SDP based methods is very high, especially when the signal length or number of observations is large, since the PhaseLift and PhaseCut involve matrix variables with $\mathcal{O}(N^{2})$ and $\mathcal{O}(M^{2})$ elements, respectively, where $M$ is the measurement number. As a result, the convex relaxation approach cannot deal with large-scale problems.

To circumvent the high computational requirement, Wirtinger flow (WF) , which is essentially a gradient descent technique for complex-valued variables, is developed for minimizing the nonconvex quartic polynomial. In general, the gradient method is only guaranteed to converge to a stationary point of a nonconvex objective function. In other words, it can trap in a saddle point or local minimum. That is to say, convergence to the global solution is not guaranteed for general nonconvex optimization problems using the gradient descent. Surprisingly, when initiated via a spectral method and the sample size is $\mathcal{O}(N\log N)$, Candès $et$ $al.$ prove that the WF algorithm converges to the global solution at a geometric rate with high probability. The truncated WF further enhances the recovery performance by adaptively selecting a portion of measurements at each iteration while the optimal stepsize for convergence rate acceleration has been derived in . Still, the convergence speed of the gradient-based WF approach is not fast.

In many applications, the SOI is sparse or only contains a few nonzero entries in some basis. Recovering a sparse signal from the intensity-only measurements is called quadratic compressed sensing. Like classic compressed sensing based on linear measurements , the sampling complexity of phase retrieval can be reduced by exploiting sparsity. Several above-mentioned phase recovery schemes for non-sparse signals have been adapted to handle sparse SOIs. For example, performing hard-thresholding at each iteration of the GS or Fienup scheme yields the so-called the sparse Fienup algorithm . Similarly, applying a thresholding operationThe thresholding operator can be soft or hard. to the WF method elicits the thresholded WF algorithm . Both will yield desirable solution because either soft- or hard-thresholding enforces the signal to be sparse. Furthermore, by borrowing the idea from orthogonal matching pursuit , a greedy algorithm is designed for sparse phase retrieval in .

In this work, we develop effective and computationally efficient algorithms with faster convergence rate for minimizing the nonconvex quartic polynomial in phase retrieval. Our approach is based on coordinate descent (CD), which adopts the strategy of “one at a time” . That is, CD solves a multivariate minimization problem by successively finding a single unknown at each iteration while keeping the remaining variables fixed. According to different rules for coordinate selection, our scheme includes three variants, namely, cyclic, randomized, and greedy CDs. One motivation using CD for phase retrieval is that the exact minimizer of each coordinate is easily obtained by finding the roots of a univariate cubic equation. It is believed that the proposed methodology provides a new path to solve phase retrieval and related problems.

We summarize the contributions of this paper as follows.

An algorithmic framework including cyclic, randomized, and greedy CDs, is proposed to solve the quartic polynomial minimization for phase retrieval. The CD algorithm is computationally simple and converges much faster than the gradient descent methods such as WF and its variants .

Theoretically, we prove that the CD globally converges to a stationary point of the nonconvex problem, where the gradient is non-Lipschitz continuous. It is worth pointing out that the proof is nontrivial because the existing convergence analyses of CD assuming convexity and Lipschitz continuity are not applicable to our problem.

It is proved that the randomized CD locally converges to the global minimum at a geometric rate with high probability using $\mathcal{O}(N\log N)$ measurements.

Currently, the applications of phase retrieval mainly focus on imaging. Here, we open up a new use of phase retrieval for blind equalization in digital communications, i.e., removing the adverse effect induced by channel propagation.

II CD for Phase Retrieval

When the noise is independent and identically distributed (i.i.d.) and Gaussian, the LS estimate given by (2) is equivalent to the maximum likelihood solution. Nevertheless, the optimization problem of (2) is not easy to solve because it is not only nonlinear but also nonconvex.

II-B Outline of CD

To derive the CD, we first analyze the structure of the objective function in (2). The $m$th ($m=1,\cdots,M$) term in (2) is

Note that $\bar{\boldsymbol{A}}_{m}$ is symmetric due to ${\rm Re}(\boldsymbol{A}_{m})={\rm Re}(\boldsymbol{A}_{m})^{T}$ and ${\rm Im}(\boldsymbol{A}_{m})=-{\rm Im}(\boldsymbol{A}_{m})^{T}$ because $\boldsymbol{A}_{m}$ is Hermitian. It is also not difficult to see $\boldsymbol{x}^{H}\boldsymbol{A}_{m}\boldsymbol{x}=\bar{\boldsymbol{x}}^{T}\bar{\boldsymbol{A}}_{m}\bar{\boldsymbol{x}}$. Denoting the quadratic form as

and the original optimization problem of (2) becomes

The objective function $f(\bar{\boldsymbol{x}})$ is a multivariate quartic polynomial of $\bar{\boldsymbol{x}}=[\bar{x}_{1},\cdots,\bar{x}_{2N}]^{T}$ since $q_{m}(\bar{\boldsymbol{x}})$ is quadratic. Minimizing multivariate fourth-order polynomial is known to be NP-hard in general . In this work, we exploit the coordinate update strategy to minimize $f(\bar{\boldsymbol{x}})$. CD is an iterative procedure that successively minimizes the objective function along coordinate directions. Denote the result of the $k$th iteration as $\bar{\boldsymbol{x}}^{k}=[\bar{x}_{1}^{k},\cdots,\bar{x}_{2N}^{k}]^{T}$. In the $k$th iteration, we minimize $f$ with respect to the $i_{k}$th ($i_{k}\in\{1,\cdots,2N\}$) variable while keeping the remaining $2N-1$ variables $\{\bar{x}_{i}^{k}\}_{i\neq i_{k}}$ fixed. This is equivalent to performing a one-dimensional search along the $i_{k}$th coordinate, which can be expressed as

where $\boldsymbol{e}_{i_{k}}$ is the unit vector with the $i_{k}$th entry being one and all other entries being zero. Then $\bar{\boldsymbol{x}}$ is updated by

which implies that only the $i_{k}$th component is updated:

while other components remain unchanged. Since $\bar{\boldsymbol{x}}^{k}$ is known, $f\left(\bar{\boldsymbol{x}}^{k}+\alpha\boldsymbol{e}_{i_{k}}\right)$ is a univariate function of $\alpha$. Thus, (9) is a one-dimensional minimization problem. We will detail how to solve it in the next subsection. Now one reason why we convert the complex-valued problem into real is clear: this makes the scalar minimization problem of (9) real-valued and easier to solve. Otherwise, we still face a problem with a complex number, which in fact is a two-dimensional optimization on the complex plane. The CD is outlined in Algorithm 1.

There are several fashions to select the coordinate index $i_{k}$. The following three selection rules are considered in this paper.

Cyclic rule: $i_{k}$ first takes 1, then 2 and so forth through $2N$. The process is then repeated starting with $i_{k}=1$ again. That is, $i_{k}$ takes value cyclically from $\{1,\cdots,2N\}$. Every $2N$ iterations are called one cycle or sweep. The cyclic rule is similar to the Gauss-Seidel iterative method for solving linear systems of equations , where each coordinate is updated using a cyclic order.

Random rule: $i_{k}$ is randomly selected from $\{1,\cdots,2N\}$ with equal probability.

is the partial derivative of $f(\bar{\boldsymbol{x}})$ with respect to $\bar{x}_{i}$, i.e., the $i$th component of the full gradient

The greedy rule is also called Gauss-Southwell rule . Obviously, it chooses the coordinate with the largest (in absolute value) partial derivative. Hence, computing the full gradient is required at each iteration while there is no need for the cyclic and random rules. We refer the three CD methods with cyclic, random, and greedy rules to as CCD, RCD, and GCD, respectively. It will be seen later that the GCD converges faster than CCD and RCD at the expense of the extra full gradient calculation.

We call every $2N$ iterations of the CD as one cycle. Based on Wirtinger calculus , the gradient of $f$ with respect to the complex vector $\boldsymbol{x}$ is computed as

with a complexity of $\mathcal{O}(MN)$. The gradient of $f$ with respect to the real vector $\bar{\boldsymbol{x}}$ is an expanded form of $\nabla f(\boldsymbol{x})$:

II-C Closed-Form Solution of Coordinate Minimization

Employing (8), $\varphi(\alpha)$ is expressed as

where $c_{2,i}^{m}$, $c_{1,i}^{m}$, and $c_{0}^{m}$ are the coefficients of the univariate quadratic polynomial. Note that the constant $c_{0}^{m}$ has no relation to $i$. According to (4), the coefficients of the quadratic term can be simplified to

Using (4) and recalling $\boldsymbol{A}_{m}=\boldsymbol{a}_{m}\boldsymbol{a}_{m}^{H}$, it is revealed that

Hence, the coefficients of the linear term are computed as

Since $q_{m}(\bar{\boldsymbol{x}}+\alpha\boldsymbol{e}_{i})$ is quadratic, $\varphi_{m}(\alpha)$ of (19) is a univariate quartic polynomial of $\alpha$, which is expressed as

where $\{d_{j,i}^{m}\}_{j=1}^{4}$ and $d_{0}^{m}$ are the coefficients of the polynomial. Note that $d_{0}^{m}$ is not related to $i$. Plugging (20) into (19), we obtain

Since $\varphi(\alpha)=\sum_{m=1}^{M}\varphi_{m}(\alpha)$, it is clear that the coefficients of the quartic polynomial

correspond to the sums of those of $\{\varphi_{m}(\alpha)\}_{m=1}^{M}$, i.e.,

The minimum point of $\varphi(\alpha)$ must be one of stationary points, i.e., the roots of the derivative

Equation (29) refers to finding the roots of a univariate cubic polynomial, which is easy and fast because there is a closed-form solution . Since the coefficients of the cubic equation are real-valued, there are only two possible cases on the roots. The first case is that (29) has a real root and a pair of complex conjugate roots. In this case, the minimizer is the unique real root because the optimal solution of a real-valued problem must be real-valued. The second case is that (29) has three real roots. Then the optimal $\alpha$ is the real root associated with the minimum objective. Once the coefficients of (29) are obtained, the complexity of calculating the roots of a cubic polynomial is merely $\mathcal{O}(1)$. Herein, the second reason why we recast the complex-valued problem into real is clear: by this fashion, it results in root finding of a cubic equation with real coefficients, which has a closed-form solution and is much simpler than the case with complex coefficients.

Computational Complexity: The leading computational cost at each iteration of the CD is calculating the coefficients $\{d_{j,i}\}_{j=1}^{4}$, or equivalently, computing $c_{2,i}^{m}$, $c_{1,i}^{m}$, and $c_{0}^{m}$ with $m=1,\cdots,M$.This is because $\{d_{j,i}\}_{j=1}^{4}$ can be easily calculated from $c_{2,i}^{m}$, $c_{1,i}^{m}$, and $c_{0}^{m}$ according to (26) and (28). From (21), $c_{2,i}^{m}$ is just the squared modulus of $[\boldsymbol{a}_{m}]_{i}$ and can be pre-computed in advance before iteration, which requires $\mathcal{O}(M)$ multiplications for determining all $M$ coefficients $\{c_{2,i}^{m}\}_{m=1}^{M}$. According to (23) and (24), we need to compute $\{\boldsymbol{a}_{m}^{H}\boldsymbol{x}\}_{m=1}^{M}$ in order to obtain $\{c_{1,i}^{m}\}_{m=1}^{M}$ and $\{c_{0}^{m}\}_{m=1}^{M}$. This involves a matrix-vector multiplication $\boldsymbol{A}\boldsymbol{x}$, where the $m$th row of the matrix $\boldsymbol{A}$ is $\boldsymbol{a}_{m}^{H}$, i.e.,

At first glance, the matrix-vector multiplication requires a complexity of $\mathcal{O}(MN)$. However, this complexity can be reduced to $\mathcal{O}(M)$ per iteration for the CD. By observing (9), we know that only one single element changes in two consecutive iterations. Specifically, we have

where $\boldsymbol{A}_{:,j}$ represents the $j$th column of $\boldsymbol{A}$. It is only required to compute $\boldsymbol{A}\boldsymbol{x}^{0}$ before iteration. After that, this matrix-vector product can be efficiently updated from that of the previous iteration by a cheap computation of a scalar-vector multiplication $(x_{j}^{k+1}-x_{j}^{k})\boldsymbol{A}_{:,j}$, which merely costs $\mathcal{O}(M)$ operations. In summary, the complexity of CCD and RCD is $\mathcal{O}(M)$ per iteration. Therefore, the complexity of $2N$ iterations, i.e., a cycle for CCD, is the same as that of the WF method using full gradient descent. While for GCD, an extra cost for computing the full gradient is needed, which results in a complexity of $\mathcal{O}(MN)$.

Initialization and Termination: The spectral method in provides a good initial value for phase retrieval. For Gaussian measurement model and in the absence of noise, we have

Since $\boldsymbol{x}$ is the principal eigenvector of $\boldsymbol{I}+2\boldsymbol{x}\boldsymbol{x}^{H}$ associated with the largest eigenvalue, the principal eigenvectorThe squared norm of the eigenvector is set to $(N\|\boldsymbol{b}\|_{1})/(\sum_{m}\|\boldsymbol{a}_{m}\|^{2})$. of the matrix $\frac{1}{M}\sum_{m=1}^{M}b_{m}\boldsymbol{a}_{m}\boldsymbol{a}_{m}^{H}$, which is an estimate of $\boldsymbol{I}+2\boldsymbol{x}\boldsymbol{x}^{H}$, is taken as the initial value $\boldsymbol{x}^{0}$. More details of the spectral method for initialization can be found in . There are several measures for terminating the CD algorithm. For example, the reduction of the objective function can be used to check for convergence. Specifically, the iteration is terminated when

holds, where $\texttt{TOL}>0$ is a small tolerance parameter. Note that the CD monotonically decreases the objective function, implying $f(\bar{\boldsymbol{x}}^{k})-f(\bar{\boldsymbol{x}}^{k+1})>0$.

III Convergence Analysis

Most existing convergence analyses for CD assume that the objective function is convex and the gradient is Lipschitz continuous . However, the objective function for phase retrieval is quartic and hence nonconvex. As shown in (15) and (16), the gradient is not Lipschitz continuous. Therefore, the available convergence analyses are not applicable to the CD for phase retrieval. In this section, We first prove that the three CD algorithms globally converge to a stationary point from any initial value. Then, it is proved that the sequence of the iterates generated by the RCD locally converges to the global minimum point in expectation at a geometric rate under a mild assumption. This implies that in the absence of noise, the RCD achieves exact phase retrieval under a moderate condition.

We first present two lemmas used in the proof.

is compact, viz. bounded and closed. The iterates of the three CD algorithms, i.e., $\bar{\boldsymbol{x}}^{k}$, $k=0,1,\cdots$, are in the compact set $\mathcal{S}_{f_{0}}$.

Proof: If $\|\bar{\boldsymbol{x}}\|\rightarrow\infty$, then $f(\bar{\boldsymbol{x}})\rightarrow\infty$ since $f(\bar{\boldsymbol{x}})$ is quartic. The converse-negative proposition implies that $f(\bar{\boldsymbol{x}})\leq f_{0}<\infty$ guaranteeing $\|\bar{\boldsymbol{x}}\|<\infty$ for all $\bar{\boldsymbol{x}}\in\mathcal{S}_{f_{0}}$. Hence, $\mathcal{S}_{f_{0}}$ is bounded. Now it is clear that all the points in the sublevel set satisfy $f(\bar{\boldsymbol{x}})\in[0,f_{0}]$ as we also have $f(\bar{\boldsymbol{x}})\geq 0$. Since the mapping $f(\bar{\boldsymbol{x}})$ is continuous and the image $[0,f_{0}]$ is a closed set, the inverse image $\{\bar{\boldsymbol{x}}|0\leq f(\bar{\boldsymbol{x}})\leq f_{0}\}$ is also closed. This completes the proof that $\mathcal{S}_{f_{0}}$ is compact. The CD monotonically decreases $f(\bar{\boldsymbol{x}})$, meaning that $f(\bar{\boldsymbol{x}}^{k})\leq f(\bar{\boldsymbol{x}}^{k-1})\leq\cdots\leq f(\bar{\boldsymbol{x}}^{0})=f_{0}$. Therefore, all the iterates $\{\bar{\boldsymbol{x}}^{k}\}$ must be in the compact set $\mathcal{S}_{f_{0}}$. $\Box$

Lemma 2: On the compact set $\mathcal{S}_{f_{0}}$, the gradient $\nabla f(\bar{\boldsymbol{x}})$ is component-wise Lipschitz continuous. That is, for each $i=1,\cdots,2N$, we have

for all $\bar{\boldsymbol{x}},\bar{\boldsymbol{x}}+t\boldsymbol{e}_{i}\in\mathcal{S}_{f_{0}}$, where $L_{i}>0$ is referred to as the component-wise Lipschitz constant on $\mathcal{S}_{f_{0}}$. Further, it follows

Proof: For $t=0$, both sides of (36) are equal to 0 and (36) holds. For $t\neq 0$, it means that $\|\bar{\boldsymbol{x}}+t\boldsymbol{e}_{i}-\bar{\boldsymbol{x}}\|=t\neq 0$. Since the function

is continuous on the compact set $\mathcal{S}_{f_{0}}$, its minimum over $\mathcal{S}_{f_{0}}$, namely, $L_{i}$, is attained by Weierstrass’ theorem . Then

immediately elicits (36). The following second-order partial derivative

is well defined since $f(\bar{\boldsymbol{x}})$ is twice continuously differentiable, which represents the $(i,i)$ entry of the Hessian matrix:

Noting that $\nabla_{i,i}^{2}f(\bar{\boldsymbol{x}})$ is the partial derivative of $\nabla_{i}f(\bar{\boldsymbol{x}})$ with respect to $\bar{x}_{i}$ and by (36), we obtain

which holds for all $\bar{\boldsymbol{x}}\in\mathcal{S}_{f_{0}}$. Applying Taylor’s theorem and (42), there exists a $\gamma\in$ with $\bar{\boldsymbol{x}}+\gamma t\boldsymbol{e}_{i}\in\mathcal{S}_{f_{0}}$ such that

The component-wise Lipschitz constant $L_{i}$ is not easy to compute or estimate because the partial derivatives are complicated multivariate polynomials. However, our CD algorithms do not require $L_{i}$. This quantity is just used for theoretical convergence analysis. The minimum and maximum of all the component-wise Lipschitz constants, respectively, are:

Employing similar steps of Lemma 2, we can prove that the full gradient $\nabla f(\bar{\boldsymbol{x}})$ is Lipschitz continuous:

with $L$ being the “full” Lipschitz constant. It is not difficult to show $L\leq\sum_{i}L_{i}$ and thus we further have $L\leq 2NL_{\max}$.

Theorem 1: The CCD, RCD, and GCD globally converge to a stationary point of the multivariate quartic polynomial from an arbitrary initialization.

Proof: Based on the component-wise Lipschitz continuous property of (37), it is derived that:

from which we obtain a lower bound on the progress made by each CD iteration

For different rules of index selection, the right-hand side of (47) will differ. We discuss the GCD, RCD, and CCD, one by one as follows. For GCD, it chooses the index with the largest partial derivative in magnitude. With the use of (12), we then have:

Substituting (48) into (47) leads to the following lower bound of the progress of one GCD iteration

This means that one GCD iteration decreases the objective function with an amount of at least $\frac{\|\nabla f(\bar{\boldsymbol{x}}^{k})\|^{2}}{4NL_{\max}}$. Setting $k=0,\cdots,j$, in (49) and summing over all inequalities yields

where we use $f(\bar{\boldsymbol{x}}^{j+1})\geq 0$. Taking the limit as $j\rightarrow\infty$ on (50), we get a convergent series

If a series converges, then its terms approach to zero, which indicates

i.e., the GCD converges to a stationary point.

For RCD, since $i_{k}$ is a random variable, $f(\bar{\boldsymbol{x}}^{k+1})$ is also random and we consider its expected value:

where the fact that $i_{k}$ is uniformly sampled from $\{1,\cdots,2N\}$ with equal probability of $1/(2N)$ is employed. Then the RCD at least obtains a reduction on the objective function in expectation

Following similar steps in the GCD, it is easy to prove that the expected gradient of the RCD approaches to the zero vector and thus it converges to a stationary point in expectation. For CCD, $i_{k}$ takes value cyclically from $\{1,\cdots,2N\}$. Applying Lemma 3.3 in for cyclic block CD, we can derive a lower bound of the decrease of the objective function after $2N$ CCD iterations:

Setting $k=0,2N,\cdots,2jN$, in (55) and summing over all the inequalities yields

Taking the limit as $j\rightarrow\infty$ of (56) yields a convergent series. Thus, the gradient approaches to zero, indicating that the CCD converges to a stationary point. $\Box$

We emphasize that the “global” convergence to a stationary point means that CD converges from an arbitrary initial value. Unlike local convergence, it does not require the initial value to be close enough to the stationary point.

Remark 1: Several existing convergence analyses of (block) CD, e.g., Proposition 2.7.1 of Bertsekas’ book and page 153 of , assume that the minimum of each block/coordinate is uniquely attained. However, our analysis in Theorem 1 does not require this assumption.

Remark 2: Theorem 6.1 of provides a convergence result for a descent method using update formula $\bar{\boldsymbol{x}}^{k+1}=\bar{\boldsymbol{x}}^{k}+\delta_{k}\boldsymbol{t}^{k}$, where $\boldsymbol{t}^{k}$ is a descent direction and $\delta_{k}>0$ is the stepsize. Theorem 6.1 of has proved

if $\delta_{k}>0$ is determined by an inexact line search procedure to ensure sufficient decrease at each iteration. Since the full gradient descent method adopts $\boldsymbol{t}^{k}=-\nabla f(\bar{\boldsymbol{x}}^{k})$, (57) becomes $\|f(\bar{\boldsymbol{x}}^{k})\|\rightarrow 0$ and hence $f(\bar{\boldsymbol{x}}^{k})\rightarrow\boldsymbol{0}$. Then Theorem 6.1 of proves that the full gradient descent converges to a stationary point. For CDs, it has $\boldsymbol{t}^{k}=\boldsymbol{e}_{i_{k}}$ and (57) becomes $\nabla_{i_{k}}f(\bar{\boldsymbol{x}}^{k})\rightarrow 0$. Clearly, we can only conclude a single partial derivative approaches zero and cannot conclude other partial derivatives approach zero. Therefore, Theorem 6.1 of cannot be used to prove the convergence to a stationary point for CDs. Moreover, the proof of Theorem 6.1 of requires the gradient is globally Lipschitz continuous, which results in that it is not applicable to our problem. In addition, uses an inexact line search for stepsize while the CDs adopt exact coordinate minimization. The self-contained convergence analysis of CDs is totally different from .

Remark 3: Even when there are enough samples, the Hessian matrix $\nabla^{2}f(\bar{\boldsymbol{x}})$ close to the minimizer $\bar{\boldsymbol{x}}^{\star}$ has $2N-1$ positive eigenvalues, and the remaining eigenvalue can be zero, positive, or negative. This implies that $f(\bar{\boldsymbol{x}})$ can never be locally convex no matter how small the local region around $\bar{\boldsymbol{x}}^{\star}$ is. Therefore, the established results for convergence rate using convexity are not applicable for our nonconvex problem.

III-B Local Convergence to Global Minimum

Theorem 1 just shows that the CD algorithm converges to a stationary point. A further question is: can the CD converge to the global minimizer and hence exactly recovers the original signal? At first glance, it seems impossible because even finding a local minimum of a fourth-order polynomial is known to be NP-hard in general . However, the answer is yes under the condition that the sample size is large enough. The backbone of the proof is based on a statistical analysis of the gradient of the nonconvex objective function established by Candès et al. . It is worth mentioning that the convergence analysis of WF is for the complex-valued full gradient method and cannot be directly applied to our real-valued problem using coordinate minimization.

Recall that if $\boldsymbol{x}^{\star}$ is an optimal solution of (2), then all the elements of the following set

and the minimum of (59) attains at $\phi=\phi(\boldsymbol{z})$. Similarly, the set of all optimal solutions of the real-valued problem (8) is defined as

where $\bar{\boldsymbol{x}}^{\star}=[{\rm Re}(\boldsymbol{x}^{\star})^{T},{\rm Im}(\boldsymbol{x}^{\star})^{T}]^{T}$ is a global minimizer of (8). That is, $T_{\phi}(\bar{\boldsymbol{x}}^{\star})$ denotes the effect of a phase rotation to $\bar{\boldsymbol{x}}^{\star}$. The projection of $\bar{\boldsymbol{x}}^{k}$ onto $\mathcal{P}$ is the point in $\mathcal{P}$ closest to $\bar{\boldsymbol{x}}^{k}$, which is denoted as $T_{\phi_{k}}(\bar{\boldsymbol{x}}^{\star})$ where

Then the distance of $\bar{\boldsymbol{x}}^{k}$ to $\mathcal{P}$ is

Our goal is to prove ${\rm dist}(\bar{\boldsymbol{x}}^{k},\mathcal{P})\rightarrow 0$. The following lemma of , which essentially states that the gradient of the objective function is well behaved, is crucial to our proof.

where $\rho>0$ and $\eta>0$, holds with high probability if the number of measurements satisfies $M\geq C_{0}N\log N$ with $C_{0}>0$ being a sufficiently large constant.

The detailed proof of Lemma 3 can be found in Condition 7.9, Theorem 3.3, and Sections 7.5–7.7 of .

Although the regularity condition of Lemma 3 corresponds to the complex-valued case, we at once obtain the real-valued version according to (16), (60), and (61). For $\bar{\boldsymbol{x}}^{k}$ satisfying ${\rm dist}(\bar{\boldsymbol{x}}^{k},\mathcal{P})\leq\epsilon$, we have

Theorem 2: Assume that the sample size satisfies $M\geq C_{0}N\log N$ with a sufficiently large $C_{0}$ and ${\rm dist}(\bar{\boldsymbol{x}}^{0},\mathcal{P})\leq\epsilon$. The iterates of the RCD with a slight modification, in which the one-dimensional search is limited to a line segment, i.e.,

satisfy ${\rm dist}(\bar{\boldsymbol{x}}^{k},\mathcal{P})\leq\epsilon$ for all $k$ and converge to $\mathcal{P}$ in expectation with high probability at a geometric rateThe geometric convergence rate is also called linear convergence rate in the optimization literature. It indicates that the logarithm of the error decreases linearly.

Proof: The updating equation of the CD, i.e., $\bar{\boldsymbol{x}}^{k+1}=\bar{\boldsymbol{x}}^{k}+\alpha_{k}\boldsymbol{e}_{i_{k}}$, is equivalently expressed as

where $\gamma_{k}=-\alpha_{k}/\nabla f_{i_{k}}(\bar{\boldsymbol{x}}^{k})$. It requires $\gamma_{k}>0$ to ensure $f(\bar{\boldsymbol{x}}^{k+1})<f(\bar{\boldsymbol{x}}^{k})$. Hence, $|\alpha_{k}|\leq 2\eta|\nabla f_{i}(\bar{\boldsymbol{x}}^{k})|$ means $0<\gamma_{k}\leq 2\eta$. Employing the development starting from (62), it follows

where the last line follows from (64). Combining (68) and (69) yields

where the last inequality follows from $0<\gamma_{k}\leq 2\eta$. Successively applying (70), we get

where $\gamma_{\min}=\min\limits_{1\leq j\leq k}\gamma_{j}$. $\Box$

Remark 4: To guarantee convergence to the globally optimal solution, it requires $|\alpha|\leq 2\eta|\nabla f_{i_{k}}(\bar{\boldsymbol{x}}^{k})|$ or equivalently $0<\gamma_{k}\leq 2\eta$. If $\eta$ is known or can be estimated, we can perform the one-dimensional search of (65) limited to a line segment. Note that (65) is on minimizing a univariate quartic polynomial in an interval. This problem is easy to solve because its solution belongs to the stationary points in the interval (if there indeed exists such a stationary point in the interval) or the endpoints of the interval. However, $\eta$ is always not easy to estimate in practice. From simulations, we find that dropping the box constraint $0<\gamma_{k}\leq 2\eta$ will not destroy the convergence. This implies that the box constraint is automatically satisfied. We conjecture $\eta$ is large enough such that $\gamma_{k}\leq 2\eta$ is always guaranteed when there are enough samples. Therefore, this empirical observation ensures us to ignore the constraint $\gamma_{k}\leq 2\eta$ at each coordinate minimization.

Remark 5: We only prove convergence to the global minimizer for RCD. For CCD and GCD, theoretical proof of the convergence remains open and constitutes a future research. Nonetheless, it is observed from the numerical simulations that the GCD converges faster than the RCD, and CCD has comparable performance to RCD. Therefore, empirically, the GCD and CCD also converge to the global minimum point with high probability if the sample size is large enough.

By ignoring the terms independent to $\alpha$, (73) is equivalent to

Making a change of variable $\beta=\alpha+\bar{x}_{i}$, substituting $\alpha=\beta-\bar{x}_{i}$ into (27), and ignoring the constant term, we obtain an equivalent scalar minimization problem

where the coefficients of the quartic polynomial $\{u_{j}\}_{j=1}^{4}$ are calculated as

It is interesting that the solution of (75) reduces to the well-known soft-thresholding operator in compressed sensing if $u_{4}=u_{3}=0$, where the quartic polynomial reduces to a quadratic function. Therefore, (75) is a generalization of the soft-thresholding operator from quadratic to fourth-order functions. We call it fourth-order soft-thresholding (FOST). Although $\psi(\beta)$ is non-smooth due to the absolute term, the closed-form solution of its minimum can still be derived. We study the minimizer of $\psi(\beta)$ in two intervals, namely, $[0,\infty)$ and $(-\infty,0)$. Define the set $\mathcal{S}^{+}$ containing the stationary points of $\psi(\beta)$ in the interval $[0,\infty)$. That is, $\mathcal{S}^{+}$ is the set of real positive roots of the cubic equation

The $\mathcal{S}^{+}$ can be empty, or has one or three elements because (77) may have none, one, or three real positive roots. Similarly, $\mathcal{S}^{-}$ is the set that contains the stationary points of $\psi(\beta)$ in $(-\infty,0)$, i.e., real negative roots of

Again, $\mathcal{S}^{-}$ can be empty, or has one or three entries. The minimizer of $\psi(\beta)$ in $\beta\in[0,\infty)$ must be the boundary, i.e., 0, or one element of $\mathcal{S}^{+}$. The minimizer in $(-\infty,0)$ must be an element of $\mathcal{S}^{-}$. In summary, the minimizer of (75) is limited to the set $\{0\cup\mathcal{S}^{+}\cup\mathcal{S}^{-}\}$ which has at most seven elements, i.e.,

V Application to Blind Equalization

We illustrate the application of phase retrieval to blind equalization, which is a fundamental problem in digital communications. Consider a communication system with discrete-time complex baseband signal model

where $r(n)$ is the received signal, $s(n)$ is the transmitted data symbol, $h(n)$ is the channel impulse response, $\nu(n)$ is the additive white noise, and $*$ denotes convolution. The received signal is distorted due to the inter-symbol interference (ISI) induced by the propagation channel. Channel equalization is such a technique to mitigate the ISI. Blind equalization aims at recovering the transmitted symbols without knowing the channel response. Define the equalizer with $P$ coefficients $\boldsymbol{w}=[w_{0},\cdots,w_{P-1}]^{T}$ and $\boldsymbol{r}_{n}=[r(n),\cdots,r(n-P+1)]^{T}$, the equalizer output is

As many modulated signals in communications such as phase shift keying (PSK), frequency modulation (FM), and phase modulation (PM), are of constant modulus (CM), we apply the CM criterion to obtain the equalizer:

where $\kappa>0$ is the dispersion constant defined as :

If $s(n)$ is of strictly constant modulus, e.g., for PSK signals, then $\kappa$ equals the square of modulus. It is obvious that the problem of CM based blind equalization in (82) has the same form as the phase retrieval of (2). Both of them are multivariate quartic polynomials. The only difference between phase retrieval and blind equalization is that the decision variable of the former is the unknown signal $\boldsymbol{x}$ while that of the latter is the equalizer $\boldsymbol{w}$. Therefore, the WF and CD methods can be applied to solve (82). By defining the composite channel-equalizer response as $v(n)=h(n)*w(n)$, the quantified ISI, which is expressed as

reflects the equalization quality. Smaller ISI implies better equalization. If ISI $=0$, then the channel is perfectly equalized and the transmitted signal is exactly recovered up to a delay and a scalar. Perfect equalization is only possible when there is no noise and the equalizer length $P$ is infinite for finite impulse response (FIR) channelThe equalizer is the inverse system of the channel. If the channel is of FIR, then its inverse has infinite impulse response (IIR). Hence, an equalizer with infinite length is required for perfectly equalizing an FIR channel.. Otherwise, only approximate equalization can be achieved, which results in a residual ISI.

VI Simulation Results

In our simulation study, all methods use the same initial value obtained from the spectral method . The sampling vectors $\{\boldsymbol{a}_{m}\}$ satisfy a complex standard i.i.d. Gaussian distribution.

We first investigate the convergence behavior of the three CD algorithms. The signal $\boldsymbol{x}$ and noise $\nu_{m}$ are i.i.d. Gaussian distributed. In this test, we set $N=64$ and $M=6N$. The WF and WFOS that uses optimal stepsize for accelerating the convergence speed of WF, are employed for comparison. Note that it is fair to compare $2N$ iterations (one cycle) for the CD with one WF or WFOS iteration because the computational complexity of the CCD and RCD per cycle is the same as the WF per iteration. The GCD has a higher complexity for every $2N$ iterations than WF, CCD, and RCD. But still, we plot the results of GCD per cycle. Two quantities are plotted to evaluate the convergence rate. The first quantity is the reduction of the objective function normalized with respect to $\|\boldsymbol{b}\|^{2}$:

where $f(\bar{\boldsymbol{x}}^{\star})=0$ if there is no noise. For the noisy case, $f(\bar{\boldsymbol{x}}^{\star})$ can be computed in advance to the machine accuracy using the CD or WF method. The second quantity is the relative recovery error, i.e.,

which reflects the convergence speed to the original signal. Fig. 1 plots the objective reduction while Fig. 2 shows the recovery error, versus the number of iterations (cycles for CD) in the absence of noise with 50 independent trials. The averaged results are also provided with thick lines. We see that all methods converge to the global minimum point at a linear rate. They exactly recover the true signal. For the noisy case, the signal-to-noise ratio (SNR) in (1) is defined as

where $\sigma_{\nu}^{2}$ is variance of $\nu_{m}$. Figs. 3 and 4 show the normalized objective reduction and recovery error, respectively, at SNR $=$ 20 dB. We clearly see that the three CD algorithms converge faster than the WF and WFOS schemes. Among them, the convergence speed of the GCD is the fastest.

VI-B Statistical Performance

The experiment settings are the same as in Section VI-A except $M$ and SNR vary. The performance of the GS algorithm is also examined here. We use the empirical probability of success and normalized mean square error (NMSE), which is the mean of the relative recovery error in (86), to measure the statistical performance. All results are averaged over 200 independent trials. In the absence of noise, if the relative recovery error of a phase retrieval scheme is smaller than $10^{-5}$, we call it success in exact recovery. Fig. 5 plots the empirical probability of success versus number of measurements $M$. It is observed that the GCD is slightly better than WF while CCD and RCD are slightly inferior to the WF. Fig. 6 shows the NMSE versus SNR from 6 dB to 30 dB. We see that the three CD algorithms and WF have comparable NMSEs and they are superior to the GS algorithm.

VI-C Phase Retrieval of Sparse Signal

VI-D Blind Equalization

We investigate the application of the CD and WF methods to blind equalization in the presence of white Gaussian noise. The results of the super-exponential (SE) algorithm are also included. The transmitted signal adopts quadrature PSK (QPSK) modulation, namely, $s(n)\in\{1,-1,{\rm j},-{\rm j}\}$. A typical FIR communication channel with impulse response $\{0.4,1,-0.7,0.6,0.3,-0.4,0.1\}$ is adopted . Fig. 9 shows the constellations of the received signal and equalizer outputs of 1000 samples at SNR = 20 dB. We observe that the received signal is severely distorted due to the channel propagation. The SE, WF and CD methods succeed in recovering the transmitted signal up to a global phase rotation. We clearly see that the CCD and GCD have a higher recovery accuracy. Fig. 10 plots the ISI versus the number of iterations/cycles at SNR = 25 dB with 2000 samples. The ISI is averaged over 100 independent trials. In addition to faster convergence than the WF, WFOS, and SE, the CDs (especially CCD) arrive at a lower ISI. This means that the CDs also achieve a more accurate recovery.

VII Conclusion

Acknowledgment

W.-J. Zeng would like to thank Prof. Emmanuel J. Candès at Stanford University for his suggestions and encouragement. He is also grateful to Prof. Xiaodong Li at University of California at Davis for his explanations on some questions of the WF method for phase retrieval.

References