On Fienup Methods for Regularized Phase Retrieval

TL;DR

This paper provides new theoretical insights into Fienup methods for phase retrieval, showing their convergence properties and proposing an accelerated version that performs well with sparse signals.

Contribution

It offers a novel analysis of Fienup methods as regularized nonconvex optimization and introduces an accelerated algorithm for sparse phase retrieval.

Findings

Fienup methods can be viewed as alternating minimization on a regularized nonconvex problem.

Under semi-algebraic priors, the algorithm converges smoothly to a critical point.

The proposed accelerated method with an $ ext{l}_1$ prior outperforms existing approaches in sparse phase retrieval.

Abstract

Alternating minimization, or Fienup methods, have a long history in phase retrieval. We provide new insights related to the empirical and theoretical analysis of these algorithms when used with Fourier measurements and combined with convex priors. In particular, we show that Fienup methods can be viewed as performing alternating minimization on a regularized nonconvex least-squares problem with respect to amplitude measurements. We then prove that under mild additional structural assumptions on the prior (semi-algebraicity), the sequence of signal estimates has a smooth convergent behaviour towards a critical point of the nonconvex regularized least-squares objective. Finally, we propose an extension to Fienup techniques, based on a projected gradient descent interpretation and acceleration using inertial terms. We demonstrate experimentally that this modification combined with an $\ell_1$ prior constitutes a competitive approach for sparse phase retrieval.

Equations

$$\mathbf{c}=|\mathcal{F}({\mathbf{x}}_0)|+\mathbf{ϵ},$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^n,\mathbf{z}\in Z_{\mathbf{c}}}{\min }\frac{1}{2}\Vert \mathbf{x}-\mathbf{z}{\Vert }_2^2+g(\mathbf{x}),$$

$$\underset{\mathbf{x}\in S}{\min }f(\mathbf{x})$$

$$\underset{\begin{array}{c}x\to \bar{x}\\ x\ne \bar{x}\\ x\in S\end{array}}{\mathrm{liminf}}$$

$${\mathrm{prox}}_g(\mathbf{x})\equiv \underset{\mathbf{y}\in {\mathbb{R}}^n}{\mathrm{argmin}}\left\{\frac{1}{2}\Vert \mathbf{x}-\mathbf{y}{\Vert }_2^2+g(\mathbf{y})\right\}.$$

$${\mathbf{z}}^{k+1}$$

$${\mathbf{x}}^{k+1}$$

$$\hat{\mathbf{z}}[j]=\{\begin{array}{ll}\mathbf{c}[j]\frac{\mathcal{F}(\mathbf{x})[j]}{|\mathcal{F}(\mathbf{x})[j]|},& \text{if }|\mathcal{F}(\mathbf{x})[j]|\ne 0,\\ \mathbf{c}[j]{\mathrm{e}}^{i{\theta }_j},& \text{for an arbitrary }{\theta }_j\text{ otherwise}.\end{array}$$

$${\text{argmin}}_{\mathbf{x}\in {\mathbb{R}}^n}\left\{\frac{1}{2}\Vert \mathbf{x}-\mathbf{z}{\Vert }_2^2+g(\mathbf{x})\right\}={\mathrm{prox}}_g(\mathrm{Re}(\mathbf{z})),$$

$${\mathbf{z}}^{k+1}$$

$${\mathbf{x}}^{k+1}$$

$${\mathbf{x}}^{k+1}=\mathrm{Re}(P_{Z_{\mathbf{c}}}({\mathbf{x}}^k)).$$

$${\mathbf{x}}^{k+1}={\mathcal{T}}_{\lambda }(\mathrm{Re}(P_{Z_{\mathbf{c}}}({\mathbf{x}}^k))).$$

$${\mathbf{x}}^{k+1}\in P_{C_K}(\mathrm{Re}(P_{Z_{\mathbf{c}}}({\mathbf{x}}^k))).$$

$$\underset{\mathbf{z}\in Z_{\mathbf{c}}}{\min }\frac{1}{2}\Vert \mathbf{x}-\mathbf{z}{\Vert }_2^2=\frac{1}{2n}\Vert |\mathcal{F}(\mathbf{x})|-\mathbf{c}{\Vert }_2^2.$$

$$\underset{\mathbf{z}\in Z_{\mathbf{c}}}{\min }\frac{1}{2}\Vert \mathbf{x}-\mathbf{z}{\Vert }_2^2$$

$$=\frac{1}{2}\Vert {\mathcal{F}}^{-1}(\mathcal{F}(\mathbf{x})-\hat{\mathbf{z}})){\Vert }_2^2$$

$$=\frac{1}{2n}\Vert \mathcal{F}(\mathbf{x})-\hat{\mathbf{z}}{\Vert }_2^2.$$

$$|\mathcal{F}(\mathbf{x})[j]-\hat{\mathbf{z}}[j]|=\{\begin{array}{ll}||\mathcal{F}(\mathbf{x})[j]|-\mathbf{c}[j]|,& \text{if }|\mathcal{F}(\mathbf{x})[j]|\ne 0,\\ \mathbf{c}[j],& \mathrm{otherwise}.\end{array}$$

$$\underset{\mathbf{z}\in Z_{\mathbf{c}}}{\min }\frac{1}{2}\Vert \mathbf{x}-\mathbf{z}{\Vert }_2^2=\frac{1}{2n}\Vert \mathcal{F}(\mathbf{x})-\hat{\mathbf{z}}{\Vert }_2^2=\frac{1}{2n}\Vert |\mathcal{F}(\mathbf{x})|-\mathbf{c}{\Vert }_2^2,$$

$$\min \left\{F(\mathbf{x})\equiv \frac{1}{2n}\Vert |\mathcal{F}(\mathbf{x})|-\mathbf{c}{\Vert }_2^2+g(\mathbf{x})\right\},$$

$$h(\mathbf{x},\mathbf{y})\equiv \frac{1}{2}\Vert \mathbf{y}-\mathbf{z}(\mathbf{x}){\Vert }_2^2+g(\mathbf{y}).$$

$$h(\mathbf{x},\mathbf{y})$$

$$\ge \frac{1}{2}\Vert \mathbf{y}-\mathbf{z}(\mathbf{y}){\Vert }_2^2+g(\mathbf{y})=F(\mathbf{y}),\forall \mathbf{y}\in {\mathbb{R}}^n,$$

$$h(\mathbf{x},\mathbf{x})$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }\frac{1}{2n}\Vert |\mathcal{F}(\mathbf{x})|-\mathbf{c}{\Vert }_2^2+g(\mathbf{x}).$$

$${\mathbf{x}}^{k+1}$$

$$={\mathrm{prox}}_g(\mathrm{Re}(\mathbf{z}({\mathbf{x}}^k)))$$

$$={\mathrm{prox}}_g(\mathrm{Re}(P_{Z_{\mathbf{c}}}({\mathbf{x}}^k))),$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^n,({\mathbf{w}}_1,{\mathbf{w}}_2)\in {\tilde{Z}}_{\mathbf{c}}}{\min }\left\{\frac{1}{2}\Vert \mathbf{x}-{\mathbf{w}}_1{\Vert }_2^2+\frac{1}{2}\Vert {\mathbf{w}}_2{\Vert }_2^2+g(\mathbf{x})\right\}.$$

$$\underset{({\mathbf{w}}_1,{\mathbf{w}}_2)\in {\tilde{Z}}_{\mathbf{c}}}{\min }\left\{H({\mathbf{w}}_1,{\mathbf{w}}_2)\equiv G({\mathbf{w}}_1)+\frac{1}{2}\Vert {\mathbf{w}}_2{\Vert }^2\right\},$$

$$G({\mathbf{w}}_1)\equiv \underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }\left\{\frac{1}{2}\Vert {\mathbf{w}}_1-\mathbf{x}{\Vert }_2^2+g(\mathbf{x})\right\}.$$

$$\nabla H({\mathbf{w}}_1,{\mathbf{w}}_2)=\left({\mathbf{w}}_1-{\mathrm{prox}}_g({\mathbf{w}}_1),{\mathbf{w}}_2\right).$$

$$\Vert \nabla H({\mathbf{w}}_1,{\mathbf{w}}_2)-\nabla H({\tilde{\mathbf{w}}}_1,{\tilde{\mathbf{w}}}_2){\Vert }_2^2$$

Full text

I Introduction

Phase retrieval is an old and fundamental problem in a variety of areas within engineering and physics [1, 2]. Many applications of the phase retrieval problem involve estimation of a signal from the modulus of its Fourier measurements. This problem is ill posed in general, so that uniqueness and recovery typically require prior knowledge on the input, particularly in one-dimensional problems. Here we focus on the estimation of real sparse signals from their Fourier magnitude, a problem which has been treated in several recent works [3, 4, 5, 6].

A longstanding line of algorithms to tackle the phase retrieval problem involve application of the alternating minimization method which alternate between the constraints in time and the Fourier magnitude constraints [7, 8, 9]. These methods were pioneered by the work of Gerchberg and Saxton and later extended by Fienup; see [10] for an optimization point of view on these techniques and a rich historical perspective. Alternating minimization approaches have also been recently applied to phase retrieval from random measurements [11]. The main advantage of this class of algorithms is their simplicity and scalability.

A more recent approach to phase retrieval is to formulate the recovery as a smooth nonconvex least-squares estimation problem and use dedicated techniques to estimate the signal using continuous optimization algorithms that guarantee convergence to stationary points. The GESPAR algorithm [4] is an example of this approach which is based on the Gauss-Newton method coupled with sparsity priors. For phase retrieval with random measurements, gradient descent methods have been proposed and analyzed such as Wirtinger flow [12] and truncated amplitude flow [13]. Both treat least-squares objectives where Wirtinger flow measures the loss with respect to the squared-magnitude of the measurements while the amplitude flow approach performs a truncated gradient descent on an amplitude objective. Another line of work suggests the use of matrix lifting and semidefinite programming based relaxations [14, 15, 16, 17, 6]. These techniques are limited by the size of problems that can be tackled using available numerical solvers.

Our main contribution is to propose a new look at alternating minimization algorithms for phase retrieval in the context of Fourier measurements and convex priors. We refer collectively to these techniques as Fienup methods. The use of Fourier measurements is less flexible than general measurements and is less suited for statistical analysis. On the other hand, the Fourier transform has very strong structure which allows for richer algorithmic constructions and analysis.

As a first step we provide two new interpretations of Fienup algorithms. First we show that these techniques are naturally linked to a nonsmooth nonconvex least-squares problem with respect to an amplitude objective. Fienup approaches can then be understood as majorization-minimization methods for solving this problem. Second, we demonstrate that Fienup algorithms can be viewed as a projected gradient descent scheme to minimize a smooth convex objective function over a nonconvex constraint set. This observation allows to characterize the behaviour of the algorithm and develop extensions based on known ideas for accelerating gradient methods using inertial terms [18]. We then specialize these results to the case of $\ell_{0}$ and $\ell_{1}$ priors, leading to a new inertial gradient scheme, which we refer to as FISTAPH: FISTA for PHase retrieval.

On the theoretical side, we show that if the convex prior is well structured (semi-algebraic or more generally representable), then the sequence of signal estimates produced by Fienup has a smooth convergence behaviour. Recall that, broadly speaking, an object is said to be semi-algebraic if it can be represented by systems of polynomial inequalities. The notion of smooth convergence is a very desirable property, even more in nonconvex settings where it is usually not possible to obtain global convergence estimates. The convergence analysis follows well established techniques from tame optimization [19, 20]. These techniques build upon the Kurdyka-Łojasiewicz (KL) property which holds for many classes of functions [21, 22, 23, 24]. We then provide numerical experiments based on synthetic problems to compare Fienup with $\ell_{0}$ and $\ell_{1}$ priors, GESPAR [4], Wirtinger flow (or gradient) methods [12] with $\ell_{0}$ and $\ell_{1}$ priors and FISTAPH. Numerical results suggest that the latter combined with an $\ell_{1}$ prior constitutes a very competitive alternative for sparse phase retrieval.

The rest of the paper is organized as follows. Section II introduces our notation and states the problem of interest more formally. We also introduce several mathematical definitions that are required for the rest of the paper and review the numerical algorithms that are used in subsequent sections. Section III describes our characterization of Fienup methods in the context of phase retrieval from Fourier measurements with convex priors. We detail the relation of Fienup with a nonsmooth nonconvex least-squares problem as well as its interpretation as projected gradient descent. Our main convergence result and our new FISTAPH algorithm are presented in Section IV. Simulation results are provided in Section V.

II Problem Formulation and Mathematical Background

II-A Notation

Throughout the paper vectors are denoted by boldface letters. For a vector ${\bf x}\in\mathbb{C}^{n}$, ${\bf x}[i]$ is the $i$-th entry of ${\bf x}$, $i=1,2,\ldots,n$ and ${\rm supp}({\bf x})$ is the support of ${\bf x}$, namely, the set $\left\{i=1,2,\ldots n;\;{\bf x}[i]\neq 0\right\}$. Furthermore, $\|{\bf x}\|_{0}$ counts the number of nonzero entries of the vector ${\bf x}$: $\|{\bf x}\|_{0}=|{\rm supp}({\bf x})|$ and $\|{\bf x}\|_{p}$ denotes the $\ell_{p}$ norm of ${\bf x}$ for $p\in\mathbb{R}_{+}$. The notations $|\cdot|$, ${\rm Re}(\cdot)$, ${\rm Im}(\cdot)$ and $\bar{\cdot}$ describe the modulus, real part, imaginary part and complex conjugate, respectively, defined over the field of complex numbers. If their argument is a vector, then they should be understood component-wise. Similarly, basic operations, e.g. powers, are taken component-wise when applied to vectors. For ${\bf x}\in\mathbb{C}^{n}$ and $N\in\mathbb{N}$, $\mathcal{F}({\bf x},N)\in\mathbb{C}^{N}$ is the vector composed of the $N$ first coefficients of the discrete Fourier transform of ${\bf x}$ (obtained by zero padding if $n<N$). For simplicity, we use the shorthand notation $\mathcal{F}({\bf x})=\mathcal{F}({\bf x},n)$ to denote the standard discrete Fourier transform of ${\bf x}\in\mathbb{C}^{n}$ and $\mathcal{F}^{-1}$ to denote its inverse. For a set $S$, $\delta_{S}\colon S\to\mathbb{R}\cup\{+\infty\}$ is the indicator function of $S$ ([math] if its argument is in $S$, $+\infty$ otherwise) and $P_{S}$ denotes the Euclidean orthogonal projection onto the set $S$.

II-B Phase Retrieval

Given ${\bf x}_{0}\in\mathbb{R}^{n}$, we consider the data acquisition process

[TABLE]

where $\boldsymbol{\epsilon}\in\mathbb{R}^{n}$ is an unknown vector of errors. In the rest of the paper, we actually assume that ${\bf c}$ has positive entries (it is always possible to set the potential negative entries of ${\bf c}$ to zero). The phase retrieval problem consists of producing an estimate $\hat{{\bf x}}\in\mathbb{R}^{n}$ of ${\bf x}_{0}$ based solely on the knowledge of ${\bf c}$ given by (1).

As mentioned in the introduction, phase retrieval of one-dimensional vectors from Fourier measurements requires the use of prior knowledge. We focus on support and sparsity inducing priors. For $J\subseteq\left\{1,2,\ldots,n\right\}$, we define the set $X_{J}=\left\{{\bf x}\in\mathbb{R}^{n};\;{\rm supp}({\bf x})\subseteq J\right\}$. The prior function that we use will be denoted by $g\colon\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$. We focus on the following priors (for a given $J$):

• •

$g\colon{\bf x}\mapsto\|{\bf x}\|_{0}+\delta_{X_{J}}({\bf x})$, or $\ell_{0}$-based nonconvex prior.

• •

$g\colon{\bf x}\mapsto\|{\bf x}\|_{1}+\delta_{X_{J}}({\bf x})$, or $\ell_{1}$-based convex prior.

In the experimental section, we compare between these two classes of priors. The algorithmic derivations in this paper will be made under the assumption that $g$ is proper and lower semicontinuous, and the main convergence result (c.f. Theorem IV.1) will require in addition convexity of $g$. In order to efficiently implement the proposed algorithm, we need to focus on priors for which the proximity operator [25] is easy to compute. We provide several examples of such priors in Section II-D.

In the rest of the paper, ${\bf c}\in\mathbb{R}^{n}_{+}$ denotes modulus measurements which are assumed to be given, fixed and obtained through (1). Given ${\bf c}\in\mathbb{R}^{n}_{+}$, we define $Z_{{\bf c}}=\left\{{\bf z}\in\mathbb{C}^{n};\;|\mathcal{F}({\bf z})|={\bf c}\right\}$ as the set of values ${\bf z}$ that could have produced ${\bf c}$ (ignoring the noise). To estimate ${\bf x}_{0}$, we consider the regularized least-squares problem

[TABLE]

where $g$ encodes our prior knowledge. Our algorithmic approach consists of employing an alternating minimization method, or one of its variants, to solve the above formulation.

II-C Prior Algorithms for Phase Retrieval

We briefly review several existing algorithms for phase retrieval that will be used in our experiments in Section V.

One approach to sparse phase retrieval is the GESPAR algorithm which is based on the damped Gauss-Newton method in conjunction with an $\ell_{0}$ prior [4]. Damped Gauss-Newton allows to solve smooth, nonlinear least-squares problems. The work of [12] is based on the notion of Wirtinger derivatives to treat the same smooth least-squares problem as GESPAR. The notion of Wirtinger derivative is needed since the objective is not differentiable (holomorphic) as a function of complex variables (see [12] for details). In the case of real valued functions of real variables, the Wirtinger derivative reduces to a standard gradient (up to a constant multiplicative factor). An obvious extension of these types of methods is the use of proximal decomposition, or forward-backward methods which consist in alternating a gradient step on the smooth part of the objective with a proximal step on the nonsmooth part [18, 26]. This is the approach that we use in the numerical experiments to treat phase retrieval with priors.

Finally, we consider alternating minimization methods that are the main focus of this work. This approach consists of solving (2) by applying the alternating minimization algorithm. The special structure of the problem allows to perform each partial minimization efficiently. In particular, the projection onto $Z_{{\bf c}}$ is easy, as described below in (5). These types of methods are also referred to as Fienup algorithms. A deeper interpretation of this approach is given in Section III.

II-D Tools from Convex and Nonsmooth Analysis

Throughout the paper, our results will be based on tools from convex and nonsmooth analysis which we review here.

The gradient of a differentiable function $f$ is denoted by $\nabla f$. This concept admits extensions to nonsmooth analysis; the subgradient of a nonsmooth function $g$ is denoted by $\partial g$. For convex functions, subgradients correspond to tangent affine lower bounds. This definition no longer holds for nonconvex functions. In this case, the proper understanding of subgradients involves much more machinery which will not be discussed here. We only consider the notion of a Fréchet critical point which generalizes classical first order criticality for differentiable functions (see [27]).

Definition 1** (Fréchet critical point)**

Let $S\subseteq\mathbb{R}^{n}$ be a closed set and $f\colon\mathbb{R}^{n}\to\mathbb{R}$ be a lower semicontinuous function. We say that $\bar{{\bf x}}\in S$ is a Fréchet critical point of the problem

[TABLE]

if

[TABLE]

In other words, the negative variations of $f$ in $S$ around $\bar{x}$ are negligible at the first order.

We will also heavily use the notion of the proximity operator of a function.

Definition 2** (Proximity operator)**

For a nonsmooth function $g:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$, the (potentially multivalued) proximity operator is denoted by $\operatorname{prox}_{g}$ and defined by

[TABLE]

Note that when $g$ is proper lower semicontinuous and convex, this operator is single valued.

We next provide a few examples of such functions with their proximity operators; many more can be found, for example, in [28].

Example 1** (Proximity operators)**

• •

Support prior:* If $C\subseteq\mathbb{R}^{n}$ is a closed convex set, then $\operatorname{prox}_{\delta_{C}}$ is the Euclidean projection onto $C$. This can be used for example to encode knowledge about the support of the signal ${\bf x}_{0}$ by choosing $C=X_{J}$ for some $J\subseteq\left\{1,2,\ldots,n\right\}$. In this case, the projection simply consists in setting the coefficients ${\bf x}[i]$ to [math] for $i\not\in J$.*

• •

Sparsity prior:* If $g$ is the $\ell_{1}$ norm, then the proximal operator is the soft thresholding operator.111 The soft thresholding operator is given by ${\mathcal{T}}_{\alpha}({\bf x})_{i}={\rm sgn}(x_{i})\max\{|x_{i}|-\alpha,0\}$. If $g({\bf x})=\lambda\|{\bf x}\|_{1}$ for some $\lambda>0$, then $\operatorname{prox}_{g}({\bf x})={\mathcal{T}}_{\lambda}({\bf x})$. This can be combined with support information prior by first setting the coefficients outside of the support to [math] and then applying the soft thresholding operator.*

• •

Change of basis:* Suppose that $D$ is an $n\times n^{\prime}$ real matrix such that its columns form an orthonormal family, that is $D^{T}D$ is the identity in $\mathbb{R}^{n^{\prime}}$. Suppose that $\tilde{g}\colon\mathbb{R}^{n^{\prime}}\to\mathbb{R}$ is a lower semicontinuous convex function and let $g({\bf x})=\tilde{g}(D^{T}{\bf x})$. In this case, we have $\operatorname{prox}_{g}({\bf x})={\bf x}+D\left(\operatorname{prox}_{\tilde{g}}(D^{T}{\bf x})-D^{T}{\bf x}\right)$ (see [28, Table 1]). This allows to express priors in different orthonormal bases, such as wavelets for example.*

It is also worth mentioning that the proximity operator is efficiently computable for some nonconvex priors. For example, if $g=\delta_{C}$ where $C=\left\{{\bf x}\in\mathbb{R}^{n};\;\|{\bf x}\|_{0}\leq k\right\}$, then the proximity operator is obtained by setting the $n-k$ lowest coefficients (in absolute value) to [math]. This can also be combined with support information.

III Fienup, Majorization-Minimization and Projected Gradient

In this section we expand on the alternating minimization approach to (2) leading to the Fienup family of algorithms. For this section, the prior term $g$ in (2) is taken to be a general proper lower semicontinuous function. We begin by describing the algorithm and then provide two interpretations of it.

III-A Alternating Minimization Algorithm

The alternating minimization algorithm applied to problem (2) is explicitly written below.

Alternating Minimization (Fienup)

Initialization. ${\bf x}^{0}\in\mathbb{R}^{n}$

General Step. For $k\in\mathbb{N}$,

$\displaystyle{\bf z}^{k+1}$ $\displaystyle\in\operatorname*{argmin}_{{\bf z}\in Z_{\bf c}}\frac{1}{2}\|{\bf x}^{k}-{\bf z}\|_{2}^{2},$

$\displaystyle{\bf x}^{k+1}$ $\displaystyle\in\operatorname*{argmin}_{{\bf x}\in\mathbb{R}^{n}}\frac{1}{2}\|{\bf x}-{\bf z}^{k+1}\|_{2}^{2}+g({\bf x}).$

(4)

The main interest in this scheme is that both partial minimization steps in (III-A) can be carried out efficiently whenever $g$ is “proximable”, meaning that its prox (or a member in its prox) is easily computed. First consider, in (III-A), the partial minimization in ${\bf z}$ with ${\bf x}\in\mathbb{R}^{n}$ being arbitrary but fixed. This minimization amounts to computing $P_{Z_{{\bf c}}}({\bf x}^{k})$, the orthogonal projection of ${\bf x}^{k}$ onto $Z_{\bf c}$. For a given ${\bf x}\in\mathbb{R}^{n}$, all the members in $P_{Z_{{\bf c}}}({\bf x})$ are of the form ${\bf z}=\mathcal{F}^{-1}(\hat{{\bf z}})$, where for $j=1,2,\ldots,n$, we have ($i=\sqrt{-1}$ in the equation below)

[TABLE]

Next, we treat the subproblem in (III-A) of minimizing with respect to ${\bf x}$ where ${\bf z}\in\mathbb{C}^{n}$ is arbitrary but fixed. The partial minimization in ${\bf x}$ is given by the expression

[TABLE]

where ${\rm Re}$ is the real part taken component-wise. We have used the definition of the proximity operator of $g$ given in (3). When this operator is easy to compute, each step of the algorithm can be carried out efficiently.

The iterations of the alternating minimization method are summarized as follows:

[TABLE]

We now consider several special cases of (III-A):

• •

If $g=0$, then $\operatorname{prox}_{g}$ is the identity and we recover the original algorithm from Fienup [9], or alternating projection [10].

Fienup

Initialization. ${\bf x}^{0}\in\mathbb{R}^{n}$.

General Step. For $k\in\mathbb{N}$,

${\bf x}^{k+1}={\rm Re}(P_{Z_{\bf c}}({\bf x}^{k})).$

The convergence result given in Theorem IV.1 also holds in this case since constant functions are convex and continuous.

• •

If $g({\bf x})=\lambda\|{\bf x}\|_{1}$ for some $\lambda>0$, then $\operatorname{prox}_{g}={\mathcal{T}}_{\lambda}$, where ${\mathcal{T}}_{\lambda}$ is the soft thresholding operator (see footnote on page 1). We refer to the resulting algorithm as “AM L1”.

AM L1

Initialization. ${\bf x}^{0}\in\mathbb{R}^{n}$, $\lambda>0$

General Step. For $k\in\mathbb{N}$,

${\bf x}^{k+1}={\mathcal{T}}_{\lambda}({\rm Re}(P_{Z_{\bf c}}({\bf x}^{k}))).$

• •

If $g=\delta_{C_{K}}$, where $C_{K}$ is the set of all $K$-sparse vectors, $C_{K}=\{{\bf x}\in\mathbb{R}^{n}:\|{\bf x}\|_{0}\leq K\}$, then $\operatorname{prox}_{g}=P_{C_{K}}$ is the so-called hard thresholding operator. This operator outputs a vector which is all zeros except for the largest $K$ components (in absolute values) of its input vector which are kept the same. The hard thresholding operator is multivalued and the resulting algorithm, which we term “AM L0” picks an arbitrary point in its range.

AM L0

Initialization. ${\bf x}^{0}\in\mathbb{R}^{n}$, $K\in\mathbb{N}$

General Step. For $k\in\mathbb{N}$,

${\bf x}^{k+1}\in P_{C_{K}}({\rm Re}(P_{Z_{\bf c}}({\bf x}^{k}))).$

III-B Majorization-Minimization Interpretation

In this section, we focus on partial minimization in ${\bf z}$. We show that the value of this partial minimization leads to a least-squares objective. This allows us to interpret the Fienup algorithm as a majorization-minimization process on this least-squares function. For the rest of this section, for any ${\bf x}\in\mathbb{R}^{n}$, we denote by ${\bf z}({\bf x})$ an arbitrary but fixed member of $P_{Z_{{\bf c}}}({\bf x})$.

III-B1 Partial Minimization in ${\bf z}$

The following lemma provides a connection between partial minimization in ${\bf z}$ and the evaluation of a nonsmooth least-squares objective.

Lemma III.1

For any ${\bf x}\in\mathbb{R}^{n}$, we have

[TABLE]

Proof:

An optimal solution of the minimization problem is given by ${\bf z}=\mathcal{F}^{-1}(\hat{{\bf z}})$ where $\hat{{\bf z}}$ has the form (5). Now,

[TABLE]

Using the expression of $\hat{{\bf z}}$ in (5), we have for all $j=1,2,\ldots,n$,

[TABLE]

Putting everything together,

[TABLE]

which completes the proof. ∎

A direct consequence of Lemma III.1 is the following corollary that connects between problem (2) and a regularized nonlinear least-squares problem.

Corollary III.2

The pair $({\bf x},{\bf z})$ is an optimal solution of problem (2) if and only if ${\bf x}$ is an optimal solution of

[TABLE]

and ${\bf z}=\mathcal{F}^{-1}(\hat{{\bf z}})$, where $\hat{{\bf z}}$ is of the form given in (5).

Note that in (10), the least-squares objective is defined with respect to the amplitude $|\mathcal{F}({\bf x})|$ and not the magnitude-squared $|\mathcal{F}({\bf x})|^{2}$. For random measurements, it has been shown in [13] that the amplitude objective leads to superior performance over the standard magnitude-squared approach.

III-B2 Fienup as Majorization-Minimization

In order to understand further the connection with the Fienup algorithm, we define the following auxiliary function:

[TABLE]

Now, for any ${\bf x}\in\mathbb{R}^{n}$, using Lemma III.1, we have the following properties (recalling the definition of $F$ in (10)):

[TABLE]

In other words, using the convexity of $g$, we have that $h({\bf x},\cdot)$ is a $1$-strongly convex global upper bound on the objective $F$. Computing this upper bound amounts to performing partial minimization over ${\bf z}$ in (2). Minimizing the upper bound $h({\bf x},{\bf y})$ in ${\bf y}$ corresponds to partial minimization over ${\bf x}$ in (2). The upper bound is tight in the sense that we recover the value of the objective at the current point, $h({\bf x},{\bf x})=F({\bf x})$. Therefore the alternating minimization algorithm is actually a majorization-minimization method for the nonsmooth least-squares problem

[TABLE]

The steps presented in (III-A) can then be summarized as follows:

[TABLE]

which is exactly the mapping given in (III-A).

III-C Projected Gradient Descent Interpretation

We now provide an additional interpretation of the alternating minimization algorithm as a projected gradient method for an optimization problem related to (2) which consists of a smooth convex objective and a nonconvex constraint set. This interpretation is valid whenever $g$ is assumed to be proper lower semicontinuous and convex.

For any ${\bf x}\in\mathbb{R}^{n}$ and ${\bf z}\in\mathbb{C}^{n}$, we can write (2) as $\|{\bf x}-{\bf z}\|_{2}^{2}=\|{\bf x}-{\rm Re}({\bf z})\|_{2}^{2}+\|{\rm Im}({\bf z})\|_{2}^{2}$. To move from complex numbers to real numbers, we set ${\bf w}_{1}={\rm Re}({\bf z})$ and ${\bf w}_{2}={\rm Im}({\bf z})$. Defining a new constraint set $\tilde{Z}_{{\bf c}}=\left\{({\bf w}_{1},{\bf w}_{2})\in\mathbb{R}^{n}\times\mathbb{R}^{n}:\;{\bf w}_{1}+i{\bf w}_{2}\in Z_{{\bf c}}\right\}$, problem (2) can be equivalently rewritten in the form

[TABLE]

Minimizing first w.r.t. ${\bf x}$, (13) reduces to the following minimization problem in ${\bf w}_{1},{\bf w}_{2}$:

[TABLE]

where

[TABLE]

The following result allows us to relate the gradient of $H$ to the optimization primitives used in the alternating minimization method.

Lemma III.3

Assume that $g$ is proper, lower semicontinuous and convex. Then the function $H$ is continuously differentiable, its gradient is $1$-Lipschitz and can be expressed as

[TABLE]

Proof:

From Moreau [25, Proposition 7.d], we know that $G$ is differentiable and $\nabla G({\bf x})={\bf x}-\operatorname{prox}_{g}({\bf x})=\operatorname{prox}_{g^{*}}({\bf x})$, where $g^{*}$ is the conjugate function of $g$, which is convex. The computation of the gradient of $H$ is then immediate. We can use the fact that proximity operators of convex functions are nonexpansive [25, Proposition 5.b] to verify that $\nabla H$ is $1$-Lipschitz. Indeed, for any $({\bf w}_{1},{\bf w}_{2})$ and $(\tilde{{\bf w}}_{1},\tilde{{\bf w}}_{2})$, we have

[TABLE]

completing the proof. ∎

Consider applying projected gradient descent to solve (14). From Lemma III.3, we can use a step size of magnitude $1$. In this case, taking into account the form of the gradient given in (15), we obtain that the general update step takes the form

[TABLE]

We now go back to the complex domain by setting ${\bf z}={\bf w}_{1}+i{\bf w}_{2}$. Note that projecting $({\bf w}_{1},{\bf w}_{2})$ onto $\tilde{Z}_{{\bf c}}$ is equivalent to projecting ${\bf z}$ onto $Z_{{\bf c}}$. With this notation, the iterations of projected gradient descent can be summarized by the following iteration mapping (on complex numbers):

[TABLE]

which is exactly the same as (III-A). Therefore, the Fienup algorithm is equivalent to projected gradient descent with unit stepsize applied to the formulation (14). Note that from the point of view of nonsmooth analysis, problem (14) is much better behaved than (12).

IV Consequences and Extensions

The interpretations of Section III can be used to analyze the convergence of alternating minimization applied to problem (2) and to offer extensions of the method.

IV-A Convergence Analysis

Our main convergence result is given in the following theorem. Recall that a function is semi-algebraic if its graph can be defined by combining systems of polynomial equalities and inequalities (for example, the $\ell_{1}$ norm is semi-algebraic).

Theorem IV.1

Assume that $g$ is proper, lower semicontinuous, convex and semi-algebraic. Then the sequence $\{{\bf x}^{k},{\bf z}^{k}\}_{k\in\mathbb{N}}$ generated by the alternating minimization algorithm satisfies the following:

• (i)

It holds that $\sum_{k\geq 0}\|{\bf x}^{k+1}-{\bf x}^{k}\|_{2}<+\infty$ and the sequence $\{{\bf x}^{k}\}_{k\in\mathbb{N}}$ converges to a point ${\bf x}^{*}\in\mathbb{R}^{n}$.

• (ii)

For any accumulation point ${\bf z}^{*}$ of $\{{\bf z}^{k}\}_{k\in\mathbb{N}}$, $({\bf x}^{*},{\bf z}^{*})$ is a Fréchet critical point of problem (2) and $({\bf w}_{1}^{*},{\bf w}_{2}^{*})=({\rm Re}({\bf z}^{*}),{\rm Im}({\bf z}^{*}))$ is a Fréchet critical point of problem (14).

The proof is quite technical and is given in Appendix A. The semi-algebraic assumption on $g$ can be relaxed to representability in o-minimal structures over the real field, see [29]. Therefore, the proposed result actually applies to much more general regularizers. For example, using boundedness of the feasible set in (14), the same result holds if $g$ is analytic (see the dicussion in [30, Section 5]). The arguments build upon a nonsmooth variant of the celebrated Kurdyka-Łojasiewicz (KL) property [21, 22, 23, 24]. Note that direct application of the results of [19, 20] to projected gradient descent or the results of [30, 31] to majorization minimization is not possible here.

The most important implication of Theorem IV.1 is that the sequence of estimated signals converges smoothly to a point which satisfies certain optimality conditions related to problems (2), and (14). This is a departure from standard convergence results that are only able to guarantee that accumulation points of the generated sequence of iterates satisfy certain optimality conditions. It is important to underline that the result is global: it holds for any initialization of the algorithm and does not require any regularity assumption beyond semi-algebraicity and convexity of $g$. This is in contrast with local convergence results which are typical for alternating projection methods [32, 33] that are applicable when the prior term $g$ is an indicator function.

IV-B Acceleration and Momentum Term

A benefit of the interpretation of alternating minimization as a projected gradient method is that it allows to propose new variants inspired by known extensions for projected gradient algorithms. In this section we focus on the incorporation of an inertial term that results in an alternating minimization scheme that includes a momentum term. This line of research has a long history in optimization, starting with the development of the heavy-ball method [34] which inspired an optimal first order scheme for convex optimization developed by Nesterov [35], and its extension to convex composite problems with the FISTA method [18]. Although this last technique was proposed and analyzed only in the context of convex optimization, we consider its application in our nonconvex constrained problem since it empirically provides interesting results. The resulting algorithm is referred to as FISTAPH, and is described as follows.

FISTAPH: FISTA for Phase retrieval

Initialization. ${\bf z}^{0}\in Z_{{\bf c}}$ and $\alpha^{k}\in\left[0,1\right)$ for all $k\in\mathbb{N}$. Set ${\bf y}^{0}={\bf z}^{0}$ and ${\bf z}^{-1}={\bf z}^{0}$.

General Step. For $k\in\mathbb{N}$,

•

${\bf z}^{k+1}\in P_{Z_{{\bf c}}}(\operatorname{prox}_{g}({\rm Re}({\bf y}^{k})))$.

•

${\bf y}^{k+1}={\bf z}^{k}+\alpha^{k}\left({\bf z}^{k}-{\bf z}^{k-1}\right)$.

If ${\bf z}^{m}$ is the last produced iteration, then the output of the algorithm is $\hat{{\bf x}}=\operatorname{prox}_{g}({\rm Re}({\bf z}^{m}))$. A typical choice for the weight sequence is $\alpha^{k}=\frac{k-1}{k+2}$. The question of the convergence of the iterates produced by this method in nonconvex settings is an interesting topic to explore in future research. We may also further consider monotone variants of similar types of methods, see e.g. [26].

In the numerical experiments we employ FISTAPH in the setting where $g({\bf x})=\lambda\|{\bf x}\|_{1}$ for some $\lambda>0$. In this case, $\operatorname{prox}_{g}={\mathcal{T}}_{\lambda}$ with ${\mathcal{T}}_{\lambda}$ being the soft thresholding operator with parameter $\lambda$ (see footnote on page 1).

V Experiments and Numerical Results

In this section, we describe experiments and numerical results comparing the different algorithms introduced in Section II-C on the task of phase retrieval.

V-A Experimental Setup

Given measurements ${\bf c}$ as in (1), our problem consists of finding the corresponding ${\bf x}_{0}$. We focus on the setting in which ${\bf x}_{0}$ is known to be sparse. We vary the signal size $n$ (with $J=\{1,2,\ldots,n/2\}$), the sparsity level $K$ and the signal to noise ratio (SNR). In the following discussion, we will refer to a recovery method $\mathcal{M}$ which can be seen as a black box which takes as input a vector of measurements ${\bf c}\in\mathbb{R}_{+}^{n}$, support information $J$, sparsity level $K$, an initial estimate ${\bf x}$ and outputs an estimate $\hat{{\bf x}}\in\mathbb{R}^{n}$ with ${\rm supp}(\hat{{\bf x}})\subseteq J$ and $\|{\bf x}\|_{0}\leq K$. One recovery experiment consists of the following:

• •

Fix a recovery method $\mathcal{M}$, a signal length $n$, a support information set $J=\{1,2,\ldots,n/2\}$, a sparsity level $K$ and an SNR.

• •

Generate ${\bf x}_{0}\in\mathbb{R}^{n}$ by the following procedure:

• –

Choose $K$ coordinates among $J$ uniformly at random.

• –

Set these coordinate values at random in $[-4,-3]\cup[3,4]$.

• –

Set all other coordinates to be [math].

• •

Generate the measurements ${\bf c}^{2}=|\mathcal{F}({\bf x}_{0})|^{2}+\boldsymbol{\epsilon}$, where $\boldsymbol{\epsilon}$ is white Gaussian noise according to the chosen SNR. Set negative entries of ${\bf c}^{2}$ to be [math] in order to take square root.

• •

Call method $\mathcal{M}$ 100 times with data $({\bf c},J,K)$ and randomly generate initial estimates to get 100 candidate solutions $\{\hat{{\bf x}}_{it}\}_{it=1,2,\ldots,100}$.

• •

Compute the best estimate $\hat{{\bf x}}_{best}$ with $best=\operatorname*{argmin}_{it=1,2\ldots,100}\{\||\mathcal{F}(\hat{{\bf x}}_{it})|-{\bf c}\|_{2}^{2}\}$.

• •

Compare $\hat{{\bf x}}_{best}$ and ${\bf x}_{0}$ (modulo Fourier invariances) with the following metric (sign is understood coordinatewise with ${\rm sign}(0)=0$):

[TABLE]

This procedure was repeated 100 times. That is, for each method, signal length, sparsity level and SNR, we have 100 signal recovery experiments, each one associated with a support recovery status. We aggregate these results by considering the recovery probability (average of $recovery(\hat{{\bf x}}_{best},{\bf x}_{0})$) and the median CPU usage for a single simulation (100 calls to the method with different initialization estimates). We use the same initialization for all methods by careful initialization of random seeds. All the experiments were performed on a desk station with two 3.2 GHz Quad Core Intel Xeon processors and 64GB of RAM.

V-B Implementation Details

In our numerical implementation, we used the following stopping criterion.

• •

For alternating minimization and Wirtinger methods: the difference in successive objective value less than $10^{-8}$.

• •

For GESPAR: no swap improvement.

• •

For FISTAPH: the norm of the gradient mapping less than $10^{-8}$.

The tuning of these criteria allows to balance accuracy and computational time to some extent.

The $\ell_{1}$ penalized problem includes a prior sparsity inducing term of the form $g(\cdot)=\lambda\|\cdot\|_{1}$. It is necessary to tune the $\lambda$ parameter in order to obtain meaningful results. We considered the following strategies for different methods.

• •

For alternating minimization, we set $\lambda=0.2$ in all experiments.

• •

For Wirtinger based method, $\lambda$ is tuned a posteriori as a function of $n$ and $k$. The experiment was conducted for $\lambda=1,2.15,4.64,10,21.5,46.4,100,215,464$, and we report only the best experiment for each setting.

An interesting feature of alternating minimization based methods is that in our experiments, recovery performance was very consistent for different values of $\lambda$ in different settings. As a result, we chose a single value of $\lambda$ for all experiments. The tuning of $\lambda$ for Wirtinger based algorithms is practically much more difficult. In particular, we found that the best $\lambda$ was a highly dependant function of the sparsity level $K$.

Finally, we note that $\ell_{0}$ based priors have the sparsity level of the estimate, $K$, as a parameter. On the other hand, $\ell_{1}$ based priors will not necessarily produce $K$-sparse estimates. We therefore use truncation and keep the $K$ largest entries in absolute value of the last iteration.

V-C Numerical Results

The performance in terms of support recovery are presented in Figure 1 with the corresponding algorithm run time in Figure 2. Each point in these plots is an average over 100 simulations of the recovery process, each simulation consisting of 100 random initializations of the method considered. AM corresponds to Fienup methods with different priors, FISTA is the accelerated variant, and WIRT stands for Wirtinger.

We make the following observations from the numerical results:

• •

For alternating minimization, there is a consistent increase in recovery performance by switching from $\ell_{0}$ to $\ell_{1}$ based regularization priors.

• •

The $\ell_{1}$ prior degrades the performance of Wirtinger based methods compared to the $\ell_{0}$ prior.

• •

FISTAPH consistently provides the best performance and is significantly faster than its competitors.

• •

Fienup with $\ell_{0}$ prior leads to lower performance compared to GESPAR, which was already reported in [4].

As described in the experimental section, we added noise on the squared measurements rather than on the measurements themselves. This noise model is closer to the optimization model considered for GESPAR and Wirtinger flow than model (2) which is related to problem (12). We tried to change the noise model on a subset of experiments (additive noise on the measurements rather than squared measurements), however, the performance of the different methods was very similar. Therefore, we only report results related to squared-measurement noise model.

VI Conclusion

The main theoretical contribution of this work is to provide a strong theoretical basis to the fact that Fienup-type methods, when used with Fourier transforms and convex priors, lead to smoothly converging sequences of estimates. This result holds under minimal assumptions and in particular, it holds globally, independently of the initialization point. Furthermore, we characterize the properties of the limiting point as Fréchet critical points of different optimization problems. These results shed light on important properties of one of the most well known algorithms used in the context of phase retrieval. Furthermore, based on an interpretation as a projected gradient method, we proposed a new variant of Fienup with the incorporation of a momentum term which we call FISTAPH.

On the practical side, we demonstrated based on numerical simulations that FISTAPH with $\ell_{1}$ regularization constitutes a very competitive alternative to other methods in the context of sparse phase retrieval.

Appendix A Proof of Theorem IV.1

The proof involves many notions of nonsmooth analysis which can be found in [27]. Throughout the proof, we only consider subgradients of subdifferentially regular functions. Each subgradient can be interpreted as a Fréchet subgradient and the subgradient set valued mapping is closed. We adopt the notation of Section III-C, letting ${\bf z}={\bf w}_{1}+i{\bf w}_{2}$ for two real vectors ${\bf w}_{1}$ and ${\bf w}_{2}$ and consider the constraint set $\tilde{Z}_{{\bf c}}=\left\{({\bf w}_{1},{\bf w}_{2})\in\mathbb{R}^{n}\times\mathbb{R}^{n};\;{\bf w}_{1}+i{\bf w}_{2}\in Z_{{\bf c}}\right\}$. We let $K({\bf x},{\bf w}_{1},{\bf w}_{2})=\frac{1}{2}\|{\bf x}-{\bf w}_{1}\|_{2}^{2}+\frac{1}{2}\|{\bf w}_{2}\|_{2}^{2}+g({\bf x})$ be the objective function of problem (2) which with this notation becomes

[TABLE]

We will denote by $\delta_{\tilde{Z}_{{\bf c}}}$, the indicator function of the set $\tilde{Z}_{{\bf c}}$ ([math] on the set and $+\infty$ outside). We set $\tilde{K}({\bf x},{\bf w}_{1},{\bf w}_{2})=K({\bf x},{\bf w}_{1},{\bf w}_{2})+\delta_{\tilde{Z}_{{\bf c}}}({\bf w}_{1},{\bf w}_{2})$ so that problem (16) is equivalent to the (unconstrained) minimization of $\tilde{K}$.

Proof of $(i)$: Using [27, Proposition 10.5 and Exercise 10.10], the subgradient of this nonsmooth function is of the form

[TABLE]

Partial minimization over iterations yields the following

[TABLE]

Combining these, we have

[TABLE]

Using (28),

[TABLE]

Finally, from strong convexity of $\tilde{K}$ with respect to its first argument, we have

[TABLE]

Since $g$ is semi-algebraic, $\tilde{K}$ is also semi-algebraic. Any semi-algebraic function satisfies the nonsmooth Kurdyka-Łojasievicz property [24]. We can now use the now well established recipe [19, Section 2.3] [20, Section 3.2] with the two conditions (37) and (A) to obtain that the sequence $\left\{\|{\bf x}^{k+1}-{\bf x}^{k}\|_{2}\right\}_{k\in\mathbb{N}}$ is summable. This proves statement $(i)$ (convergence holds by Cauchy criterion).

Proof of $(ii)$: Using the fact that $\tilde{K}$ has compact sublevel sets, the sequence $\left\{({\bf x}^{k+1},{\bf w}_{1}^{k},{\bf w}_{2}^{k})\right\}_{k\in\mathbb{N}}$ is bounded and hence has a converging subsequence. We fix an accumulation point $({\bf x}^{*},{\bf w}_{1}^{*},{\bf w}_{2}^{*})$ of the sequence (note that ${\bf x}^{*}$ is given by $(i)$). We remark that, thanks to (33) and the fact that $\|{\bf x}^{k+1}-{\bf x}^{k}\|\to 0$, any accumulation point of the sequence is a critical point of $\tilde{K}$. Furthermore, since ${\bf x}^{k}\to{\bf x}^{*}$, we have using (33) that

[TABLE]

This is actually the criticality condition for problem (14) which proves statement $(ii)$.