Phaseless Rcovery using Gauss-Newton Method

I Introduction

where $x_{k}$ is the output of the $k$-th iteration of WF method, $z$ is the true signal, $0<\rho<1$ is a constant and the definition of ${\rm dist}(\cdot)$ is given in Section I-C. The truncated WF method is introduced in , which improves the performance of WF method with showing that $O(n)$ Gaussian random measurements are enough to attain the linear convergence rate. Recently another two stage iterative algorithm, the truncated amplitude flow (TAF), was proposed by Wang, Giannakis and Eldar in . The TAF uses null initialization method to obtain an initial estimate and refines it by successive updates of truncated generalized gradient iterations. It is proved that TAF can geometrically converge to the exact signal with $O(n)$ measurements . Despite iterative algorithms to solve phaseless recovery problem, a recent approach is to recast phaseless recovery as a semi-definite programming (SDP), such as PhaseLift . PhaseLift is to lift a vector problem to a rank-1 matrix one and then one can recover the rank-1 matrix by minimizing the trace of matrices. Though PhaseLift can provide the exact solution using $O(n)$ measurements, the computational cost is large when the dimension of the signal is high.

I-B Our contribution

The aim of this paper is twofold. We first present an alternative initial guess which is the eigenvector corresponding to the largest eigenvalue of

where $x_{k}$ is the output of the $k$-th iteration and $\beta$ is a constant. Hence, to reach the accuracy $\epsilon$, re-sampled Gauss-Newton method needs $O(\log\log\frac{1}{\epsilon})$ iterations, which has an improvement over the Altmin Phase algorithm. Since many signals from real world are real, the assumption of $z$ being real is reasonable. For the case where signals are complex, we derive a revised Gauss-Newton method.

I-C Notations

As the problem setup naturally leads to ambiguous solutions, we define

I-D Organization

II Initialization

For non-convex problem (I.1), proper initial criteria is essential to avoid the iterative algorithm trapping in a local minimum. So the first step of Gauss-Newton method is to choose an initial estimation. Before giving our initialization method, we first review several other methods.

Spectral initialization method estimates the initial guess $z_{0}$ as the eigenvector corresponding to the largest eigenvalue of $\frac{1}{m}\sum_{j=1}^{m}y_{j}a_{j}a_{j}^{*}$ with norm $\sqrt{\sum_{j=1}^{m}y_{j}/m}$. In , Candès, Li and Soltanolkotabi prove that when $a_{j},\,j=1,\ldots,m$ are Gaussian random measurements with $m\geq C_{0}n\log n$, ${\rm dist}(z_{0},z)\leq 1/8\|z\|$ holds with probability at least $1-10\exp(-\gamma n)-8/n^{2}$. To reduce the number of observations, a modified spectral method is introduced in , which precludes $y_{j}$ with large magnitudes. Particularly, they select the initial value as the eigenvector corresponding to the largest eigenvalue of $\frac{1}{m}\sum_{j=1}^{m}y_{j}a_{j}a_{j}^{*}I_{\{|y_{j}|\leq\beta_{y}\lambda^{2}\}}$, where $\beta_{y}$ is an appropriate truncation criteria and $\lambda^{2}=\sum_{j}y_{j}/m$. This method only requires the number of measurements $m\geq Cn$ with a sufficient large constant $C$. The null initialization method is introduced by Chen, Fannjiang and Liu in . This method builds on the orthogonality characteristics of high-dimensional random vectors, and choose the eigenvector corresponding to the largest eigenvalue of $\frac{1}{|I|}\sum_{j\in I}a_{j}a_{j}^{*}$ as the initial guess, where $I$ is an index set selected by the $|I|$ largest magnitudes of $\frac{|\langle a_{j},z\rangle|^{2}}{\|a_{j}\|^{2}\|z\|^{2}}$, $j=1,\ldots,m$. When the number of measurements is on the order of $n$, null initialization method can guarantee a good precision. More details can be found in , . To state conveniently, we name the first method as SI (Spectral Initialization), the second method as TSI (Truncated Spectral Initialization) and the third method as NI (Null Initialization).

Next we introduce a new method for initialization, which is stated in Algorithm 1. In fact, the initial guess is chosen as the eigenvector corresponding to the largest eigenvalue of the Hermitian matrix

Noting that $\lambda^{2}$ is a good approximation to $\|z\|^{2}$. So we choose $Y$ as an approximation to $\dfrac{zz^{*}}{4\|z\|^{2}}$, whose eigenvector of the largest eigenvalue is of the form $cz$ where $c$ is a constant. Meanwhile, the eigenvector associated with the largest eigenvalue of $Y$ can be efficiently calculated by the power method (see details in ).

The new method makes full use of every observation and can obtain an alternative initial value by nearly optimal number of measurements (see Theorem II.1). Beyond theoretical results, numerical experiments also show that this method has better performance than that of SI, TSI and NI (see Example IV.1).

II-B The performance of Algorithm 1

holds with probability at least $1-4\exp(-c_{\theta}n)$, where $c_{\theta}>0$.

where $\lambda^{2}=\dfrac{1}{m}\sum\limits_{j=1}^{m}y_{j}$. Then for any $\theta>0$, we have

with probability at least $1-4\exp(-c_{\theta}n)$ provided $m\geq C_{\theta}n$. Using similar method with the proof of Theorem II.1, we can obtain

It is possible to obtain similar results with replacing $\exp(-y_{j}/\lambda^{2})$ in $Y$ by another bounded function $g(y_{j})$. For example, we can take $g(y_{j})=\exp(-y_{j}^{p}/\lambda^{2})$ where $0<p\leq 1$. We need adjust the constant $1/2$ in $Y$ when we replace the function $\exp(-y_{j}/\lambda^{2})$ in $Y$ by another bounded function $g(y_{j})$.

III Gauss-Newton Method

In this section, we present Gauss-Newton iterations which are used to refine the initial guess.

where $y_{j}=\lvert\langle a_{j},z\rangle\rvert^{2}$. To state conveniently, we set $F_{j}(x):=\frac{1}{\sqrt{m}}(\langle a_{jR},x\rangle^{2}+\langle a_{jI},x\rangle^{2}-y_{j})$ and we write (III.3) in the form of

To solve the nonlinear least square problem (III.4), our algorithm uses the well-known Gauss-Newton iteration. To make the paper self-contained, we introduce the Gauss-Newton iteration in detail (see also ). Suppose the $k$-th iteration point $x_{k}$ is real-valued, we first linearize the nonlinear term $F_{j}(x)$ at the point $x_{k}$:

We choose the next iteration point $x_{k+1}$ as the solution to (III.5), i.e.,

Thus we obtain the update rule (III.6). Note that the $x_{k+1}$ is also real-valued.

III-A2 Gauss-Newton Method with Re-sampling

The Gauss-Newton method uses Algorithm 1 to obtain an initial guess $x_{0}$ and iteratively refine $x_{k}$ by the update rule (III.6). In theoretical analysis, as we require that the current measurements are independent with the last iteration point (see Ramark III.2), we re-sample measurement matrix $A$ in every iteration step. Then Algorithm 2 is in fact a variant of Gauss-Newton method with using different measurements in each iteration. The re-sampling idea is also used in for the alternating minimization algorithm and in for the WF algorithm with coded diffraction patterns.

In step 4 of the Algorithm 2, the concrete form of matrix $J^{k+1}{(x_{k})}^{\top}J^{k+1}(x_{k})$ and vector $J^{k+1}(x_{k})^{\top}F^{k+1}(x_{k})$ is same with (III-A1) and (III-A1). Here $y_{1},\ldots,y_{m^{\prime}}$ are the entries of $y^{(k+1)}$ and $a_{1},\ldots,a_{m^{\prime}}$ are the rows of $A^{(k+1)}$.

III-A3 Convergence Property of Gauss-Newton Method with Re-sampling

We next present theoretical convergence property of Algorithm 2. Without loss of generality, we assume $\|z\|=1$. Theorem III.1 illustrates that under given conditions, Algorithm 2 has a quadratic convergence rate. Furthermore, we show that to achieve an $\epsilon$ accuracy, the Gauss-Newton method with re-sampling only needs $O(\log\log(\frac{1}{\epsilon}))$ iterations.

In Theorem III.1, the reason why we require $0<\delta\leq 1/93$ is to guarantee $\beta\cdot\delta\leq\sqrt{\delta}$. Hence the condition $\textup{dist}(x_{k+1},z)\leq\beta\cdot\delta\leq\sqrt{\delta}$ still holds and we can use Theorem III.1 at the $(k+1)$-th iteration.

According to Theorem II.1 or Remark II.1, for any $0<\delta\leq 1/93$ and $0<\theta\leq\delta/3$, when $m\geq C_{\theta}n$, it holds with probability at least $1-4\exp(-c_{\theta}n)$ that

Combining this initialization result with Theorem III.1, we have the following conclusion.

where $C$ is a constant depending on $\delta$, $\epsilon$.

According to Theorem II.1 or Remark II.1, we have

with probability at least $1-4\exp(-c_{\delta}n)$. From Remark III.1, we know

where $\beta$ is defined in Theorem III.1. In Algorithm 2, we choose $T=c\log\log\frac{1}{\epsilon}$ and $m^{\prime}\geq C_{1}n\log n$, where $C_{1}$ is a constant depending on $C,\,c$. Iterating (III.9) in Theorem III.1 $T$ times leads to

In Algorithm 2, we use different measurement vectors in each iteration. In fact, Theorem III.1 requires that the Gaussian random measurement vectors $a_{j}$ are independent with $x_{k}$. According to (III.6), $x_{k+1}$ depends on the current measurement vectors $a_{j}$. Hence, to use Theorem III.1 at the next step, we need choose different measurement vectors which are independent with the previous ones.

III-B Complex-valued signals

Using a similar argument with above, at the $k$-th iteration, we can update $x_{k}$ by solving

Noting that $A_{k}\begin{pmatrix}x_{k}\\ -\overline{x_{k}}\end{pmatrix}=0,$ we have

which implies the conclusion. Here, we use the definition of $x^{\sharp}$ (see (III.12)). ∎

We denote the solution set to (III.13) as ${\mathcal{L}}_{k}$. Our idea is to choose $\hat{x}\in{\mathcal{L}}_{k}$ so that $\|{x}_{k+1}-x_{k}\|_{2}=\|\hat{x}\|_{2}$ reaches the minimum since we already know $x_{k}$ is not far from the exact signal. Then we have

We use $A_{k}^{{\dagger}}$ to denote the moore-penrose pseudoinverse of $A_{k}$. Then

which implies (III.15) since $A_{k}A_{k}^{*}+\delta I$ is a real matrix. ∎

The numerical experiments show (III.16) has quadratic convergence rate provided the initial guess is not far from the exact signal (see Example IV.2 (b)). The analysis of the convergence property of (III.16) is the subject of our future work.

IV Numerical Experiments

V Appendix

To prove the Theorem II.1, we first recall some useful results.

holds with probability at least $1-2\exp(-c_{\zeta}m)$. Here $C_{\zeta}$, $c_{\zeta}$ depend on the constant $\zeta$ and the sub-gaussian norm $\max_{j}\|a_{j}\|_{\psi_{2}}$.

The next lemma plays an essential role in proving Theorem II.1.

holds with probability at least $1-4\exp(-c_{\eta}n)$ provided $m\geq C_{\eta}n$, where $c_{\eta}>0$, $C_{\eta}$ are constants depending on $\eta$.

holds with probability at least $1-4\exp(-c_{\eta}n)$ provided $m\geq C_{\eta}n$, where $c_{\eta}$, $C_{\eta}$ are constants depending on $\eta$ and subgassian norm of $a_{j}$, $\sqrt{\exp(-|a_{j}^{*}z|^{2}/\|z\|^{2})}a_{j}$, $j=1,\ldots,m$. The inequality (V.19) also implies that

holds with high probability. Combining (V.19) and (V.20), we obtain that

where the second inequality dues to (V.21) and the fact that $x\exp(-x)<1$ for all $x$. The second line of (V.23) uses the Lagrange’s mean value theorem with $\xi\in[(1-\frac{\eta}{4})\|z\|^{2},\,\,(1+\frac{\eta}{4})\|z\|^{2}]$ with high probability. Thus putting (V.22) and (V.23) into (V.18), we get

So the conclusion holds with probability at least $1-4\exp(-c_{\eta}n)$ provided $m\geq C_{\eta}n$, where $c_{\eta}$, $C_{\eta}$ are constants depending on $\eta$. ∎

From Lemma V.2, for any $0<\theta\leq 1$ and $m\geq C_{\theta}n$, we have

with probability at least $1-4\exp(-c_{\theta}n)$. Note that the largest eigenvalue of $\dfrac{zz^{*}}{4\|z\|^{2}}$ is $\dfrac{1}{4}$. Then according to the Wely Theorem,

holds with probability at least $1-4\exp(-c_{\theta}n)$. On the other hand,

From the proof of Lemma V.2 (see (V.21)), we have

with probability at least $1-4\exp(-c_{\theta}n)$ provided $m\geq C_{\theta}n$. Thus we get the conclusion

V-B Proof of Theorem III.1

In this section, we devote to prove the Theorem III.1. At first, we give some essential lemmas.

holds with probability at least $1-5\exp(-\gamma_{\delta}n)-4/n^{2}$.

Recall that $S_{k}=\{tz+(1-t)x_{k}:0\leq t\leq 1\}$. We set

is $L_{J}$-Lipschitz continuous on $S_{k}$ with probability at least $1-5\exp(-\gamma_{\delta}n)-4/n^{2}$, i.e, for any $x,y\in S_{k}$,

holds with $L_{J}=8(2+\frac{\delta}{4})(1+\sqrt{\delta})$.

Since the measurement vectors $a_{j}$, $j=1,\ldots,m$ are rotationally invariant and independent with $x_{k}$ and $z$, wlog, we can assume that $z=e_{1}$ and $x_{k}=\|x_{k}\|(\alpha e_{1}+\sqrt{1-\alpha^{2}}e_{2})$, where $\alpha=\langle x_{k},z\rangle/\|x_{k}\|$. As $\|x_{k}-z\|\leq\sqrt{\delta}$, so $\langle x_{k},z\rangle\geq 0$, i.e., $\alpha\geq 0$. We can write $x,y\in S_{k}$ in the form of

where $\sigma_{R,+}:=a_{jR}^{\top}(x+y)$, $\sigma_{R,-}:=a_{jR}^{\top}(x-y)$, $\sigma_{I,+}:=a_{jI}^{\top}(x+y)$, $\sigma_{I,-}:=a_{jI}^{\top}(x-y)$, $\kappa_{1}:=\sqrt{(a_{jR}^{\top}e_{1})^{2}+(a_{jR}^{\top}e_{2})^{2}}$ and $\kappa_{2}:=\sqrt{(a_{jI}^{\top}e_{1})^{2}+(a_{jI}^{\top}e_{2})^{2}}$ and the last inequality is obtained by Cauchy-Schwarz inequality. Next we set

holds with probability at least $1-5\exp(-\gamma_{\delta}n)-4/n^{2}$. So

On the other hand, as $\|x_{k}-z\|\leq\sqrt{\delta}$, we have

Putting (V.28) and (V.30) into (V.27), we obtain

So we conclude that when $m\geq Cn\log n$, $J(x)^{\top}J(x)$ is Lipschitz continuous on the line $S_{k}$ with constant $L_{J}=8(2+\frac{\delta}{4})(1+\sqrt{\delta})$ with probability at least $1-5\exp(-\gamma_{\delta}n)-4/n^{2}$. ∎

Under the same conditions as in Lemma V.4,

is Lipschitz continuous on $S_{k}$ with Lipschitz constant

with probability at least $1-5\exp(-\gamma_{\delta}n)-4/n^{2}$.

So according to Lemma V.3, for $0<\delta\leq 1/93$ and $m\geq Cn\log n$,

holds with probability at least $1-5\exp(-\gamma_{\delta}n)-4/n^{2}$. So we have

Putting (V.32) and (V.30) into (V.31), we have

So $H(x)$ is Lipschitz continuous on $S_{k}$ with constant $L_{H}=4(1+\sqrt{\delta})(3+\frac{\delta}{4})$. ∎

Next we present an estimation of the largest eigenvalue of $(J(x_{k})^{\top}J(x_{k}))^{-1}$.

According to Lemma V.3, for $0<\delta\leq 1/93$ and $m\geq Cn\log n$,

holds with probability at least $1-5\exp(-\gamma_{\delta}n)-4/n^{2}$. Then according to the Wely Theorem, we have

Here, we use (V.29) in the last inequality. Then with probability at least $1-5\exp(-\gamma_{\delta}n)-4/n^{2}$, we have

We next present the proof of Theorem III.1.

Without loss of generality, we suppose $\langle x_{k},z\rangle\geq 0$, i.e.,

Then we just need to prove when $\|x_{k}-z\|\leq\sqrt{\delta}$ and $m\geq Cn\log n$,

holds with probability at least $1-c/n^{2}$.

As $z$ is an exact solution to (III.3), we have $\nabla f(z)=H(z)=0$. The definition of $x_{k+1}$ shows that

Define $S_{k}:=\{x_{k}+t(z-x_{k}):0\leq t\leq 1\}$ and $x(t)=x_{k}+t(z-x_{k})$. Then we have

The integral in (V.35) is interpreted as element-wise. Combining (V.26) and $H(z)=0$, we obtain

According to Lemma V.4 and Corollary V.1, $J(x)^{\top}J(x)$ and $H(x)$ are Lipschitz continuous on the line $S_{k}$ with probability at least $1-5\exp(-\gamma_{\delta}n)-4/n^{2}$ provided $m\geq Cn\log n$. So using (V.36), we obtain

Thus according to Lemma V.5 and (V.34), when $m\geq Cn\log n$,

holds with probability at least $1-c/n^{2}$. Based on the discussion in Theorem II.1, we have

Then we have $\langle x_{k+1},z\rangle\geq 0$, i.e., $\text{dist}(x_{k+1},z)=\|x_{k+1}-z\|$. ∎

Acknowledgements. We are grateful to Xin Liu for discussions and comments at the beginning of this project, which contributed to the proof of Theorem III.1. We would like to thank the referees for thorough and useful comments which have helped to improve the presentation of the paper.

References