PhasePack: A Phase Retrieval Library

Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch, Christoph Studer, Tom Goldstein

I Introduction

Phase Retrieval is an example of a non-convex quadratic program with quadratic constraints. In the real-valued case, it is a combinatorial problem of determining the missing signs of $\mathbf{A}\mathbf{x}$ , which is known to be NP-hard . Despite this observation, recent years have seen the development of new algorithms that solve phase retrieval problems effectively. Unfortunately, because of the lack of publicly available real-world data, the lack of a common software interface for different algorithms, and a knowledge gap between practitioners and theoreticians, only little work has been devoted to compare and evaluate newer phase retrieval methods.

II Introducing PhasePack

PhasePack is a software package that contains implementations of many different phase retrieval methods, and tools for easily applying them to real and synthetic datasets. The purpose of PhasePack is to create a common interface for a wide range of phase retrieval schemes, and to provide a common testbed using both synthetic data and empirical imaging datasets . PhasePack can also benchmark different algorithms against one another, and generate performance comparisons with varying numbers of measurements, signal-to-noise ratio, iterations, and running time. The package handles single method testing as well as multiple method comparisons.

The methods in PhasePack differ in numerous ways from their original descriptions in the literature in order to achieve improved robustness to different measurement models and to enable faster convergence. After reviewing the kinds of methods available in PhasePack in Section III, we will discuss practical implementation consideration in Section IV.

III Methods available in PhasePack

Second, many approaches attack the non-convex formulation directly using least-squares formulations . These methods use gradient descent to minimize an objective of the form

where $\hat{}\mathbf{A}$ and $\hat{}\mathbf{b}$ contain either a subset of the measurements or a re-weighting of the system $\mathbf{A}\mathbf{x}=\mathbf{b},$ and $p\in\{1,2\}$ is an integer exponent. Because of the non-convexity of the formulation (1), these methods require careful initialization to avoid local minimizers.

The third class of methods convert phase retrieval into a convex problem. This includes the “lifting methods” PhaseLift and PhaseCut , which square the dimensionality of the problem. PhasePack also implements several low-rank approximate solvers for large instances of PhaseLift, including gauge duality methods and sketching methods . Finally, PhasePack contains the newer non-lifting relaxations PhaseMax and PhaseLamp .

IV Building practical implementations

Many recent phase retrieval methods have been developed with the goal of proving rigorous guarantees for random Gaussian measurement matrices. As a result, the implementations described in the literature are not optimized for performance, and may become unstable when used on non-Gaussian measurements. This is particularly an issue for variations of Wirtinger flow, which require the choice of stepsize parameters. This is illustrated in Figure 1, which compares a “strict” implementation of Wirtinger flow using the stepsize rules described in to PhasePack’s implementation with adaptive stepsize.

The implementations in PhasePack depart from the literature in several ways that make them more robust and efficient. We discuss several such issues below.

Wirtinger Flow and its variants rely on gradient descent methods, and are sensitive to the choice of stepsize parameters. Existing stepsize rules presented in the literature are designed for Gaussian measurement models, and may be unstable or slow when used on generic matrices. PhasePack solves least-squares formulations of the phase retrieval problem using the general gradient descent solver FASTA , which supports adaptive stepsizes, automated stopping conditions, and conjugate-gradient acceleration.

This gradient solver automates optimization in several ways. For stepsize selection, PhasePack uses the Barzilai-Borwein adaptive method . To guarantee stability, the method uses a backtracking line-search similar to the classical Armijo line search . Classical line searches enforce that the objective decreases monotonically. For non-convex problems where the local curvature changes rapidly and local minima are prevalent, monotonic searches may result in excessive backtracking, and may be more prone to getting stuck in bad local minimizers. For this reason, several authors have proposed non-monotonic line searches . Rather than requiring $f(\mathbf{x}^{(k+1)})<f(\mathbf{x}^{(k)})$ on the $k$ th iteration, these methods require the weaker condition $f(\mathbf{x}^{(k+1)})<\max_{k-w<j<k}\{f(\mathbf{x}^{(j)})\}$ for some “window” parameter $w.$ The particular adaptive strategy and line search rule used by PhasePack is described in detail in .

Finally, the solver in PhasePack is capable of using L-BFGS acceleration (which was studied for phase retrieval in ) and non-linear conjugate gradient methods to speed up convergence. However, these options need to be tuned manually.

IV-B Practical methods for spectral initialization

Non-convex methods for phase retrieval are prone to getting stuck in local minimizers if they are not initialized properly. Phasepack includes two main classes of initializers: (i) the spectral initializer and its variants and, (ii) the orthogonality promoting initializer. The initializers described in the literature are optimized for random Gaussian data, and may not perform well on non-Gaussian real-world data. To be robust to different measurement models, PhasePack initializers depart from the literature in three main ways: we use a Krylov subspace-based eigensolver, we re-scale data before pre-processing, and we re-scale the initialization vector after computation. We discuss these important differences below after briefly describing how initialization methods work.

Spectral initialization methods begin by forming the matrix

where $\mathbf{a}_{i}$ is the $i$ th row of $\mathbf{A}$ , and $T(\cdot)$ is a “pre-processing” function, which is simply the identity in the original spectral method . These methods are motivated by the observation that if $\mathbf{a}$ is a random Gaussian vector and $b$ is its corresponding measurement, then we have

In fact, the leading eigenvector of (3) is the unknown signal $\mathbf{x}.$ Spectral methods work by approximating the expectation (3) with the empirical matrix (2), computing the leading eigenvector, and then scaling it appropriately. A variety of methods have been proposed with different pre-processing functions, including the identity , a “truncation” operator that is zero for small magnitude measurements , and a “re-weighting” operator that shrinks its argument . A rigorous study of these methods appears in . The default initialization in PhaseMax uses the “optimal” preprocessing function proposed in , which is

where $\delta=m/n$ is the number of measurements divided by the signal dimension.

Orthogonality promoting initializers work by identifying measurement vectors that are not correlated with the signal $\mathbf{x}$ , and finding an initialization vector that is orthogonal to those uncorrelated vectors. This is done by forming the matrix

where $\mathcal{I}$ contains the indices of measurement vectors that produce the smallest measurements. The initializer is then the eigenvector of $\mathbf{Y}$ with smallest eigenvalue.

PhasePack initializers depart from the literature in several ways that we now describe. First the literature relies on the power method to compute leading eigenvectors. This is limiting for several reasons. The power method can be unstable if the matrix does not have a large spectral gap (which empirical matrices may have not), and it is incapable of finding the smallest eigenvector (which is needed for orthogonality promoting initializers). For this reason, we compute eigenvectors using the iterated Arnoldi method with restart, which is the default algorithm behind MATLAB’s “eigs” routine. This method exhibits faster convergence than the power method and can reliably resolve both the leading and trailing eigenvectors.

A more major departure from the literature is in how PhasePack scales data before pre-processing. The pre-processing functions $T$ defined in the literature are optimized for random Gaussian measurement matrices and normalized signals, and initializers may be highly sub-optimal (or even ill-posed) when empirical data is used. For this reason, data is re-scaled to have the same statistics as the Gaussian model before pre-processing. Note that the measurement operator may be implemented as a black-box function, and so re-scaling needs to be done without access to the matrix entries or the underlying signal.

Finally, PhasePack concludes initialization with a least-squares step to determine the optimal length of the initial vector. Most eigensolvers produce eigenvalues that have unit length by default, however, the length of the signal may effect the convergence of the phase retrieval solver. For this reason, after finding the initializer $\hat{}\mathbf{x}$ we compute the least-squares solution to

and then replace the initializer $\hat{}\mathbf{x}$ with $\alpha\hat{}\mathbf{x}.$ The spectral initialization methods described above come with their own scaling routines, however PhasePack’s least-squares scaling method is more practical for empirical data as it does not depend on the Gaussian measurement model.

V Empirical datasets

Public datasets for phase retrieval are scarce, and available datasets are often difficult to use because of lack of documentation or unusual data formats. One of PhasePacks goals it to create a simple API for testing methods on empirical datasets. PhasePack provides routines for unpacking and pre-processing datasets, and prepares a simple measurement operator $\mathbf{A}$ and measurement vector $\mathbf{b}$ that can be used for phase retrieval. The datasets currently supported by PhasePack were obtained for the purpose of imaging an object through a diffusive medium, and are described in detail in . We provide a brief description of this imaging modality below.

The imaging modality described is for reconstructing an image from light that has passed through a diffusive medium (e.g., a light-scattering material like paper or eggshell). Imaging though such a material is impossible using a conventional camera, but this can be done using phase retrieval methods. A binary mask is created by placing a transparent LCD display in front of a light source. The pixels in this mask can be made either opaque or clear. An image pattern in loaded onto the mask, and the mask is illuminated from behind using a coherent (laser) light source. The light passes through the mask (creating an image), and then passes through a diffuser. The diffused light then lands on a photodetector that measures the intensity of light (but not the phase) at many locations, producing measurements.

The light pattern (image) passing through the mask is the unknown signal $\mathbf{x}$ . The diffuser is modeled as a “transmission matrix” $\mathbf{A}$ that describes how light passing through each pixel of the mask effects each of the detectors. Under this model, the signal received at each detector is given by $\mathbf{A}\mathbf{x}.$ Because standard photodetectors measure intensity and not phase, the recorded signal is $\mathbf{b}=|\mathbf{A}\mathbf{x}|,$ and the image must be reconstructed by phase retrieval.

This imaging modality is particularly useful for benchmarking phase retrieval methods because the ground-truth image (the pattern loaded onto the mask) is known, and can be used to evaluate the quality of reconstructions. A sample mask and its reconstructions is shown in Figure 3.

VI Experimental Results

We demonstrate the capabilities of PhasePack using a range of experiments on real and synthetic data. All of these examples were produced using the scripts included in the benchmark sub-folder of the PhasePack distribution. Figure 2 shows the performance of several algorithms for recovering a random signal from random Gaussian measurements (left) and from measurements acquired using a real empirical transmission matrix (right). Figure 3 shows the reconstruction of an image from empirical measurements obtained using an optical device .

In most of our experiments, we found that the classical Fienup and Gerchberg-Saxton methods remain to be state-of-the-art. However, the re-weighted Amplitude flow method appears to be highly competitive, if not practically identical in performance.

Convex methods appear to under-perform compared to their non-convex alternatives. The convex method PhaseMax seems to outperform the lifted convex relaxation PhaseLift (and its approximate low-rank version SketchyCGM ) on real image data. However, PhaseLift outperforms PhaseMax by a small margin on random Gaussian measurments.

VII Discussion

PhasePack was created with the goal of providing a unified framework in which researchers can investigate phase retrieval. By providing a common interface for different methods, and a simple API for testing methods on empirical datasets, we hope that PhasePack can help the community to better understand the strengths and weaknesses of different methods.

Acknowledgments

We would like to thank Ashok Veeraghavan and Chris Metzler at Rice University for making their empirical datasets available in PhasePack.