Smoothed Robust Phase Retrieval

In the robust phase retrieval problem (2), we observe vectors $\{a_{i}\}$ and non-negative scalars $\{b_{i}\}$ and aim to recover the true phases $\{x_{\star},-x_{\star}\}$¹¹1Due to the quadratic relationship in (2), both $x_{\star}$ and $-x_{\star}$ are true signal vectors..

In Assumption 1, $n=n_{1}+n_{2}$ represents the total number of measurements, and $p_{\text{fail}}=\frac{n_{2}}{n_{1}+n_{2}}$ denotes the proportion of corrupted measurements. The sets $\mathbb{I}_{1}=[n_{1}]$ and $\mathbb{I}_{2}=[n]-[n_{1}]$ correspond to the inliers and outliers among the measurements, respectively. In the special case where $\gamma=p_{\text{fail}}=0$, the model reverts to the noiseless phase retrieval scenario as depicted in (1). Our goal is to recover the true signal vectors $x_{\star}$ or $-x_{\star}$ from the observed $\{(a_{i},b_{i})\}_{i=1}^{n}$, without prior knowledge of which measurements are corrupted. This setup ensures the mutual independence of $\{a_{i},i\in[n_{1}]\}$ and $\{a_{i},i\in[n_{1}+n_{2}]-[n_{1}]\}$. We will introduce stronger assumptions on the distribution of $a_{i}$ later. The bounded noise variables $\tau_{i}$s require no additional distributional assumptions, whereas no assumptions are made regarding the distribution of $\{\xi_{i}\}_{i=1}^{n}$ or its independence with respect to $\{a_{i}\}$, allowing for potentially heavy-tailed and adversarially distributed corruption. This assumption is more general then the model assumption $M_{2}$ of [11], which requires independence between $\{\xi_{i},i\in[n_{1}+n_{2}]-[n_{1}]\}$ and $\{a_{i},i\in[n_{1}]\}$. For simplicity, in subsequent discussions, we will refer to a random vector $a$ to denote the generic distribution of the $a_{i}$s.

As shown in [11], when $\gamma=0$, exact recovery is achievable in the optimization problem (3) employing an $\ell_{1}$-loss, indicating that $\{x_{\star},-x_{\star}\}=\operatorname*{arg\,min}_{x\in\mathbb{R}^{p}}F(x)$, with high probability, under certain assumptions. To refine this approach, we employ convolution-type smoothing, focusing on the problem formulated as:

$$\min_{x\in\mathbb{R}^{p}}F_{\delta}(x)=\frac{1}{n}\sum_{i=1}^{n}l_{\delta}\left((a_{i}^{\top}x)^{2}-b_{i}\right),$$

(5)

where

$$l_{\delta}(x)=(l*K_{\delta})(x)=\int_{-\infty}^{+\infty}|y|K_{\delta}(x-y)dy,$$

(6)

$K(x)$ represents a kernel function, $\delta>0$ is the bandwidth, $*$ denotes the convolution operation, and $K_{\delta}(x)=\frac{1}{\delta}K(\frac{x}{\delta})$. In this context, $l_{\delta}(x)$ serves as a smoothed loss function, substituting $l(x)$ for $F(x)$. The adoption of convolution-type smoothing is motivated by gradient approximations: $\nabla F_{\delta}(x)$ closely mirrors the subgradient of $F(x)$ because $\lim_{\delta\rightarrow 0}l^{\prime}_{\delta}(x)=\text{sign}(x)\in\partial|x|$, thereby suggesting that the statistical behavior of $F_{\delta}(x)$ will approximate that of $F(x)$. Moreover, the smoothness of $F_{\delta}(x)$ significantly benefits gradient-based optimization and landscape analysis. Provided some mild conditions on the kernel function $K(x)$, both $l^{\prime}_{\delta}(x)$ and $l^{\prime\prime}_{\delta}(x)$ are well-defined, ensuring that $F_{\delta}(x)$ is smooth and equipped with gradients and Hessians for both optimization and landscape analysis purposes.

In this subsection, we elaborate on the convolution-type smoothed loss function as defined in (6). Additionally, we introduce the notations $\tilde{K}(x)=\int_{-\infty}^{x}K(y)dy$ and $\tilde{K}_{\delta}(x)=\tilde{K}(x/\delta)$ for the integrated kernel function and its scaled version, respectively. To facilitate our analysis, we adhere to the subsequent assumption regarding the kernel function $K(x)$ throughout this paper:

Next, we give some examples of the kernel $K(x)$ that satisfies Assumption 2.

• •

  Gaussian kernel $K(x)=(2\pi)^{-1/2}e^{-x^{2}/2},$ which generates $l_{\delta}(x)=\delta l^{\mathrm{G}}(x/\delta)$ with $l^{\mathrm{G}}(x):=(2/\pi)^{1/2}e^{-x^{2}/2}+x\{1-$ $2\Phi(-x)\}$ and $\Phi(\cdot)$ being the CDF of $N(0,1)$.

• •

  Logistic kernel $K(x)=e^{-x}/\left(1+e^{-x}\right)^{2},$ which generates $l_{\delta}(x)=\delta l^{\mathrm{L}}(x/\delta)$ with $l^{\mathrm{L}}(x)=x+2\log\left(1+e^{-x}\right)$.

• •

  Epanechnikov kernel $K(x)=(3/4)\left(1-x^{2}\right)\mathbf{1}(|x|\leq 1),$ which generates $l_{\delta}(x)=\delta l^{\mathrm{E}}(x/\delta)$ with $l^{\mathrm{E}}(x)=\left(3x^{2}/4-x^{4}/8+3/8\right)\mathbf{1}(|x|\leq 1)+|x|\mathbf{1}(|x|>1)$.

• •

  Triangular kernel $K(x)=(1-|x|)\mathbf{1}(|x\leq 1)$, which generates $l_{\delta}(x)=\delta l^{\mathrm{T}}(x/\delta)$ with $l^{\mathrm{T}}(u)=(x^{2}-|x|^{3}/3+1/3)\mathbf{1}(|x|\leq 1)+|x|\mathbf{1}(|x|>1)$.

• •

  $K(x)=\frac{1}{2(x^{2}+1)^{3/2}},$ which generates the Pseudo-Huber loss $l_{\delta}(x)=\sqrt{x^{2}+\delta^{2}}$.

In practical implementations, $\delta$ ought to be a small constant so that $l_{\delta}(x)$ approximates the $\ell_{1}$-loss, following the ideas introduced in smoothed quantile regression [16]. Following Assumption 2, we can derive several properties of $l_{\delta}(x)$, which are briefly discussed here for space consideration. More details are summarized in Lemma LABEL:lemma_property_huber of Appendix LABEL:appendix:property_hubers in the supplement. Similar to [16], we prove that $l_{\delta}(x)$ is convex and has continuous first-order and second-order derivatives, namely, $l^{\prime}_{\delta}(x)=2\tilde{K}_{\delta}(x)-1$ and $l^{\prime\prime}_{\delta}(x)=2K_{\delta}(x).$ This ensures the existence of the gradient and Hessian for $F_{\delta}(x)$. Moreover, note that $l_{\delta}(x)-l_{\delta}(0)\geq\min\left\{\frac{C_{1}}{\delta}x^{2},C_{2}|x|\right\},\forall x\in\mathbb{R}^{p},$ for some positive constants $C_{1}$ and $C_{2}$. This relation suggests that $l_{\delta}(x)$ exhibits characteristics akin to the $\ell_{2}$-loss for values of $x$ near $0$, and resembles the $\ell_{1}$-loss for values of $x$ distant from $0$. This bound acts as the foundation for the generalized sharpness condition to $F_{\delta}(x)$ in the absence of noise, a topic to be elaborated upon in the subsequent subsection.

In this subsection, we discuss the properties of $F_{\delta}(x)$, which will indicate the similarities between $F(x)$ and $F_{\delta}(x)$ and how minimizing $F_{\delta}(x)$ works for robust phase retrieval. We first introduce two parallel assumptions given in [11]. Assumption 3 claims that the distribution of $a_{i}$ is non-degenerate in all the directions, and Assumption 4 describes the tail-behavior of the distribution of $a_{i}$.

With the same Assumptions, we can prove two properties for $F_{\delta}(x)$ which are comparable to properties for $F(x)$ as follows. The first one will be generalized sharpness.

In the literature, the sharpness of a function $f(x)$ bounds the difference between $f(x)$ and its global minimum by the distance between $x$ and its nearest global minimum point. This concept is prevalent in the optimization problems characterized by equation-solving aspects and nonsmooth objective functions. For more details, please see Section 1 of [29]. The generalized sharpness condition extends the concept of sharpness for $F(x)$, as outlined in Corollary 3.1 of [11], which is the special example of (7) when $\Delta(x)\geq\delta$:

$$F(x)-F(x_{\star})\geq\lambda_{s}\Delta(x),\forall x\in\mathbb{R}^{p}.$$

(8)

[11] and [12] proved that (8) holds with high probability for the noiseless robust phase retrieval. [11], [9], [10], and [38] leveraged (8) to show the convergence rate for their algorithms.

The second one is the weak convexity that relies on the Assumption 4.

This relationship also holds when substituting $F_{\delta}(x)$ with $F(x)$ and replacing the gradient with the subgradient. Weak convexity, a notion extensively recognized in optimization literature, serves as a generalization of convexity. It necessitates that $f(x)+\frac{\rho}{2}\|x\|_{2}^{2}$ be convex for a certain positive constant $\rho$. Various definitions for weak convexity are discussed throughout the literature. The formulation presented in Lemma 2 aligns with that in [10]. Furthermore, Corollary 3.2 in [11] defines weak convexity by bounding the deviation between $F(x)$ and its linear approximation, requiring the existence of a positive constant $\rho>0$ such that:

$$\left|F(x)-\frac{1}{n}\sum_{i=1}^{n}\left|(a_{i}^{\top}y)^{2}-b_{i}+2(a_{i}^{\top}y)a_{i}^{\top}(x-y)\right|\right|\leq\frac{\rho}{2}\|y-x\|_{2}^{2},\quad\forall x,y\in\mathbb{R}^{p}.$$

Next, under the case $\gamma=0$, we discuss how generalized sharpness and weak convexity lead to the robustness of using $F_{\delta}(x)$. As $l_{\delta}(x)$ is an approximation of $l(x)$, we can expect the approximation in function values, which is $\sup_{x\in\mathbb{R}}\left|l_{\delta}(x)-x\right|\leq C_{0}\delta$ for some positive constant $C_{0}$. Together with the sharpness of $F(x)$ given in (8), we can conclude that, for some positive constant $C_{1}$,

$$\Delta(\bar{x})\leq C_{1}\delta,\ \forall\bar{x}\in\operatorname*{arg\,min}_{x\in\mathbb{R}^{p}}F_{\delta}(x).$$

This inequality shows the closeness of minimums for $F(x)$ and $F_{\delta}(x)$. Moreover, our local landscape analysis, which aims at identifying and categorizing stationary points near the true signal vectors, allows for a finer characterization of local regions potentially containing a stationary point. When $p_{\text{fail}}=0$, and under the conditions specified in (7) and (9), it is deduced that:

$$-\lambda_{s}\min\{\Delta(x),\frac{1}{\delta}\Delta^{2}(x)\}\geq F_{\delta}(x_{\star})-F_{\delta}(x)\geq\nabla F_{\delta}(x)^{\top}(x_{\star}-x)-\frac{\rho}{2}\Delta^{2}(x).$$

So, when $\delta$ is sufficiently small, there exists positive constants $C_{2},C_{3}$ such that

$$\nabla F_{\delta}(x)^{\top}(x-x_{\star})\geq C_{2}\min\{\Delta(x),\frac{1}{\delta}\Delta^{2}(x)\},$$

(10)

$$\forall x\in\mathbb{R}^{p}\ s.t.\ \Delta(x)\leq C_{3},\|x-x_{\star}\|_{2}=\Delta(x).$$

Similar conclusions hold for the symmetrical area where $\Delta(z)=\|z+x_{\star}\|_{2}$. So, within $R_{0}:=\{z\in\mathbb{R}^{p}:\Delta(z)\leq C_{3}\}$, we can claim that true signal vectors are stationary points and local minima.

In scenarios involving corruptions with $\gamma=0$, we employ the decomposition

$$F_{\delta}(x)=(1-p_{\text{fail}})F_{\delta,1}(x)+p_{\text{fail}}F_{\delta,2}(x),$$

where $F_{\delta,i}(x)$, for $i=1,2$, represents the sample mean of $l_{\delta}\left((a_{i}^{\top}x)^{2}-b_{i}\right)$ for the noiseless and corrupted measurements, respectively. This allows us to express the gradient as

$$\nabla F_{\delta}(x)^{\top}(x-x_{\star})=(1-p_{\text{fail}})\nabla F_{\delta,1}(x)^{\top}(x-x_{\star})+p_{\text{fail}}\nabla F_{\delta,2}(x)^{\top}(x-x_{\star}).$$

The first term on the right-hand side can be bounded in a manner akin to that in (10). For the second term, we will show in Section 4 that, with high probability, its absolute value can be uniformly bounded by $C_{4}p_{\text{fail}}\|x-x_{\star}\|_{2}$ within $R_{0}$, where $C_{4}$ is a positive constant independent of $\delta$ and $p_{\text{fail}}$. This serves as an indication of the impact corrupted measurements have on the gradient. When $p_{\text{fail}}$ is sufficiently small, these bounds suggest that the vicinity $R_{1}$ around the true signal vectors within $R_{0}$ is defined as

$$R_{1}=\{x\in\mathbb{R}^{p}:\Delta(x)\leq r_{0}\delta p_{\text{fail}}/(1-p_{\text{fail}})\},$$

for a positive constant $r_{0}$, and no stationary points are found in $R_{0}\setminus R_{1}$. Therefore, if an algorithm minimizing $F_{\delta}(x)$ converges to a local minimum $x_{0,\star}$ within $R_{0}$, it must be located in $R_{1}$. This illustrates the robustness of employing $F_{\delta}(x)$, as the algorithm converges to a point near the true signal vectors. In a noiseless context, this point coincides with one of the true signal vectors. In cases with corruption, and when $\delta p_{\text{fail}}/(1-p_{\text{fail}})$ is small, the convergence point is close to either $\{x_{\star},-x_{\star}\}$, serving as an effective warm-start for the IPL algorithm that minimizes $F(x)$ with notable local efficiency, thereby exactly finding one of the true signal vectors.

Addressing robustness alone is not sufficient given the nonconvex nature of $F_{\delta}(x)$. We also need to demonstrate the absence of spurious local minima outside the region $R_{0}$ to ensure the efficacy of algorithms with random initialization. In Sections 3 and 4, we conduct landscape analysis to identify $R_{0}$ and to illustrate the benign landscape characteristics beyond $R_{0}$.

Remark: Inner product quantities akin to $\nabla F_{\delta}(x)^{\top}(x-x_{\star})$, which involve the interaction between the (sub)gradient and the deviation from $x$ to the true vector, are pivotal in the landscape analysis across various nonconvex formulations [10, 24, 32, 8]. Specifically, [24] leveraged such quantities in the global landscape analysis of single-index models with the essential assumption of a strictly monotone link function. Given that phase retrieval, with its squared link function not satisfying this monotonicity requirement, usage of inner product quantities in our paper, [32] and [8] is confined to local vicinities. In these works, nonlocal landscape analyses are conducted through the examination of Hessian matrix-related quantities. A distinctive aspect of our study is the incorporation of corruption into the analysis, setting it apart from existing literature. The comprehensive details of this approach will be elaborated upon in Section 4.

Finally, we provide the following normal assumption that is stronger than both 3 and 4.

Assumption 5 was used in [11, 10, 32, 8] to study initialization and landscape analysis for phase retrieval problems.

In this subsection, we discuss the algorithms for minimizing $F_{\delta}(x)$. We will show that under the noiseless situation, with high probability, gradient descent on $F_{\delta}(x)$ enjoys local linear convergence. Before we start, we first give a Lemma that finds a Lipschitz constant $\mu_{max}/\delta$ for $\nabla F_{\delta}(x)$ with high probability.

Given the above lemma, the local linear convergence for gradient descent can be shown in Theorem 1.

In (12), if we let $k=0$ and choose $\Delta(x^{0})>0$, we have $\lambda_{s}/(8\mu_{max})\leq 1$. Thus, we can define $\kappa:=(8\mu_{max})/\lambda_{s}\geq 1$ as the condition number. It is worth pointing out that $\kappa$ is not related to $\delta$ or $p$ so the convergence rate is statistically stable.

Next, we comment on the technical details we use for reaching this result. Linear convergence usually appears in a strongly convex situation. When investigating the Hessian matrix $\frac{\partial^{2}F_{\delta}(x)}{\partial x\partial x^{\top}}$, $\frac{\mu_{max}}{\delta}$ provides an upper bound for $\left\|\frac{\partial^{2}F_{\delta}(x)}{\partial x\partial x^{\top}}\right\|_{2}$. However, it is hard to find a positive lower bound for the smallest eigenvalue for the Hessian matrix even within the local vicinity of the true signal vectors. So, we turn to the Polyak-Lojasiewicz condition introduced in [21], which proves linear convergence of gradient descent algorithm for minimizing a function $f(x)$ by assuming $\frac{1}{2}\|\nabla f(x)\|_{2}^{2}\geq\mu\left(f(x)-\min_{x}f(x)\right)$ with $\mu>0$. We can prove that when (7), (9) and (11) hold with $\lambda_{s},\rho,\mu_{max}>0$, for any $\delta>0$ that is small enough,

$$\frac{1}{2}\|\nabla F_{\delta}(x)\|_{2}^{2}\geq\frac{\mu_{min}}{\delta}\left(F_{\delta}(x)-F_{\delta}(x_{\star})\right),\forall x\in\{x\in\mathbb{R}^{p}:\Delta(x)\leq C\delta\}$$

in which $\mu_{min}$ and $C$ are positive constants that are not related to the selection of $\delta$. With these conditions, we can get the convergence rate.

Linear convergence guarantees that minimizing the smoothed objective function $F_{\delta}(x)$ is not only sufficient for phase retrieval in the noiseless scenario but also offers enhanced efficiency in comparison to PL [11] and IPL [38]. This improvement stems from the methodologies adopted by [11], which solves the subproblem utilizing the POGS algorithm [28], and IPL, which solves the same subproblem through accelerated proximal gradient descent tailored for the dual problem. Both of these algorithms are characterized by their sub-linear convergence rates. Furthermore, a local linear convergence rate has been established across various phase retrieval algorithms, including Wirtinger flow [4], truncated Wirtinger flow [7], truncated amplitude flow [35], and subgradient methods [9, 10]. Smoothed robust phase retrieval thus joins this list of methodologies. In addition, Theorem 1 and (7) indicate that, for any $\varepsilon>0$, the vanilla gradient descent on the smoothed robust objective function $F_{\delta}(x)$ can find an $x$ such that $\Delta(x)\leq\varepsilon$ with at most $O(\kappa\log\frac{1}{\varepsilon})$ iterations. This is better than $O(\kappa^{2}\log\frac{1}{\varepsilon})$ for the subgradient algorithms [9, 10] using the robust $\ell_{1}-$loss. In our numerical experiments in Section 5, we employ monotone accelerated gradient descent (as outlined in Algorithm 1 of [22]) combined with line search techniques to achieve better overall efficiency. These experiments demonstrate that smoothed robust phase retrieval exhibits superior numerical efficiency on noiseless datasets when compared to both PL and IPL.

In this subsection, we provide an overview of the benign landscape for $F_{\delta}(x)$. Its gradient and Hessian matrix help us find and classify the stationary points, and we will discuss different kinds of benign landscapes in this paper.

1. (a)

   Noiseless population landscape. In this case, $p_{\text{fail}}=\gamma=0$ and we focus on the stationary point for $\mathbb{E}\left[F_{\delta}(x)\right]$. Under some assumptions, we exactly locate all the stationary points and prove that only $\{x_{\star},-x_{\star}\}$ are local minimums.

2. (b)

   Noiseless empirical landscape. In this case, $p_{\text{fail}}=\gamma=0$ and we care about the stationary point of $F_{\delta}(x)$. Under some assumptions, we show that with high probability, the stationary points should be close to those for noiseless population landscape, and only $\{x_{\star},-x_{\star}\}$ are local minimums.

3. (c)

   Empirical landscape with infrequent corruptions only. In this case, $p_{\text{fail}}>0,\gamma=0$ and we study the stationary point of $F_{\delta}(x)$. Under some assumptions, we show that the stationary points are close to those of the noiseless population landscape and all the local minimums are close to $\{x_{\star},-x_{\star}\}$ with high probability. Specifically, all the local minimums lie in the set $\left\{x\in\mathbb{R}^{p}:\Delta(x)\leq r_{0}\delta p_{\text{fail}}/(1-p_{\text{fail}})\right\}$ in which $\Delta(x)=\mbox{dist}(x,\{x_{\star},-x_{\star}\})$ and $r_{0}$ is a positive constant. Thus, if an optimization algorithm finds a local minimum of $F_{\delta}(x)$, its gap to $\{x_{\star},-x_{\star}\}$ is small when both $\delta$ and $p_{\text{fail}}$ are small.

4. (d)

   Empirical landscape with infrequent corruptions and bounded noises. When $p_{\text{fail}}>0,\gamma>0$, similar results can be obtained when $\delta,\gamma,p_{\text{fail}}$ are small enough and $\delta\geq\Omega(\sqrt{\gamma})$, with high probability, all the local minimums lie in the set $\{x\in\mathbb{R}^{p}:\Delta(x)\leq\tilde{r}_{0}\max\{\sqrt{\delta\gamma},\delta p_{\text{fail}}/(1-p_{\text{fail}})\}\}$ where $\tilde{r}_{0}$ is a positive constant.

Now, we describe the procedure of smoothed robust phase retrieval. We first apply some gradient-based algorithm to minimize $F_{\delta}(x)$ and find a local minimum $x_{0,\star}$. Due to the benign landscapes described above, for the optimization, we expect that both random initialization and modified spectral initialization work well such that $x_{0,\star}$ is close to $x_{\star}$ or $-x_{\star}$. Second, for the noiseless situation, we can claim that $x_{0,\star}$ is already a true signal vector. For the situation with noise, benign landscape guarantees that $x_{0,\star}$ is close to one of the true signal vectors and is treated as a warm-start for IPL that minimizes $F(x)$ for exact recovery.

Denote $L(x;a)=l_{\delta}((a^{\top}x)^{2}-(a^{\top}x_{\star})^{2}),x,a\in\mathbb{R}^{p}.$ The gradient vector is $\frac{\partial L(x;a)}{\partial x}=\left\{2l^{\prime}_{\delta}\left((a^{\top}x)^{2}-(a^{\top}x_{\star})^{2}\right)aa^{\top}\right\}x$ and the Hessian matrix $\frac{\partial^{2}L(x;a)}{\partial x\partial x^{\top}}=2l^{\prime}_{\delta}((a^{\top}x)^{2}-(a^{\top}x_{\star})^{2})aa^{\top}+4(a^{\top}x)^{2}l^{\prime\prime}_{\delta}((a^{\top}x)^{2}-(a^{\top}x_{\star})^{2})aa^{\top}$. The empirical and population versions of objective functions can be written as $F_{\delta}(x)=\mathbb{E}_{n}\left[L(x;a)\right]$ and $\mathbb{E}\left[F_{\delta}(x)\right]=\mathbb{E}\left[L(x;a)\right].$ They have continuous gradient vectors and Hessian matrices as $\nabla F_{\delta}(x)=\mathbb{E}_{n}\left[\nabla L(x;a)\right]$, $\nabla\mathbb{E}\left[F_{\delta}(x)\right]=\mathbb{E}\left[\nabla L(x;a)\right]$, $\frac{\partial^{2}}{\partial x\partial x^{\top}}F_{\delta}(x)=\mathbb{E}_{n}\left[\frac{\partial^{2}}{\partial x\partial x^{\top}}L(x;a)\right]$, and $\frac{\partial^{2}}{\partial x\partial x^{\top}}\mathbb{E}\left[F_{\delta}(x)\right]=\mathbb{E}\left[\frac{\partial^{2}}{\partial x\partial x^{\top}}L(x;a)\right].$ We will use the Hessian matrix to classify the stationary points and find the benign noiseless population landscapes for $\mathbb{E}\left[F_{\delta}(x)\right]$ and empirical landscapes for $F_{\delta}(x)$.