Undersampled Phase Retrieval via Majorization-Minimization

I Introduction

In this paper, we propose two simple and efficient algorithms based on the majorization-minimization framework to solve the undersampled phase retrieval problem under two different problem settings. The first setting is simple in which the unknown signal is considered to be sparse in the standard basis. The second, more general, setting considers cases where the original signals are not sparse in the standard basis (or any other bases). We share the same idea of applying the sparse coding techniques to the phase retrieval problem , inspired by the fact that a lot of image and video signals can be sparsely approximated by a linear combination of a few columns in a dictionary . Recently the authors in have shown encouraging results of exploiting sparse coding for the oversampled phase retrieval problem (DOLPHIn algorithm). In this paper, we propose efficient algorithms to recover the unknown signals with high accuracy for the undersampled phase retrieval problem by jointly designing the dictionary and the sparse codes.

Numerical methods for the undersampled phase retrieval problem.

Monotonicity and guarantee of convergence to a stationary point for the sequence of points generated by our algorithms.

Faster numerical convergence of our algorithms compared to state-of-the-art methods UPRwO and DOLPHIn.

Low complexity per iteration of our algorithms (only requiring basic matrix multiplication).

The remaining sections are organized as follows. We first provide a brief introduction on the majorization-minimization framework in section II. Later, we propose an algorithm to solve the undersampled phase retrieval problem of sparse signals using the majorization-minimization techniques. When the unknown signals are not sparse, we propose another algorithm in section III to solve the undersampled phase retrieval problem through sparse coding. Numerical results and comparisons with up-to-date benchmark methods are presented and discussed in section IV. Finally, we conclude our work in section V.

II Compressive Phase Retrieval via Majorization-Minimization

In this section, we first provide a brief overview of the general majorization-minimization (MM) framework. Later, we propose a simple iterative algorithm to solve the undersampled phase retrieval problem of sparse signals via the MM techniques.

The majorization-minimization (MM) algorithm is an iterative optimization method, which includs the well-known expectation-maximization (EM) algorithm as a special case. Instead of solving the original difficult optimization problem directly, an MM algorithm deals with a sequence of simple surrogate problems that produce a sequence of points to drive the original objective function downhill.

For a real valued function $f(\boldsymbol{\theta})$, any function $g(\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(m)})$ satisfying the following two conditions is a majorization function of $f(\boldsymbol{\theta})$ at the point $\boldsymbol{\theta}^{(m)}$:

The function $g(\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(m)})$ is a global upper bound of $f(\boldsymbol{\theta})$ and touches it at the point $\boldsymbol{\theta}^{(m)}$. In general, these majorization functions $\{g(\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(m)})\}_{m}$ are chosen to be convex and much easier to deal with than the original function, which usually is non-convex or non-differentiable.

Initialized by any feasible point $\boldsymbol{\theta}^{(0)}$, an MM algorithm generates a sequence of points $\{\boldsymbol{\theta}^{(m)}\}_{m}$ according to the updating rule:

This sequence of points $\{\boldsymbol{\theta}^{(m)}\}_{m}$ has a favorable property of driving the original objective function $f(\boldsymbol{\theta})$ downhill:

The first inequality and the third equality come from the definition of the majorization function (2). The second inequality is valid because $\boldsymbol{\theta}^{(m+1)}$ is a minimizer of $g(\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(m)})$ from (3). Therefore, one can find a stationary point for the original problem by solving the surrogate problems instead.

II-B C-PRIME

Instead of using the intensity measurements $\{y_{i}\}_{i=1}^{M}$ directly, we decide to use the modulus measurements $\{\sqrt{y_{i}}\}_{i=1}^{M}$ (we assume $y_{i}\geq 0$ otherwise we just discard this measurement). Justification on the advantage of using modulus information $\{\sqrt{y_{i}}\}_{i=1}^{M}$ over intensity information $\{y_{i}\}_{i=1}^{M}$ is provided in Appendix A. We propose to solve the following problem to balance the importance of minimizing the sum of squared error and utilizing the prior sparsity information of the original signal:

The first term is a data fitting term measuring how well the sought signal $\mathbf{x}$ fits the modulus measurements $\{\sqrt{y_{i}}\}_{i=1}^{M}$. This value should be comparable to the noise level for a successful recovery. The second term $\|\mathbf{x}\|_{1}$ is used to promote sparsity in $\mathbf{x}$. And $\rho$ is a regularization parameter to balance the weights between the sum of squared error and sparsity level to produce a desired solution.

Here the square root operator $\sqrt{\cdot}$ and the magnitude operator $|\cdot|$ are applied element-wise. This problem is not convex because of the magnitude operator. Using the majorization-minimization technique, we propose an efficient method to solve the convex surrogate problems instead. Note that

where $const.$ is a constant independent of the variable $\mathbf{x}$.

The claim is valid simply by rearranging the terms in $(\mathbf{x}-\mathbf{x}_{0})^{H}(\mathbf{M}-\mathbf{L})(\mathbf{x}-\mathbf{x}_{0})\geq 0$, cf. . ∎

According to Claim 1, the first term in (7) can be majorized as

for any constant $C\geq\lambda_{\max}(\mathbf{A}^{H}\mathbf{A})$. Further,

the second term in (7) can be majorized as

Combining these two majorization functions together, the corresponding surrogate problem for (6) is

This surrogate problem is convex in $\mathbf{x}$ and is equivalent to the following problem:

The vector $\mathbf{c}$ is a constant independent of the variable $\mathbf{x}$:

Now it is clear to see the benefits of using the majorization-minimization framework. Instead of dealing with the original non-convex non-differentiable problem (6), we only need to solve a surrogate problem (12) which has a simple closed-form solution at every iteration. We name our algorithm compressive phase retrieval via the majorization-minimization technique (C-PRIME for short) and summarize the procedure in Algorithm 1.

In the algorithm, we further adopt the SQUAREM algorithm to accelerate the convergence speed of our method. SQUAREM generally achieves a superlinear convergence rate and only requires parameter updating. Instead of updating $\mathbf{x}^{(k+1)}$ directly from $\mathbf{x}^{(k)}$ at the $k$-th iteration, SQUAREM first seeks an intermediate point $\mathbf{x}_{3}$ based on $\mathbf{x}^{(k)}$ and later updates the next point $\mathbf{x}^{(k+1)}$ from this intermediate point. Unfortunately, this updating rule may violate the descent property of the MM framework. Therefore, we add a backtracking step in our algorithm (the while loop) to maintain the descent property. In detail, we repeatedly halve the distance between $\alpha$ and $-1$ until the descent property is valid. This strategy is guaranteed to work because in the worst case where $\alpha=-1$, the intermediate point satisfies $\mathbf{x}_{3}=\mathbf{x}^{(k)}+2\mathbf{r}+\mathbf{v}=\mathbf{x}_{2}$, which ensures that the algorithm will jump out of the while loop. (Actually it only takes several updates $\alpha\leftarrow(\alpha-1)/2$ for the descent property to be maintained in the simulation.)

III Sparse Coding for Phase Retrieval

where $\mathcal{D}$ is a convex set defined as

The data fitting term in the objective is the same as in the sparse signal case discussed in the last section, the second term measures how well the unknown signal can be approximated by the dictionary, and the last term promotes sparse code so that only a few atoms are chosen to approximate the unknown signal. The two regularization parameters $\mu$ and $\rho$ are used to balance the weights on the data fitting, the dictionary representation, and the sparse code. Unfortunately, there is more than one solution for problem (15) because the unknown dictionary is considered as an additional variable. Now that one single signal is insufficient to uniquely determine the dictionary, multiple signals should be exploited to jointly recover the original signals and unknown dictionary.

The number of atoms should be less than the number of unknown signals $L<P$. Otherwise, each signal is trivially represented by a $1$-sparse vector $\mathbf{z}_{p}$ after including $\mathbf{x}_{p}/\|\mathbf{x}_{p}\|_{2}$ as an atom in the dictionary.

Problem (17) is not convex, not only because of the magnitude operator, but also because of the quadratic term $\mathbf{Dz}_{p}$. But the problem is convex with regard to $\mathbf{D}$ if $\{\mathbf{x}_{p}\}$ and $\{\mathbf{z}_{p}\}$ are fixed. Also, it is convex with regard to $\{\mathbf{z}_{p}\}$ when $\{\mathbf{x}_{p}\}$ and $\mathbf{D}$ are fixed. Another problem is that all variables are tangled together because of the shared dictionary $\mathbf{D}$. But once $\mathbf{D}$ is fixed, problem (17) can be separated into $P$ independent smaller problems. Therefore, we propose to solve this problem using the block successive upper-bound minimization method (BSUM) . BSUM is a simple and iterative algorithm framework to successively optimize upper bounds functions, instead of the original objective function, in a block by block manner. And the convergence analysis is provided in .

We first consider updating $\{\mathbf{z}_{p}^{(k+1)}\}$ at the $k+1$-th iteration. When $\{\mathbf{x}_{p}\}$ and $\mathbf{D}$ are fixed to be $\{\mathbf{x}_{p}^{(k)}\}$ and $\mathbf{D}^{(k)}$, problem (17) is equivalent to the following problem:

which can be separated into $P$ independent problems in the form:

This is a typical sparse coding problem . If the dictionary $\mathbf{D}^{(k)}$ satisfies $(\mathbf{D}^{(k)})^{H}\mathbf{D}^{(k)}=\mathbf{I}$ (semi-unitary), (19) is equivalent to the following problem:

When $(\mathbf{D}^{(k)})^{H}\mathbf{D}^{(k)}\neq\mathbf{I}$, it is difficult to find a simple closed-form solution for (19) directly. Hence we propose to solve a surrogate problem instead. Using the same technique in Claim 1, the first term in (19) can be majorized as

where $\mathbf{e}_{p}$ is independent of the variable $\mathbf{z}_{p}$ and is defined as

The scalar $E^{(k)}$ is a constant larger than $\lambda_{\max}\left((\mathbf{D}^{(k)})^{H}\mathbf{D}^{(k)}\right)$ and we prove in Appendix B that $E^{(k)}\geq L$ is sufficient for a valid majorization function. Therefore, the corresponding surrogate problem for (19) is

When $\mathbf{D}$ and $\{\mathbf{z}_{p}\}$ are fixed to be $\mathbf{D}^{(k)}$ and $\{\mathbf{z}_{p}^{(k+1)}\}$, updating $\{\mathbf{x}_{p}\}$ leads to solving the following problem:

which also can be separated into $P$ independent problems:

This problem is not convex due to the magnitude operator. Hence, we choose to solve a surrogate problem instead. According to Claim 1 and (10), choosing a constant $F\geq\lambda_{\max}(\mathbf{A}^{H}\mathbf{A})$, the objective function can be majorized as

where $\mathbf{f}_{p}$ is a constant with regard to the variable $\mathbf{x}_{p}$:

Therefore, the surrogate problem for (27) is

Note that this constant $\mathbf{f}_{p}$ is similar to the constant $\mathbf{c}$ in (14), which was used to update the signal in the last section. The additional third term $\mathbf{D}^{(k)}\mathbf{z}_{p}^{(k+1)}$ in $\mathbf{f}_{p}$ is due to the second approximation over a dictionary term in (17).

III-C Updating the Dictionary 𝐃𝐃\mathbf{D}

The final step is to update the dictionary $\mathbf{D}$ fixing the other two blocks of variables $\{\mathbf{x}_{p}\}$ and $\{\mathbf{z}_{p}\}$ as $\{\mathbf{x}_{p}^{(k+1)}\}$ and $\{\mathbf{z}_{p}^{(k+1)}\}$. Since the regularization parameter $\mu\geq 0$, we need to solve the following problem:

which can be formulated in a more compact form:

where $\mathbf{z}_{l,T}^{(k+1)}$ is a row vector denoting the $l$-th row in matrix $\mathbf{Z}^{(k+1)}$. The objective function in (34) satisfies

where $\mathbf{g}_{l}$ is a constant with regard to the variable $\mathbf{d}_{l}$:

Note that we only need to calculate $\mathbf{X}^{(k+1)}-\mathbf{D}^{(k)}\mathbf{Z}^{(k+1)}$ once to update $\{\mathbf{d}_{l}\}_{l=1}^{L}$ in parallel.

III-D Convergence Analysis

III-E Computational Complexity

The updating procedures of our two algorithms are quite straightforward, only requiring basic matrix and vector operations. To recover a signal that is sparse in the standard basis, C-PRIME has a time complexity $O(MN)$ under a general measurement matrix setting and $O(M\log M)$ under a DFT measurement matrix setting by exploiting fast Fourier transform and inverse fast Fourier transform. When the unknown signal is not sparse, SC-PRIME utilizes the sparse coding technique to approximate the unknown signal by a linear combination of a few columns in a dictionary. The time complexity of SC-PRIME is $O(LNP)$ to solve the undersampled phase retrieval task.

IV Numerical Results

In this section, we investigate the numerical performance of our MM-based algorithms, C-PRIME and SC-PRIME, and compare them with two up-to-date benchmark methods: UPRwO and DOLPHIn , respectively. First, we compare C-PRIME with UPRwO on the same randomly generated data. Later, we compare SC-PRIME with DOLPHIn on practical test images. Experimental results validate that C-PRIME and SC-PRIME outperform their corresponding benchmark method in terms of successful recovery probability and accuracy. All experiments are conducted on a personal computer with a $3.20$ GHz Intel Core i$5$-$4570$ CPU and $8.00$ GB RAM running Matlab R$2014$b.

We first investigate the performance of C-PRIME and compare it with the benchmark method UPRwO . To implement the UPRwO algorithm, we use the code provided on the authors’ homepagehttp://people.virginia.edu/~dsw8c/sw.html. In this subsection, we consider the clean measurements case and therefore we set the number of outliers to and signal-to-noise ratio (SNR) to infinity in the code. All other parameters are set as the default value.

The initialization steps of the UPRwO method are summarized below:

Generating the intensity measurements $\mathbf{y}=\left|\mathbf{Ax}_{o}\right|^{2}$.

For a fair comparison, in the simulation, we run the UPRwO code first with a fixed $(N,M,K)$ value. Besides the final results, we also store the original signal $\mathbf{x}_{o}$, the measurement matrix $\mathbf{A}$, and the intensity measurements $\mathbf{y}$. Later, we run our C-PRIME code using the same measurement matrix $\mathbf{A}$ and intensity measurements $\mathbf{y}$ from the UPRwO simulation.

In detail, the length of the original signal $N$ is set as the default value $128$. Since we consider the undersampled phase retrieval problem, the number of measurements is limited to be $M\in\{128,64,32,16,8\}$ and the sparsity level is set to be $K\in\{16,8,4,2\}$ (a value larger than $16$ ends up with unsuccessful recovery). For each of these possible $(M,K)$ pairs, experiments are conducted to test the performance of both algorithms provided with the same measurement matrix and intensity measurements. Note that under the DFT measurement matrix setting, any individual or combination of the following three trivial ambiguities conserve the Fourier magnitude:

Global constant phase shift: $\mathbf{x}\rightarrow\mathbf{x}\cdot e^{j\phi}$,

Circular shift: $[\mathbf{x}]_{i}\rightarrow[\mathbf{x}]_{(i+i_{0})\mod N}$,

Conjugate inversion: $[\mathbf{x}]_{i}\rightarrow\overline{[\mathbf{x}]_{N-i}}$.

Hence a disambiguation step is necessary to find the unique solution. For each solution $\mathbf{x}^{\star}$ returned by UPRwO and C-PRIME, we check all the possible candidates within the trivial ambiguities and choose the one with least normalized squared error (NSE) with regard to the original signal $\mathbf{x}_{o}$ as the final solution. The NSE between $\mathbf{x}^{\star}$ and $\mathbf{x}_{o}$ is calculated as

where the set $\mathcal{S}(\mathbf{x}^{\star})$ contains all the possible signals within the trivial ambiguities of $\mathbf{x}^{\star}$. Furthermore, since the original signal is generated as a random vector, experiments are repeated $100$ times for every $(M,K)$ pair using different and independent original signals with everything else fixed. The normalized mean squared error is calculated as the average of these $100$ independent NSE values. And among these $100$ independent trials, an algorithm is considered to successfully recover the original signal if the corresponding NSE is less than $10^{-4}$.

Final experimental results of UPRwO and C-PRIME are plotted in Fig. 1 on the successful recovery probabilities, and Fig. 2 on the normalized mean squared error. For sparse signals under different $(M,K)$ settings, our MM-based algorithm C-PRIME has a larger successful recovery probability and less normalized mean squared error than the benchmark algorithm UPRwO. The average CPU times of UPRwO and C-PRIME over these $100$ independent trials under all $(M,K)$ settings are presented in Table I. Both algorithms have a similar computational time.

IV-B SC-PRIME vs. DOLPHIn

We now investigate the performance of SC-PRIME and compare it with the benchmark method DOLPHIn on practical test images. To implement DOLPHIn, we use the code provided on the authors’ homepagehttp://www.mathematik.tu-darmstadt.de/~tillmann/#software. The test images are also downloaded from the same website. We choose the Gaussian measurement matrix setting and change the sampling rate from $4$ to $0.5$ ($M=0.5N$) to set up a valid undersampled phase retrieval problem. All other parameters are kept as the default value.

At the initialization step, the DOLPHIn algorithm takes the $2$D image as the original signal, thereby generating a random complex Gaussian measurement matrix, and generating the noisy intensity measurements with additive white Gaussian noise. The default SNR is $15$ dB. First, we run the DOLPHIn code and store the measurement matrix as well as the noisy intensity measurements. Later, these same noisy intensity measurements and the measurement matrix are provided as the input for SC-PRIME. Different from our problem setting in Section III where the patch signals are directly treated as the target signal, there is a sorting step between the $2$D image signal and the $2$D patch-based signal in . The corresponding change in the implementation of SC-PRIME is easy and trivial since changing the order of the elements in a matrix conserves its Frobenius norm.

To evaluate the quality of the reconstructed images, two standard image quality metrics are considered in this paper, namely the peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM). PSNR is the ratio between the maximum possible power of the original image and the mean squared error between the reconstructed image and the original image. It is usually expressed in terms of the logarithmic decibel scale, and the larger the value the better the image quality. SSIM reflects the structural similarities between the reconstructed image and the original image. It is on a scale from to $1$, and a larger value represents more similarities in the structure to the original image.

Final reconstruction results of DOLPHIn and SC-PRIME are presented in Fig. 3 on the $512\times 512$ color mandrill image, and Fig. 4 on the $2816\times 2112$ color waldspirale image. In both cases, our MM-based algorithm SC-PRIME can reconstruct the image with a larger PSNR and SSIM value as well as an impressively better visual quality than the benchmark algorithm DOLPHIn. Moreover, we summarize in Table II the PSNR and SSIM value of the reconstructed images for both algorithms on the rest of the test images. Besides the results of the reconstructed images ($\mathbf{X}^{\star}$), we also show the results of images approximated by the dictionary ($\mathbf{D}^{\star}\mathbf{Z}^{\star}$). All numbers in the table are averaged over $100$ Monte Carlo simulations using different and independent additive white Gaussian noise. The average CPU times over these $100$ independent trials are presented in Table III. It is interesting that the images approximated by the dictionary have a slightly larger PSNR and SSIM value than those reconstructed directly by the algorithms. Nevertheless, our MM-based algorithm SC-PRIME outperforms the benchmark method DOLPHIn in terms of PSNR and SSIM on all of the test images at the cost of slightly more CPU time.

V Conclusion

The undersampled phase retrieval problem draws great attention in various imaging applications. The difficulty lies in both theoretical analysis and practical algorithm design. Provided that only undersampled intensity measurements are available, one needs to solve a non-linear, non-convex, and under-determined inverse problem. In this paper, we have proposed two efficient algorithms based on the majorization-minimization framework that outperform the up-to-date benchmark methods in terms of successful recovery probability and accuracy under different problem settings. The first algorithm C-PRIME can uniquely recover a signal that is sparse in the standard basis. When the signal is not sparse itself, the second algorithm SC-PRIME utilizes the sparse coding technique to approximate the target signal by a linear combination of a few columns in a dictionary. Experimental results on randomly generated data and practical test images are also provided in the paper to further validate the efficiency of our algorithms with impressive results.

Appendix A Justification for Using Modulus Information

Recall the $M$ noisy intensity measurements

and we assume $y_{i}\geq 0$ (otherwise we just discard this measurement). And the noise $n_{i}$ is assumed to be independent of the measurements. Therefore,

Usually the noise level is much smaller than the value of the clean intensity measurements, $|n_{i}|\ll\left|\mathbf{a}_{i}^{H}\mathbf{x}\right|^{2}$. It is sufficient to make the following approximation taking the first two terms in the Taylor series:

The first term $\left|\mathbf{a}_{i}^{H}\mathbf{x}\right|$ is the actual clean modulus measurement. And the second term can be regarded as the additive noise, with expectation and variance

Therefore, the additive noise to the modulus information $\left|\mathbf{a}_{i}^{H}\mathbf{x}\right|$ has a lesser expectation and variance value than the additive noise to the intensity information $\left|\mathbf{a}_{i}^{H}\mathbf{x}\right|^{2}$ when $\left|\mathbf{a}_{i}^{H}\mathbf{x}\right|>\frac{1}{2}$.

Note that the dictionary $\mathbf{D}^{(k)}=[\mathbf{d}_{1}^{(k)},\ldots,\mathbf{d}_{L}^{(k)}]$ satisfies $\|\mathbf{d}_{l}^{(k)}\|_{2}\leq 1,\forall l=1,\ldots,L$, so

The equality is achieved when all $\{\mathbf{d}_{l}^{(k)}\}$ lie on the same line and $\|\mathbf{d}_{l}^{(k)}\|_{2}=1,\forall l$.

References