1. INTRODUCTION

Recovery of a signal from several measured intensity patterns, also known as the phase retrieval problem, is of great interest in optics and imaging. Recently, it was shown in [1] that the problem of estimating the wavefront aberration from measurements of the point spread functions can be formulated as a phase retrieval problem.

In this paper, we consider the general phase retrieval problem [2],

The sparse variant of the phase retrieval problem corresponds to the case that the unknown parameter $\mathbf{a}$ is a sparse vector. A special case of this problem is when the measurements are the magnitude of the Fourier transform of multiples of $\mathbf{a}$ with certain phase diversity patterns. A number of algorithms utilizing the Fourier transform have been proposed for solving this class of phase retrieval problems [3–5].

The fundamental nature of Eq. (1) has given rise to a wide variety of solution methods that have been developed for specific variants of this problem since the observation of Sayre in 1952 that phase information of a scattered wave may be recovered from the recorded intensity patterns at and between Bragg peaks of a diffracted wave [6]. Direct methods [7] usually use insights about the crystallographic structure and randomization to search for the missing phase information. The requirement of such a priori structural information and the expensive computational complexity often limit the application of these methods in practice.

A second class of methods first devised by Gerchberg and Saxton [8] and Fienup [3] can be described as variants of the method of alternating projections on certain sets defined by the constraints. For an overview of these methods and latter refinements, we refer the reader to [4,9].

In [10], Eq. (1) is relaxed to a convex optimization problem. The inclusion of the sparsity constraint in the same framework of convex relaxations has been considered in [11]. However, as reported in [5], the combination of matrix lifting and semidefinite programming (SDP) makes this method not suitable for large-scale problems. To deal with large-scale problems, the authors of [5] have proposed an iterative solution method, called greedy sparse phase retrieval (GESPAR), which appears to yield the promising recovery of very sparse signals. However, this method consists of a heuristic search for the support of $\mathbf{a}$ in combination with a variant of the Gauss–Newton method, whose computational complexity is often expensive. These algorithmic features are potential drawbacks of GESPAR.

In this paper, we propose a sequence of convex relaxations for the phase retrieval problem in Eq. (1). Contrary to existing convex relaxation schemes such as those proposed in [10,11], matrix lifting is not required in our strategy. The obtained convex problems are affine in the unknown parameter vector $\mathbf{a}$. Contrary to [12], our strategy does not require the tuning of regularization parameters when the measurements are corrupted by noise. We then present an algorithm based on the alternating direction method of multipliers (ADMM) that can solve the resulting optimization problems effectively. This potentially addresses the restriction of current SDP-based methods to only relatively small-scale problems.

In Section 2 we formulate the estimation problem of our interest for both zonal and modal forms. In Section 3 we propose an algorithm for solving this problem. The algorithm is based on iteratively minimizing a nuclear norm. The nuclear norm of a matrix is the sum of its singular values. Its benefit in optimization is that it is used as a convex relaxation to the rank function [13]. The convexity enables direct use of standard software libraries for solving convex optimization problems. However, since it is a computationally heavy minimization problem, we suggest an ADMM-based algorithm in Section 4 that exploits the problem structure and is therefore more efficient in practical cases. This ADMM algorithm features two minimization problems whose solutions can be computed exactly and with complexity $\mathcal{O}(n_y{n}_a)$, where $n_y$ is the number of measurements, and $n_a$ is the number of unknown variables. To find these solutions, either a least-squares problem has to be solved or the singular value decompositions of $2\times 2$ matrices have to be computed. Analytic solutions for the ADMM algorithm update steps will be presented in Subsections A and B. The convergence behavior of the algorithm proposed in Section 3 is analyzed in Section 5. Compared to the other sections, the mathematical analysis in this section is more involved, which is often the case for convergence analyses. In Section 6 we describe and discuss the results of a number of numerical experiments that demonstrate the promising performances of our algorithms. We end with concluding remarks in Section 7.

2. WAVEFRONT ESTIMATION FROM INTENSITY MEASUREMENTS

The problem of phase retrieval from the point spread function images can be approached from two directions. We take the opportunity to present them in a unified way. We first describe the problem in zonal form and then in modal form. The modal form approach used in this paper seems less popular than that of the zonal form.

A. Problem Formulation in Zonal Form

In [1] it was shown that reconstructing the wavefront from charge-coupled device (CCD) recorded images of a point source may also be formulated as a phase retrieval problem. These recorded images are called point spread functions (PSFs). As such approaches avoid the requirement of extra hardware to sense the wavefront, such as a Shack–Hartmann wavefront sensor, the problem is relevant and summarized here.

The PSF is derived from the magnitude of the Fourier transform of the generalized pupil function (GPF). For an aberrated optical system, the GPF is defined as the complex-valued function [14],

The aim of the wavefront reconstruction problem is to estimate $\varphi (\rho ,\theta )$. Once this phase aberration of an optical system has been estimated, it can be corrected by using phase modulating devices such as deformable mirrors.

In order to estimate $\varphi (\rho ,\theta )$, a known phase diversity pattern ${\varphi }_d(\rho ,\theta )$ can be introduced (e.g., by using a deformable mirror) to transform the GPF in a controlled manner into the aberrated GPF,

Let us define $\mathrm{vect}(·)$, the vectorization operator, such that $\mathrm{vect}(Z)$ yields the vector obtained by stacking the columns of matrix $Z$ into a column vector. The inverse operator ${\mathrm{vect}}^{-1}(·)$, which maps a column vector of size $m^2$ to a square matrix of size $m\times m$, is also well defined. Let in particular the matrix $Z$ and the vector $\mathbf{a}$ be defined as

By stacking the vectors ${\mathbf{y}}_d$ and the matrices $U_d$ obtained from the $n_d$ images with $n_d$ different phase diversities, correspondingly, into the vector $\mathbf{y}$ and the matrix $U$ (of size $n_d{m}^2\times m^2$), the problem of finding $\mathbf{a}$ from noise-free intensity measurements can be formulated as in Eq. (1), and that from noisy measurements can be formulated as in Eq. (2) for $n_a=m^2$ and $n_y=n_d{m}^2$.

It is worth noting that the dimension of the unknown $\mathbf{a}$ with $m$ in the range of a couple of hundreds turns this problem into a non-convex large-scale optimization problem. For such a problem, the implementation of PhaseLift [12] using standard SDP using libraries like MOSEK [15] will not be tractable because of the large matrix dimensions of the unknown quantity. If we assume that the computational complexity of semidefinite programming with matrix constraints of size $n\times n$ increases with $\mathcal{O}(n^6)$ [16], then a naive implementation of the PhaseLift method applied to Eq. (2) involving a single image has a worst-case computational complexity of $\mathcal{O}(m^{12})$.

B. Problem Formulation in Modal Form

In general, only approximate solutions can be expected for a phase retrieval problem. In the modal form of the phase retrieval problem, also considered in [1] for extended Nijboer–Zernike (ENZ) basis functions, the GPF is assumed to be well approximated by a weighted sum of basis functions. We make use of real-valued radial basis functions [17] with complex coefficients to approximate the GPF. These are studied in the scope of wavefront estimation in [18], and an illustration of these basis function on a $4\times 4$ grid in the pupil plane is given in Fig. 1.

 

Fig. 1. 16 radial basis functions with centers in a $4\times 4$ grid, with circular aperture support.

Download Full Size | PDF

Switching from the polar coordinates $(\rho ,\theta )$ to the Cartesian coordinates $(x,y)$ in the pupil plane, let us consider the radial basis functions and the approximate GPF given by

The aberrated GPF corresponding to the introduction of phase diversity ${\varphi }_d$ is

The normalized complex PSF is the two-dimensional Fourier transform of the GPF [19,20]. The aberrated PSF corresponding to the aberrated GPF in Eq. (8) is given as

We now drop the dependency on the coordinates and vectorize expression Eq. (9) for all $n_d$ diversities that have been applied to obtain the following compact form of a single matrix-vector multiplication,

Let the intensity of the PSFs be recorded on the corresponding grid of pixels of size $m\times m$, and let the vectorization of this intensity pattern for different phase diversities be concatenated into the vector $\mathbf{y}$. We can again formulate the problem of finding $\mathbf{a}$ from noise-free intensity measurements as in Eq. (1) and from noisy measurements as in Eq. (2) for $n_y=m^2{n}_d$.

It is worth noting that the dimension of $\mathbf{a}$ is not dependent on the size of the sample grid (the size of the problem). This is the fundamental advantage of the modal form formulation over the zonal form one, for which the size of $\mathbf{a}$ directly depends on the size of the problem, i.e., $n_a=m^2$.

In this paper two steps are combined to deal with the large-scale nature of optimization Eq. (2):

In the following we assume that the problem is normalized such that all entries of $\mathbf{y}$ have values between 0 and 1.

3. CONVEX OPTIMIZATION-BASED PHASE RETRIEVAL ALGORITHM

Equation (1) is equivalent to a rank constraint. Define the matrix-valued function

Our proposed algorithm in this paper relies on the following fundamental result.

Lemma 1. [21] For any $\mathbf{b}\in {\mathbb{C}}^{n_a}$, the constraint $\mathbf{y}={|U\mathbf{a}|}^2$ is equivalent to the constraint

For addressing problem Eq. (2), Lemma 1 suggests a consideration of the following approximate problem for a user-selected parameter vector $\mathbf{b}$:

Since Eq. (12) is a non-convex problem, and to anticipate the presence of measurement noise, we propose to solve the following convex optimization problem:

Note that for $\mathbf{b}=-\mathbf{a}$,

Since the result of the optimization in Eq. (13) might not produce a desired solution sufficiently fitting the measurements, we propose the iterative convex optimization-based phase retrieval (COPR) algorithm, outlined in Algorithm 1.

Algorithm 1. Convex Optimization-Based Phase Retrieval (COPR)

View Table | View all tables in this article

The nuclear norm is a convex function, and standard software like YALMIP [23] or CVX [24] can be used to concisely implement Algorithm 1. However, the nuclear norm minimization in Algorithm 1 is the main computational burden for an implementation. Usual implementations of the nuclear norm involve semidefinite constraints and require a semidefinite optimization solver. If we assume that their computational complexity increases with $\mathcal{O}(n^6)$ [16] with constraint on matrices of size $n\times n$, then minimizing the nuclear norm of the matrix $M(U,\mathbf{a},\mathbf{b},\mathbf{y})$ of size $2n_y\times 2n_y$ is computationally infeasible, even for relatively small-scale problems. Therefore, we propose a tailored ADMM algorithm of which the computational complexity of the iterations scales $\mathcal{O}(n_y{n}_a)$ and which requires the inverse of a matrix of size $2n_a\times 2n_a$ for every iteration of Algorithm 1.

4. EFFICIENT COMPUTATION OF THE SOLUTION TO EQ. (13)

The minimization problem Eq. (13) can be reformulated as

Applying the ADMM optimization technique [25,26] to the constraint optimization problem Eq. (16), we obtain the steps in Algorithm 2.

Algorithm 2. ADMM Algorithm for Solving Eq. (16)

View Table | View all tables in this article

The advantage of using this ADMM formulation is that both of the update steps in Eqs. (17) and (18) have solutions that can be computed analytically. The efficient computation of the solutions is described in the following two subsections.

A. Efficient Computation of the Solution to Eq. (17)

Upon inspection of Eq. (17), we see that this is a complex-valued standard least-squares problem, since $M(U,\mathbf{a},\mathbf{b},\mathbf{y})$ is parameterized affinely in $\mathbf{a}$. Let $\mathcal{R}(·)$ and $\mathcal{I}(·)$ denote the real and the imaginary parts of a complex object, respectively. Let the subscripts ${(·)}_1$, ${(·)}_2$, and ${(·)}_3$ denote the top-left, top-right, and bottom-left submatrices, respectively, according to Eq. (11). Define

Reordering the elements in Eq. (17), separating the real and the imaginary parts, removing all matrix elements in the argument of the Frobenius norm that do not depend on $\mathbf{a}$, and vectorizing the result gives the following least-squares problem:

Once this inverse matrix is obtained, the optimal solution to the least-squares problem in Eq. (19) can be computed by a simple matrix-vector multiplication, whose complexity scales with $\mathcal{O}(n_y{n}_a)$.

Note that in the case that the objective term includes regularization as in Eq. (14), the optimization in Eq. (19) should be modified appropriately to include the additive regularization term $\lambda {\Vert \mathbf{a}\Vert }_1$.

B. Efficient Computation of the Solution to Eq. (18)

The optimization in Eq. (18) is of the form

Lemma 2. The solution $\mathit{X}$ to Eq. (21) has singular vectors $U_C$ and $V_C$.

proof. Let $X=U_X{\mathrm{\Sigma }}_X{V}_X^T$ be a singular value decomposition of $X$. Then

Denote the singular values of $C$ in descending order as ${\sigma }_{C,1},\dots ,{\sigma }_{C,2n_y}$, and those of $X$ similarly. Thanks to Lemma 2, Eq. (21) can be simplified to

This problem is completely decoupled in ${\sigma }_{X,i}$, and the optimal solution to Eq. (22) is computed with

This shows that a valid SVD can be computed very efficiently in $\mathcal{O}(1)$, that is, in theory, in a computation time independent of the number of pixels in the image, the number of images taken, or of the number of basis functions.

5. CONVERGENCE ANALYSIS OF ALGORITHM 1

Algorithm 1 can be reformulated as a Picard iteration ${\mathbf{a}}_{k+1}\in T({\mathbf{a}}_k)$, where the fixed point operator $T\text{:}{\mathbb{C}}^{n_a}\to {\mathbb{C}}^{n_a}$ is given by

Our subsequent analysis will show that the set of fixed points, $\text{Fix}T$, the set of $\mathbf{a}$s for which $\mathbf{a}=T(\mathbf{a})$, of $T$ is in general non-convex and as a result, iterations generated by $T$ cannot be Fejér monotone (Definition 5.1 of [27]) with respect to $\text{Fix}T$. That is, each new iterate is not guaranteed to be closer to all fixed points in $\text{Fix}T$. Therefore, the widely known convergence theory based on the properties of Fejér monotone operators and averaging operators is not applicable to the operator $T$ given in Eq. (23).

In this section, we make an attempt to prove the convergence of Algorithm 1, which has been observed from our numerical experiments, via a relatively newly developed convergence theory based on the theory of pointwise almost averaging operators [28]. It is worth mentioning that we are not aware of any other analysis schemes addressing the convergence of Picard iterations generated by general non-averaging fixed point operators. Our discussion consists of two stages. Based on the convergence theory developed in [28], we first formulate a convergence criterion for Algorithm 1 (Proposition 5.1) under rather abstract assumptions on the operator $T$. Due to the highly complicated structure of the nuclear norm of a general complex matrix, we are unable to verify these mathematical conditions for general matrices $U$. However, we will verify that they are well satisfied in the case that $U$ is a unitary matrix (Theorem 5.2). From the latter result, we heuristically hope that Algorithm 1 still enjoys the convergence result when the matrix $U$ is close to being unitary in a certain sense. In Section 6 we demonstrate that convergence is obtained in practice for the imaging case.

It is a common prerequisite for analyzing the local convergence of a fixed point algorithm that the set of solutions to the original problem is non-empty. That is, there exists $\mathbf{a}\in {\mathbb{C}}^{n_a}$ such that $\mathbf{y}={|U\mathbf{a}|}^2$. Before stating the convergence result, we need to verify that the fixed point set of $T$ is non-empty.

Lemma 3. The fixed point operator $T$ defined at Eq. (23) holds

The next proposition provides an abstract convergence result for Algorithm 1. $\text{Fix}T$ is supposed to be closed. In the sequel, the metric projection associated with a set $\mathrm{\Omega }$ is denoted $P_{\mathrm{\Omega }}$,

The pointwise property, instead of the standard averaged property, imposed in (i) of Propositon 5.1 allows us to deal with the intrinsic non-convexity of the fixed point set $\text{Fix}T$. The metric subregularity assumption imposed in (ii) technically ensures adequate progression of the iterates relative to the distance from the current iterate to the fixed point set. This is not only a technical assumption but also a necessary condition for local linear convergence of a fixed point algorithm, Theorem 3.12 of [29]. Condition (iii) is, on one hand, a technical assumption and becomes redundant when $S=\text{Fix}T$. On the other hand, the set $S$ allows one to exclude from the analysis possible inhomogeneous fixed points of $T$, at which the algorithm often exposes weird convergence behavior (see Example 2.1 of [28]).

The size of neighborhood $W$ appearing in Proposition 5.1 indicates the robustness of the algorithm in terms of erroneous input (the distance from the starting point to the nearest solution).

We now apply the abstract result of Proposition 5.1 to the following special but important case.

Theorem 5.2. Let $U\in {\mathbb{C}}^{n_a\times n_a}$ be unitary and ${\mathbf{a}}^{*}\in {\mathbb{C}}^{n_a}$ be such that ${|U{\mathbf{a}}^{*}|}^2=\mathbf{y}$. Then, every Picard iteration generated by Algorithm 1 ${\mathbf{a}}_{k+1}\in T({\mathbf{a}}_k)$ starting sufficiently close to ${\mathbf{a}}^{*}$ converges linearly to a point $\tilde{\mathbf{a}}\in \text{Fix}T$, satisfying ${|U\tilde{\mathbf{a}}|}^2=\mathbf{y}$.

proof. See Appendix B. □

6. NUMERICAL EXPERIMENTS

Four important numerical aspects of the COPR algorithm, including convergence, flexibility, complexity, and robustness, are tested on relevant problems. First we discuss convergence and the number of iterations of the COPR and the ADMM algorithms. Second, we demonstrate the flexibility of the convex relaxation by comparing the COPR algorithm with an added ${\ell }_1$-regularization to the PhaseLift method [12] and to the compressive sensing phase retrieval (CPRL) method in [11] on an under-determined sparse estimation problem. Then we compare the practically observed computational complexity of COPR and a naive implementation of PhaseLift [12]. Finally, we investigate the robustness of COPR relative to noise in a Monte Carlo simulation for 25 and 100 basis functions. We compare four algorithms: COPR, PhaseLift [12], a basic alternating projections method (Section 4.3 in [12]), and an averaged projections method based on [30]. We note that the latter method fundamentally employs the Fourier transform at every iteration and hence is in general not applicable for phase retrieval in the modal form.

A. Convergence

The while loops in Algorithms 1 and 2 can be run for a fixed number of iterations. Figure 2 shows four such combinations for a typical problem with five images of size $256\times 256$, of which a subset of $25\times 25$ pixels per image is used, and 64 basis functions. All cases are identically initialized with coefficients that best approximate a flat wavefront. As can be seen from the figure and the line with a square marker, only one COPR iteration is necessary here, as the ADMM algorithm slowly converges towards 0. However, stopping the ADMM algorithm after a limited number of iterations and having more than one COPR iteration can have a clear benefit, since faster convergence is achieved this way.

 

Fig. 2. Convergence plot for four different combinations of COPR iterations and ADMM iterations. Denoted in the legend are first the number of COPR iterations and then the number of ADMM iterations used to solve Eq. (16) in each COPR iteration. Markers denote each new COPR iteration.

Download Full Size | PDF

B. Application of COPR to Compressive Sensing Problems

The first problem is to estimate 16 coefficients from 8 measurements, where the optimal vector is known to be sparse.

We generate a sparse coefficient vector $\mathbf{a}$ with two randomly generated non-zero complex elements. We generate two images ($n_d=2$, $m=128$) by applying two different amounts of defocus with Zernike coefficients $-\frac{\pi }{8}$ and $\frac{\pi }{8}$, respectively. From each image we use the center $2\times 2$ pixels, resulting in a total of $n_y=8$ measurements.

The applied algorithms are the COPR algorithm, the COPR algorithm with an additional ${\ell }_1$-regularization, the PhaseLift algorithm [12], and the CPRL algorithm of [11]. The results are displayed in Fig. 3. As can be seen from the figure, COPR and PhaseLift fail to retrieve the correct solution. The CPRL method and the regularized COPR algorithm compute the correct solution.

 

Fig. 3. Absolute values of 16 estimated coefficients according to four different algorithms.

Download Full Size | PDF

C. Computational Complexity

The second problem demonstrates the trends of the required computation time when the number of estimated coefficients increases. The underlying estimation problem consists of seven images with different amounts of defocus applied as phase diversity, where each image is of size 64 by 64 pixels. A subset of 20 by 20 pixels of each image is used in the estimation. We compare in Fig. 4 the COPR algorithm to the PhaseLift algorithm, which is implemented according to the optimization problem (2.5) in [12]. For PhaseLift, the reported time is the time it takes the MOSEK solver [15] to solve the optimization problem. This does not include the time taken by YALMIP [23] to convert the problem as given to the solver-specific form. For COPR, the initial guesses for the coefficients are drawn randomly from a Gaussian distribution, the number of iterations is set beforehand according to the convergence to the correct solution, and the total time is recorded. The implementation of COPR does not exploit the parallelism referred to in Section 4.B. By convergence, we mean that the estimated vector $\widehat{\mathbf{a}}$ satisfies the tolerance criterion

The minimization over parameter $c$ ensures that the (unobservable) piston mode in the phase is canceled (Let $(\begin{array}{cc}\widehat{\mathbf{a}}& {\mathbf{a}}^{*}\end{array})=\mathrm{QR}$ be the QR decomposition. Then $\angle c^{*}=\angle \frac{R_{12}}{R_{11}}$). The computational complexity of PhaseLift is, as implemented, approximately $\mathcal{O}(n^4)$. The MOSEK solver ran into numerical issues for more than 25 estimated parameters. The COPR algorithm’s computational complexity is approximately $\mathcal{O}(n)$. The better complexity is offset by a longer computation time for very small problems.

D. Robustness to Noise

When estimating an unknown phase aberration, it is more logical to evaluate the performance of the algorithm on its ability to estimate the phase and not the coefficients of the basis functions.

We assume the phase is randomly generated with a deformable mirror. Let $H\in {\mathbb{R}}^{m^2\times n_u}$ be the mirror’s influence matrix and $\mathbf{u}\in {\mathbb{R}}^{n_u}$ be the input to the mirror’s actuators, such that

The input values $u_i$ are drawn from the uniform distribution between 0 and 1. The mirror has $n_u=44$ actuators, and the images have sides $m=128$. The aperture radius is 0.4.

Five different defocus diversities are applied with Zernike coefficients uniformly spaced between $-\frac{\pi }{2}$ and $\frac{\pi }{2}$. Gaussian noise is added to the obtained images such that

The phase is estimated from $\mathbf{y}$ using four different algorithms. The first is the COPR algorithm. The second is the averaged projections (AvP) algorithm [30]. The AvP algorithm is an extension of the well-known Gerchberg–Saxton algorithm [31] for solving problems with multiple images and is in the same class of algorithms as the hybrid-input-output algorithm and the difference map algorithm [32,33]. This makes this algorithm relevant for comparison. The third is the alternating projections (AlP) method ([12], Section 4.3), and the fourth algorithm is the PhaseLift method [12].

The COPR, PhaseLift, and AlP methods are applied to estimate the phase using 25 basis functions, where the initial guesses for the coefficients are those coefficients that best approximate a flat wavefront. The AvP method is not based on the use of basis functions but on projection and the Fourier transform.

We make use of the Strehl ratio as a measure of optical quality. The Strehl ratio $S$ is the ratio of the maximum intensity of the aberrated PSF and that of the unaberrated one and can be approximated with the expression of Mahajan,

For every noise level, 100 different phases were generated with the deformable mirror model Eq. (27). The results are presented in Fig. 5. The resulting Strehl ratios are plotted with a trend line connecting the 50% quantiles. Figure 6 gives a qualitative comparison of the estimates for a single case.

 

Fig. 4. Computation time comparison between PhaseLift and COPR for different numbers of coefficients.

Download Full Size | PDF

 

Fig. 5. Strehl ratio of the estimated phase aberration as a function of SNR.

Download Full Size | PDF

 

Fig. 6. Example PSF and phase estimates of the COPR, Alternating Projection [12], Averaged Projection [30], and PhaseLift [12] algorithms for three PSF measurements with an SNR of approximately 36 dB.

Download Full Size | PDF

In the case of PhaseLift, the tuning parameter that trades off measurement fit and the rank of the “lifted” matrix is tuned once and applied to all problems. This has the effect that the reported performance is not as high as it could be with optimal tuning for individual problems. This points to another advantage of COPR: the absence of tuning parameters aside from the choice of basis functions.

The two figures show that COPR is robust to noise and gives accurate phase estimates for a wide range of noise levels.

7. CONCLUDING REMARKS

The convex relaxations in solving the phase retrieval problem as proposed in Eq. (13) have the advantage over current convex relaxation methods, such as PhaseLift, that our strategy is affine in the coefficients that are to be estimated. This allows for easy extension of the proposed method to phase retrieval problems that incorporates prior knowledge on the coefficients by regularization of the objective function. One such successful extension is the regularization with the ${\ell }_1$-norm to find sparse solutions, as demonstrated in Fig. 3.

In Section 4, an ADMM algorithm was proposed for efficient computation of the solution to Eq. (13). The result is that for the COPR algorithm a better computational complexity is observed compared to that of PhaseLift; see Fig. 4. COPR is also able to solve phase estimation problems with larger numbers of parameters.

The required computations are favorable both in computation time and accuracy (they have simple analytic solutions) and in worst-case scaling behavior $\mathcal{O}(n_y{n}_a)$ for every ADMM iteration, where $n_y$ is the number of pixels, and $n_a$ is the number of basis functions.

We discussed the convergence properties of the COPR algorithm in Section 5 and showed that for selected problems this convergence is linear or faster.

Finally, COPR has been shown to be robust against measurement noise and outperform the two projection-based methods whose naive forms are often sensitive to noise, as expected.

We are aware that in practice the performance of projection methods can be substantially better than what we have observed in this study, provided that appropriate denoising techniques are also applied. Keeping aside from the matter of using denoising techniques, we have chosen to compare the algorithms in their very definition forms.