Robust Phase Retrieval by Alternating Minimization

Phase retrieval refers to the recovery of unknown signals $\bm{x}_{\star}\in\mathbb{R}^{d}$ (or $\mathbb{C}^{d}$) from the magnitudes of its linear measurements, which are formulated as

where $\bm{a}_{1},\dots,\bm{a}_{m}\in\mathbb{R}^{d}$ (or $\mathbb{C}^{d}$) and are known measurement vectors. Solving the set of nonlinear equations in (1) arises in numerous applications including X-ray crystallography, diffraction and array imaging, and optics (e.g. [2, 3, 4, 5]). We consider the robust phase retrieval from the noisy amplitude measurements in (1) corrupted with sparse noise, i.e.

where $I_{\mathrm{out}}\subset[m]$ and $I_{\mathrm{in}}=[m]\setminus I_{\mathrm{out}}$ collect the unknown indices of outliers and inliers respectively, and $\{\xi_{i}\}_{i\in I_{\mathrm{out}}}$ is an arbitrary sequence in $\mathbb{R}$. For example, such a scenario arises in phase retrieval imaging applications [6] due to various reasons including detection failures and recording errors.

A suite of methods designed for the plain phase retrieval [7] has been adapted to address the outliers. These methods provide not only empirically successful performances but also theoretical analyses under random measurement models. For instance, anchored regression [8] and PhaseMax [9] formulate phase retrieval given an initial estimate as a linear program. RobustPhaseMax [10] modifies these methods to offer robust estimation by introducing an auxiliary variable to describe the outliers. In another example, Reshaped Wirtinger Flow (RWF) [11] and Amplitude Flow [12] follow a subgradient descent approach for a least squares estimator (LSE). Median-RWF [13] is a variant of these methods tailored to robust phase retrieval. Specifically, Median-RWF uses a truncation type of regularization that identifies and excludes outliers in each iteration by median-based thresholding on the consistency of the current estimate to the measurements. Median-RWF significantly improves the empirical performance of RobustPhaseMax by tolerating a higher fraction of outliers. However, the regularization of Median-RWF involves algorithmic parameters that have been tuned specifically for the Gaussian measurement model. However, it has not been discussed how to generalize the tuning parameters to other measurement models.

A recent work proposed an approach to robust phase retrieval in the classical robust regression framework in statistics [14]. Instead of the least squares, they adopted the least absolute deviation (LAD) [15] to enforce the consistency to the squared amplitude measurements with outliers. The parameter estimation is then cast as a nonconvex optimization problem. They proposed a prox-linear method that updates the estimate iteratively through local linearization of the forward model. This algorithm can be viewed as a variant of the Gauss-Newton method that regularizes the updates with the proximity to the previous iterate. The prox-linear algorithm iteratively refines the estimate through a sequence of quadratic programs for prox-linear problems and provides comparable performance to Median-RWF. Importantly, the Gauss-Newton method does not involve any tuning parameter. However, for large-scale applications such as those in astronomical or medical imaging, further acceleration of this iterative method is desired. They developed the proximal operator graph splitting (POGS) solver for this purpose.

In this paper, we propose an optimization approach to robust phase retrieval that shares strong theoretical guarantees (high tolerance of outlier ratio and no tuning parameters) with the prox-linear algorithm and further improves its computational cost. The objective is achieved by a simple unconstrained Gauss-Newton method for LAD. The resulting optimization is equivalent to an alternating minimization algorithm for LAD, as described in [16], which is solved by a sequence of linear programs. Since this alternating minimization approach is robust in the presence of outliers, we refer to the optimization as Robust-AM. Since this alternating minimization is a robust estimator in the presence of outliers, we refer to the optimization as Robust-AM Our main theoretical result demonstrates that a suitably initialized Robust-AM converges to the ground-truth signal linearly from $m=\mathcal{O}(d)$ random amplitude-only measurements including up to $25\%$ outliers. The desired initialization can be obtained by the existing robust spectral estimators [13, 14]. We verified through comprehensive numerical simulations that Robust-AM empirically outperforms the existing methods for robust phase retrieval. Particularly, it can tolerate a higher fraction of outliers and provide exact recovery with fewer observations. Furthermore, due to its unconstrained optimization formulation with the absolute amplitude measurement model, Robust-AM admits a computationally efficient ADMM algorithm, which runs faster than POGS for the prox-linear method. As shown in Figure 1, ADMM for Robust-AM converges faster than POGS for the prox-linear method. In this experiment, the fraction of outliers is set to $\eta:={|I_{\mathrm{out}}|}/{m}=0.3$, with outlier entries generated following zero and a Cauchy distribution with median $0$ and mean-absolute-deviation $1$. The convergence is measured by the metric $\mathrm{dist}(\mathbf{x},\mathbf{x}_{\star}):=\min_{\alpha\in\{\pm 1\}}\|\mathbf{x}-\alpha\mathbf{x}_{\star}\|_{2}$ for $\mathbf{x},\mathbf{x}_{\star}\in\mathbb{R}^{d}$. Figure 1 shows that the unconstrained Gauss-Newton method, without any explicit control over the proximity to previous iterates, converges to the ground truth signal $\mathbf{x}_{\star}$ without overshooting.

We consider the minimization of the composite function $\ell=h\circ F$ where $h:\mathbb{R}^{m}\to\mathbb{R}$ is a convex function and $F:\mathbb{R}^{d}\to\mathbb{R}^{m}$ is a nonlinear mapping. In the special case when $F$ is differentiable, Burke and Ferris [20] proposed a constrained Gauss-Newton method where the amount of the update is upper-bounded by a threshold. Duchi and Ruan [14] considered a variant where the constraint on the proximity on consecutive iterates is substituted by regularization with an additive penalty. We consider a more challenging case where $F$ is non-differentiable and propose an unconstrained Gauss-Newton method where the variable sequence $(\bm{x}_{k})_{k\in\mathbb{N}\cup\{0\}}$ is iteratively updated by

where $F^{\prime}(\bm{x}_{k})\in\mathbb{R}^{m\times d}$ denotes the Clarke’s generalized Jacobian matrix at $\bm{x}_{k}$ [21]. Due to the local linear approximation of $F$ at $\bm{x}_{k}$ in (3), $\bm{x}_{k+1}$ is obtained as a solution to a convex program. In a special case where $h:\mathbb{R}^{m}\to\mathbb{R}$ and $F:\mathbb{R}^{d}\to\mathbb{R}^{m}$ are respectively given by

and

their composition reduces to

Then the minimization of $\ell$ corresponds to the LAD approach to robust phase retrieval with the absolute amplitude measurement model. Furthermore, given $h$ and $F$ as in (4) and (5), the update rule in (3) is explicitly written as

The resulting algorithm (7), derived from an unconstrained Gauss-Newton method of robust phase retrieval, is equivalent to an alternating minimization approach to the LAD formulation of robust phase retrieval when noisy measurements with a negative sign are discarded. An analogous alternating minimization for least-squares phase retrieval has been studied in the literature [16, 22]. Due to the robustness of LAD, we refer to the iterative algorithm by (7) as a robust alternating minimization (Robust-AM).

Duchi and Ruan [14] considered a similar robust phase retrieval with the squared amplitude measurement model via their regularized Gauss-Newton method.

This section discusses numerical algorithms for Robust-AM. First, we note that the optimization in (7) is equivalent to a linear program

where $\bm{1}_{m}=[1,\dots,1]^{\mathsf{T}}\in\mathbb{R}^{m}$. There exist various computationally efficient numerical methods to solve linear programs. For example, the derandomized algorithm by van den Brand [19] finds an exact solution to a linear program with $d$ variables and $m$ constraints at the cost of $\widetilde{\mathcal{O}}\left((m+d)^{c}\right)$ multiplications where $c\approx 2.38$.

To further accelerate the convergence of Robust-AM, we also adopt iterative numerical algorithms that provide an approximate solution to the inner optimization in (7). In particular, we consider two alternating direction method of multipliers (ADMM) algorithms and a subgradient descent algorithm for the inner optimization. We refer to the Robust-AM with approximate solutions to the inner optimization by these ADMM algorithms as fast Robust-AM since they provide a significantly lower computational cost for the entire convergence of Robust-AM to an $\epsilon$-accurate estimate.

In this section, we present the convergence analysis of the Robust-AM algorithms under the following assumptions. First, we adopt the standard random linear measurements and outliers with arbitrary support and adversarial values [14].

Additionally, to provide the convergence analysis of the fast Robust-AM, we introduce an extra assumption that quantifies the suboptimality of solving (13) by ADMM.

The Robust-AM algorithm updates iterates with an exact solution to (7). Therefore, setting $\epsilon_{\max}$ to $0$ in Theorem IV.1 provides a sufficient condition for the exact recovery of $\bm{x}_{\star}$ by Robust-AM. We compare the specification of Theorem IV.1 to this scenario to the analogous results for competing methods: RobustPhaseMax [10], Median-RWF[13], and prox-linear [14]. Theorem IV.1 as well as the previous results achieve the exact recovery when the number of observations $m$ exceeds a multiple of the signal dimension $d$. Earlier theoretical results on RobustPhaseMax and Median-RWF showed that there exists an unspecified numerical constant so that the algorithms provide the exact recovery if the outlier fraction is below this constant. In contrast, the analyses of the prox-linear [14] and Robust-AM (Theorem IV.1) demonstrate that these methods can tolerate outliers up to $1/4$ of the total observations. These theoretical guarantees consider different degrees of adversary for their outlier models. The performance guarantee of RobustPhaseMax by Hand [10] assumed the highest adversary so that both the support and values of sparse noise are adversarial. The performance guarantees of Median-RWF by Zhang et al. [13] considered the same outlier model as in Assumption IV, but they also introduced additive noise of a bounded norm in addition to sparse noise. Duchi and Ruan [14] used the lowest adversary so that the support of sparse noise is random but the nonzero values of sparse noise can depend on the measurements. Despite providing performance guarantees under the highest adversary, as shown in Section V, RobustPhaseMax showed significantly inferior empirical performance relative to the other methods in terms of the tolerable outlier ratio.

Theorem IV.1 establishes a local linear convergence of the Robust-AM algorithms. As discussed in Section II, Robust-AM has no explicit control over the amount of the update in each iteration unlike the constrained or regularized versions of the Gauss-Newton method [20, 14]. However, despite its simple form, Robust-AM provides the monotone decrease of the estimation error toward zero without any overshooting for robust phase retrieval in the setting of Theorem IV.1. All convergence analyses by Theorem IV.1 and previous work [13, 14] require an initialization within a neighborhood of the ground truth. The size of the basin of convergence was determined with an explicit numerical constant only in [10] and Theorem IV.1. Various initialization methods with theoretical performance guarantees have been developed to obtain the desired initial estimate [13, 14]. The sample complexity for these initialization methods does not exceed those for the subsequent estimators in order.

Next, we discuss the computational costs for the robust estimators. First, RobustPhaseMax is formulated as a linear program and thus it can be exactly solved with $\widetilde{\mathcal{O}}((m+d)^{2.38}\log(1/\epsilon))$ multiplications by derandomized algorithm [19]. Furthermore, as we discussed in Section III-B, there exists an ADMM algorithm for the linear program that costs $\mathcal{O}(m^{3}+(m+d)^{2}\log(1/\epsilon))$ for an $\epsilon$-accurate solution. Due to the term $\log(1/\epsilon)$, if the desired accuracy decreases in proportion to the size of the problem, it is preferable to use ADMM. Otherwise, the derandomized algorithm will be computationally efficient. The other estimators are given as an iterative algorithm with a proven convergence rate. Therefore, we compare their computational costs to obtain an $\epsilon$-accurate solution. Median-RWF is a truncated gradient descent with the per-iteration cost of $\mathcal{O}(md)$. Since the linear convergence of Median-RWF has been established, the total cost is $\mathcal{O}(md\log(1/\epsilon))$. Unlike Median-RWF, the updates in prox-linear and Robust-AM involve a nontrivial inner optimization, respectively cast as a quadratic program and a linear program. One may use an exact solver for these sub-problems. For example, there exists an interior point method for quadratic programs with the cost $\mathcal{O}((m+d)^{4})$ [27]. Since it has been shown that prox-linear converges quadratically, the total cost with this exact inner solver is $\mathcal{O}((m+d)^{4})\log\log(1/\epsilon)$. The inner optimization in Robust-AM can be exactly solved at the cost $\widetilde{\mathcal{O}}((m+d)^{2.38}\log(1/\epsilon))$ by the derandomized algorithm [19]. Due to its linear convergence, the total cost of Robust-AM is $\widetilde{\mathcal{O}}((m+d)^{2.38}\log(1/\epsilon))$. However, as shown in Theorem IV.1, the linear convergence of Robust-AM remains valid when the inner optimization problems are solved only approximately. The fast Robust-AM with the ADMM solver for linear programs has the per-iteration cost of $\mathcal{O}(m^{3}+(m+d)^{2}\log(1/\epsilon_{\max}))$ as shown in Section III. Due to its linear convergence in Theorem IV.1, the total cost to obtain the $\epsilon+\lambda_{\eta}\epsilon_{\max}$ accuracy is $\mathcal{O}(m^{3}+(m+d)^{2}\log(1/\epsilon_{\max})\log(1/\epsilon))$. In contrast, the convergence rate of POGS for the inner optimization in prox-linear has not been established. We summarize the comparison for the computational costs of algorithms in Table I.

Lastly, we elaborate on the dependence of the parameters $\nu_{\eta}$ and $\lambda_{\eta}$ in Theorem IV.1 on the outlier ratio $\eta$. The linear convergence parameter $\nu_{\eta}$ in (15) is explicitly specified as an increasing function of $\eta$ in the proof of Theorem IV.1 and illustrated in Figure 2(a). Therefore, smaller $\eta$ implies faster convergence. The final error bound by (15) with $k$ going to infinity is given as the amplification of the sub-optimality parameter $\epsilon_{\max}$ in the inner optimization by a factor of $\lambda_{\eta}$. First, similar to $\nu_{\eta}$, the parameter is also explicitly given as an increasing function of $\eta$ in the proof (see Figure 2(b)). However, the final estimation can still be sufficiently small, as one can set the accuracy parameter to a sufficiently low value (less than $10^{-10}$) using linear program packages in readily available software such as CPLEX and Gurobi. Hence, the assumption on $\{\epsilon_{i}\}_{i=1}^{k}$ in Theorem IV.1 is easily satisfied.

This section compares the empirical performances of Robust-AM to its theoretical analysis in Theorem IV.1. Robust-AM is also compared against the competing methods for robust phase retrieval, which are RobustPhaseMax, Median-RWF, and the prox-linear. Recall that all these methods require an initial estimate. For this purpose, we adopt the spectral method by Zhang et al. [13].