A New Inexact Proximal Linear Algorithm with Adaptive Stopping Criteria for Robust Phase Retrieval

Phase retrieval aims to recover a signal from intensity-based or magnitude-based measurements. It finds various applications in different fields, including X-ray crystallography [1], optics [2], array and high-power coherent diffractive imaging [3], astronomy [4] and microscopy [5]. Mathematically, phase retrieval tries to find the true signal vectors $x_{\star}$ or $-x_{\star}$ in $\mathbb{R}^{n}$ from a set of magnitude measurements:

$$b_{i}=(a_{i}^{\top}x_{\star})^{2},\text{ for }i=1,2,\ldots,m,$$

(1)

where $a_{i}\in\mathbb{R}^{n}$ and $b_{i}\geq 0,i=1,2,\ldots,m$. Directly solving the equations leads to an NP-hard problem [6], and nonconvex algorithms based on different designs of objective functions have been well studied in the literature, including Wirtinger flow [7], truncated Wirtinger flow [8], truncated amplitude flow [9], reshaped Wirtinger flow [10], etc.

In this paper, we focus on the robust phase retrieval (RPR) problem [11], which considers the case where $b_{i}$ contains noise due to measurement errors of equipment. That is,

$$b_{i}=\begin{cases}(a_{i}^{\top}x_{\star})^{2},&i\in\mathbb{I}_{1},\\ \xi_{i},&i\in\mathbb{I}_{2},\end{cases}$$

(2)

in which $\mathbb{I}_{1}\bigcup\mathbb{I}_{2}=\{1,2\ldots,m\}$, $\mathbb{I}_{1}\cap\mathbb{I}_{2}=\emptyset$, and $\xi_{i}$ denotes a random noise. [11] proposed to formulate RPR as the following optimization problem:

$$\min_{x\in\mathbb{R}^{n}}F(x):=\frac{1}{m}\sum_{i=1}^{m}\left|(a_{i}^{\top}x)^{2}-b_{i}\right|.$$

(3)

It is demonstrated in [11] that using (3) for RPR possesses better recoverability compared to the median truncated Wirtinger flow algorithm [12] based on the $\ell_{2}$-loss.

Solving (3) is challenging because it is a nonconvex and nonsmooth optimization problem. In [13], the authors proposed the subgradient method to solve it. This method requires geometrically decaying step size, and it is unclear how to schedule this kind of step size in practice. [11] proposed to use the proximal linear (PL) algorithm to solve (3). For ease of presentation, we rewrite (3) as

$$\min_{x\in\mathbb{R}^{n}}F(x)=h(c(x))=\frac{1}{m}\left\||Ax|^{2}-b\right\|_{1},$$

(4)

where $Ax=[\langle a_{1},x\rangle,\ldots,\langle a_{m},x\rangle]^{\top}$, $b=[b_{1},\ldots,b_{m}]^{\top}$, $h(z):=\frac{1}{m}\|z\|_{1}$, and $c(x):=|Ax|^{2}-b$ is a smooth map in which $|\cdot|^{2}$ is element-wise square. One typical iteration of the PL algorithm is

$$x^{k+1}\approx\operatorname*{arg\,min}_{x\in\mathbb{R}^{n}}\ F_{t}(x;x^{k}),$$

(5)

where $t>0$ is the step size,

$$F(z;y):=h(c(y)+\nabla c(y)(z-y)),$$

(6)

$$F_{t}(z;y):=F(z;y)+\frac{1}{2t}\|z-y\|_{2}^{2},$$

(7)

$\nabla c$ denotes the Jacobian of $c$, and “$\approx$” means that the subproblem is solved inexactly. The subproblem (5) is convex and can be solved by various methods such as the proximal operator graph splitting (POGS) algorithm used in [11]. The PL method has drawn lots of attention recently. It has been studied by [14, 15, 16] and applied to solving many important applications such as RPR [11], robust matrix recovery [17, 18], and sparse spectral clustering [19]. The subproblem (5) is usually solved inexactly for practical concerns. As pointed out in [11], the PL implemented in [11] is much slower than the median truncated Wirtinger flow method. We found that this is mainly due to their stopping criterion for solving the subproblem (5), which unnecessarily solves (5) to very high accuracy in the early stage of the algorithm. Moreover, we found that the POGS algorithm used in [11] is ineffective in solving the subproblem (5). In this paper, we propose adaptive stopping criteria for inexactly solving (5) with the fast iterative shrinkage-thresholding algorithm (FISTA) [20, 21, 22]. We found that our new inexact PL (IPL) with the adaptive stopping criteria greatly outperforms existing implementations of PL methods [11] for solving RPR (4).

Our Contributions. In this paper, we propose two new adaptive stopping criteria for inexactly solving (5). The first one ensures that (5) is solved to a relatively low accuracy:

	(LACC)	$\displaystyle\quad F_{t}(x^{k+1};x^{k})-\min_{x\in\mathbb{R}^{n}}\ F_{t}(x;x^{k})$		(8)
		$\displaystyle\leq\rho_{l}\left(F_{t}(x^{k};x^{k})-F_{t}(x^{k+1};x^{k})\right),\quad\rho_{l}>0,$

and the second one ensures that (5) is solved to a relatively high accuracy:

	(HACC)	$\displaystyle\quad F_{t}(x^{k+1};x^{k})-\min_{x\in\mathbb{R}^{n}}\ F_{t}(x;x^{k})$		(9)
		$\displaystyle\leq\frac{\rho_{h}}{2t}\|x^{k+1}-x^{k}\|_{2}^{2},\quad 0<\rho_{h}<1/4.$

Here, $\rho_{l}$ and $\rho_{h}$ are given constants. Similar to the proximal bundle method [23] for nonsmooth convex problems, (LACC) and (HACC) are designed to ensure the sufficient decrease of the objective function for the nonsmooth and nonconvex RPR problem. Note that both (LACC) and (HACC) are only used theoretically because $\min_{x\in\mathbb{R}^{n}}\ F_{t}(x;x^{k})$ is not available. Later we will propose more practical stopping criteria that can guarantee (LACC) and (HACC). The connections of our approach to existing work are listed below.

• (a)

  Our (LACC) condition coincides with the inexact stopping criterion proposed in [24, 25, 26] for the proximal gradient method. In these papers, the authors focus on a different optimization problem

  $$\min_{x\in\mathbb{R}^{n}}f_{0}(x):=f_{1}(x)+f_{2}(x),$$

  in which $f_{1}$ is a smooth function, and $f_{2}$ is a proper, convex, and lower semi-continuous function. One typical iteration of their algorithms can be written as

  	$\displaystyle y^{k+1}\approx$	$\displaystyle\ \min_{x\in\mathbb{R}^{n}}f_{0k}(x)=f_{1}(x^{k})+(x-x_{k})^{\top}\nabla f_{1}(x^{k})$
  		$\displaystyle+f_{2}(x)+\frac{1}{2}(x-x^{k})^{\top}H_{k}(x-x^{k}),$		(10a)
  	$\displaystyle x^{k+1}=$	$\displaystyle\ x^{k}+\lambda_{k}(y^{k+1}-x^{k}),$		(10b)

  where $H_{k}\in\mathbb{R}^{n\times n}$ is a positive semi-definite matrix and $\lambda_{k}\in[0,1]$ is a step size. The stopping criterion for inexactly solving (10a) proposed in [24, 25, 26] is

  $$f_{0k}(y^{k+1})-\tilde{f}_{0k}\leq\eta\left(f_{0k}(x^{k})-\tilde{f}_{0k}\right),$$

  (11)

  where $\tilde{f}_{0k}=\min_{x\in\mathbb{R}^{n}}f_{0k}(x)$ and $\eta\in(0,1)$. We note that this is the same as our (LACC). Therefore, our (LACC) is essentially an extension of (11) from the proximal gradient method to the proximal linear method.

• (b)

  To the best of our knowledge, our (HACC) criterion is new and serves as a good alternative to (LACC). From our numerical experiments, we found that (HACC) works comparably with (LACC), and we believe that it can be useful for other applications.

• (c)

  We analyze the overall complexity and the local convergence of our IPL algorithm for solving RPR under the sharpness condition. To the best of our knowledge, this is the first time such results have been obtained under the sharpness condition.

• (d)

  We propose to solve (5) inexactly using FISTA [20, 21, 22], which uses easily verifiable stopping conditions that can guarantee (LACC) and (HACC). Through extensive numerical experiments, we demonstrate that our IPL with the new stopping criteria significantly outperforms existing algorithms for solving RPR.

Organization. The rest of this paper is organized as follows. In Section 2, we propose the main framework of our inexact proximal linear algorithm with two new adaptive stopping criteria for the subproblem. We establish its iteration complexity for obtaining an $\epsilon$-stationary point and its local convergence under the sharpness condition. Connections with some existing methods are also discussed. In Section 3, we discuss how to adapt the FISTA to solve the subproblem inexactly. We also establish the overall complexity of FISTA – the total number of iterations of the FISTA – in order to obtain an $\epsilon$-optimal solution under the sharpness condition. In Section 4, we show the numerical results on both synthetic and real datasets to demonstrate the advantage of the proposed methods over some existing methods. The proofs for all the theorems and lemmas are given in Section 5. Finally, we include some concluding remarks in Section 6.

In this section, we introduce our IPL algorithm for solving the RPR (4) with the inexact stopping criteria (LACC) and (HACC) for the subproblem (5) and analyze its convergence. We will discuss the FISTA for solving (5) that guarantees (LACC) and (HACC) in the next section.

We first follow [15] to introduce some notation. Let

$$S_{t}(y):=\operatorname*{arg\,min}_{x\in\mathbb{R}^{n}}\ F_{t}(x;y),$$

$$\varepsilon(x;y):=F_{t}(x;y)-F_{t}(S_{t}(y);y),$$

where $F_{t}(z;y)$ is defined in (7). We will also use the notation

$$L=\frac{2}{m}\|A\|_{2}^{2}=\frac{2}{m}\left\|\sum_{i=1}^{m}a_{i}a_{i}^{\top}\right\|_{2}.$$

A Meta Algorithm of our IPL is summarized in Algorithm 1. We again emphasize that Algorithm 1 cannot be implemented because $\min_{x\in\mathbb{R}^{n}}\ F_{t}(x;x^{k})$ is not available, and we will discuss practical versions of it in the next section.

In this section, we propose to use the FISTA to inexactly solve (5) with more practical stopping criteria that guarantee (LACC) and (HACC). Therefore, the convergence results (Theorems 1 and 2) in Section 2 still apply here.

For simplicity, we let $t=1/L$ throughout this section and rewrite (5) as follows.

$$\min_{z\in\mathbb{R}^{n}}\ H_{k}(z)=\frac{1}{2t}\|z\|_{2}^{2}+\|B_{k}z-d_{k}\|_{1},$$

(17)

where we denote $z=x-x^{k},$ $B_{k}=\frac{2}{m}\mbox{diag}(Ax^{k})A,$ and $d_{k}=\frac{1}{m}\left(b-(Ax^{k})^{2}\right)$. As a result, (LACC) and (HACC) can be rewritten respectively as

$$H_{k}(z_{k})-\min_{z\in\mathbb{R}^{n}}H_{k}(z)\leq\rho_{l}\left(H_{k}(0)-H_{k}(z_{k})\right),\rho_{l}\geq 0,$$

(18)

and

$$H_{k}(z_{k})-\min_{z\in\mathbb{R}^{n}}H_{k}(z)\leq\frac{\rho_{h}}{2t}\|z_{k}\|_{2}^{2},0\leq\rho_{h}<1/4.$$

(19)

In IPL, we set $x^{k+1}=x^{k}+z_{k}$, where $z_{k}$ satisfies either (18) or (19). The dual problem of (17) is

$$\max_{\lambda\in R^{m},\|\lambda\|_{\infty}\leq 1}D_{k}(\lambda)=-\frac{t}{2}\left\|B_{k}^{\top}\lambda\right\|_{2}^{2}-\lambda^{\top}d_{k}.$$

(20)

From weak duality, we know that

$\displaystyle D_{k}(\lambda)\leq H_{k}(z),\forall z\in\mathbb{R}^{n},\mbox{ and }\|\lambda\|_{\infty}\leq 1.$

(21)

Therefore, $D_{k}(\lambda)$ can serve as a lower bound for $\min_{z}H_{k}(z)$, and we can obtain verifiable stopping criteria that are sufficient conditions for (18) and (19). This leads to our inexact FISTA for solving (17), which is summarized in Algorithm 2. Here we define $z_{k}(\lambda)=-tB_{k}^{\top}\lambda$.

We now discuss the overall complexity of the IPL (Algorithm 1) with the subproblem (5) solved inexactly by FISTA (Algorithm 2). For ease of presentation, we denote this algorithm as IPL+FISTA. We assume that IPL is terminated after $K_{\epsilon}$ iterations, when an $\epsilon$-optimal solution is found, i.e., $\Delta(x^{K_{\epsilon}})\leq\epsilon,\Delta(x^{K_{\epsilon}-1})>\epsilon$. We use $J_{k},k\geq 0$ to denote the number of iterations of FISTA when it is called in the $k$-th iteration (getting $x^{k+1}$ from $x^{k}$) of IPL. The overall complexity of IPL+FISTA for obtaining an $\epsilon$-optimal solution is thus given by $J(\epsilon)=\sum_{k=0}^{K_{\epsilon}-1}J_{k}$, which equals the times that we call (22). Now we are ready to give the overall complexity of IPL+FISTA in Theorem 3.

Theorem 3 shows that under the (LACC-FISTA) stopping criterion, we need $O(1/\epsilon)$ iterations to find an $\epsilon$-optimal solution, and under the (HACC-FISTA) stopping criterion, we have the same rate with regard to a countable positive sequence that decreases to zero. Theorem 3 provides better theoretical rates compared to $O(1/\epsilon^{3})$ in [15]. Moreover, the results in Theorem 3 are about the convergence to $\epsilon$-optimal solution, while the results in [15] are for convergence to $\epsilon$-stationary point. Our results require the sharpness condition, which was not assumed in [15].

For (b) in Theorem 3, we can only find a countable sequence of diminishing $\epsilon_{i}$’s to show the $O(1/\epsilon_{i})$ rate. We cannot show the $O(1/\epsilon)$ rate for any fixed $\epsilon>0$. This is because of the local quadratic convergence under (HACC) shown in (b) of Theorem 3. For instance, if our initial point $x^{0}$ satisfies $\Delta(x^{0})=2\epsilon>\epsilon$, under the (HACC), the quadratic convergence result in Theorem 3 (b) implies that $\Delta(x^{1})\leq C\epsilon^{2}<\epsilon$ for some constant $C>0$ when $\epsilon$ is sufficiently small. Therefore, IPL finds an $\epsilon$-optimal stationary point with only one main iteration. However, Lemma 15 (b) indicates that $J(\epsilon)=J_{0}=O(1/\epsilon^{2})$. Therefore, it is possible that the overall complexity becomes $O(1/\epsilon^{2})$, which is higher than $O(1/\epsilon)$.

In this section, we conduct numerical experiments to compare our IPL method with existing methods for solving the RPR problem (4). Readers can find the code and datasets to replicate the experiments in this section via https://github.com/zhengzhongpku/IPL-code-share. The algorithms that we test include the following ones.

• (i)

  PL: The original proximal linear algorithm proposed by [11] where the subproblem (5) is solved by POGS [27]. POGS terminates when both the primal residual and the dual residual are small enough. In their code, the authors [11] implemented a two-stage trick that uses a relatively larger tolerance in early iterations and a smaller tolerance in later iterations to terminate the POGS. In our comparison, we use all the default parameters set by the authors in their code³³3The code of [11] can be downloaded from https://web.stanford.edu/~jduchi/projects/phase-retrieval-code.tgz.

• (ii)

  Subgradient method. The subgradient method with geometrically decaying step sizes was proposed by [13], and they used this algorithm to solve the RPR (4). One typical iteration of this algorithm is

  $\displaystyle x^{k+1}=x^{k}-\lambda_{0}q^{k}\xi_{k}/\|\xi_{k}\|_{2},k\geq 0,$

  (27)

  in which $\lambda_{0}>0,q\in(0,1)$ are hyper-parameters and $\xi_{k}=\frac{1}{n}\sum_{i=1}^{m}2a_{i}^{\top}x^{k}\mbox{sign}((a_{i}^{\top}x^{k})^{2}-b_{i}).$

• (iii)

  IPL-FISTA-Low, IPL-FISTA-High: our IPL+FISTA algorithm with stopping criteria (LACC-FISTA) and (HACC-FISTA) in Algorithm 2, respectively, and we also used the Armijo backtracking line search discussed in Remark 1.

The initial point for all the tested algorithms is generated by the spectral initialization given in Algorithm 3 in [11]. All the code is run on a server with Intel Xeon E5-2650v4 (2.2GHz). Each task is limited to a single core – no multi-threading is used.