On phase retrieval via matrix completion and the estimation of low rank PSD matrices

TL;DR

This paper introduces a new non-convex optimization algorithm for estimating low-rank PSD matrices from underdetermined measurements, with applications to phase retrieval and analysis of oversampling effects on stability.

Contribution

It proposes a novel algorithm for low-rank PSD matrix recovery and provides theoretical insights into the stability of phase retrieval under oversampling.

Findings

The algorithm effectively estimates low-rank PSD matrices from limited data.

Oversampling improves the stability of the phase retrieval problem.

Theoretical analysis links oversampling to enhanced robustness in matrix recovery.

Abstract

Given underdetermined measurements of a Positive Semi-Definite (PSD) matrix $X$ of known low rank $K$, we present a new algorithm to estimate $X$ based on recent advances in non-convex optimization schemes. We apply this in particular to the phase retrieval problem for Fourier data, which can be formulated as a rank 1 PSD matrix recovery problem. Moreover, we provide theory for how oversampling affects the stability of the lifted inverse problem.

Tables

Table 1. Table 1 : Reweighted nuclear norm eigenvalues table

	Exp 1	Exp 2	Exp 3	Exp 4	Exp 5
${\lambda }_1$	0.9402	0.9419	0.9360	0.9461	0.9635
${\lambda }_2$	0.0489	0.0579	-0.0487	0.0341	0.0439
${\lambda }_3$	-0.0263	-0.0327	0.0339	0.0318	0.0369
${\lambda }_4$	0.0260	-0.0212	0.0227	-0.0298	0.0274
${\lambda }_5$	0.0173	0.0158	-0.0172	-0.0177	0.0155
${\lambda }_6$	-0.0119	0.0104	0.0100	0.0090	-0.0118
${\lambda }_7$	0.0051	0.0046	0.0073	0.0061	0.0079
${\lambda }_8$	0	0.0003	0.0022	0.0050	-0.0021
${\lambda }_9$	0	0	0	0.0031	0
${\lambda }_{10}$	0	0	0	0	0

Equations

$$b_k=|⟨{\mathbf{a}}_k,\mathbf{x}⟩{|}^2,k=1,\dots ,M$$

$$\underset{\mathbf{x}\in {\mathbb{C}}^N}{\min }\sum _{k=1}^M(|⟨{\mathbf{a}}_k,\mathbf{x}⟩{|}^2-b_k{)}^2,$$

$$\left(\frac{1}{c^2}\frac{\partial }{\partial t^2}-{\nabla }^2\right)\Psi (x,y,z,t)=0,$$

$$\Psi (x,y,z,t)=\frac{1}{\sqrt{2\pi }}{\int }_0^{\infty }{\psi }_{\omega }(x,y,z)\exp (-i\omega t)d\omega$$

$$({\nabla }^2+k^2){\psi }_{\omega }(x,y,z)=0,k=\omega /c$$

$${\psi }_{\omega }(x,y,z=\Delta )⟶-\frac{ik\exp ik\Delta }{\Delta }\exp \left[\frac{ik}{2\Delta }(x^2+y^2)\right]\widehat{{\psi }_{\omega }}\left(k_x=\frac{kx}{\Delta },k_y=\frac{ky}{\Delta },z=0\right)$$

$$\{\psi \in C_c^{\infty }:|\hat{\psi }|=|\widehat{{\psi }_O}|\}\text{ and }\{\psi \in C_c^{\infty }:\text{supp}(\psi )\subseteq C\}.$$

$$X\mapsto \left(\begin{array}{c}⟨A_1,X{⟩}_F\\ ⟨A_2,X{⟩}_F\\ ⋮\\ ⟨A_M,X{⟩}_F\end{array}\right)$$

$$\mathcal{A}(X)=\mathbf{b}.$$

$$\min \Vert \mathcal{A}(X)-\mathbf{b}{\Vert }^2\text{subject to }\text{rank}(X)=1,X\succcurlyeq 0.$$

$$\mathrm{argmin}\lambda \Vert X{\Vert }_{\ast }+\frac{1}{2}\Vert \mathcal{A}(X)-\mathbf{b}{\Vert }_2^2\text{subject to }X\succcurlyeq 0.$$

$$\mathcal{R}(\mathbf{x})=f(\mathbf{x})+\frac{1}{2}\Vert \mathcal{A}(\mathbf{x})-\mathbf{b}{\Vert }^2.$$

$${\mathcal{Q}}_{\gamma }(f)(\mathbf{x}):=\underset{\alpha \in \mathbb{R},\mathbf{y}\in \mathcal{H}}{\sup }\left\{\alpha -\frac{\gamma }{2}\Vert \mathbf{x}-\mathbf{y}{\Vert }^2:\alpha -\frac{\gamma }{2}\Vert \cdot -\mathbf{y}{\Vert }^2\le f\right\}$$

$${\mathcal{R}}_{reg}(\mathbf{x})={\mathcal{Q}}_{\gamma }(f)(\mathbf{x})+\frac{1}{2}\Vert \mathcal{A}(\mathbf{x})-\mathbf{b}{\Vert }^2$$

$$f(\mathbf{x})={\iota }_K(\mathbf{x})=\{\begin{array}{ll}0& \text{if }\Vert \mathbf{x}{\Vert }_0\le K\\ +\infty & \text{otherwise}.\end{array}$$

$$f(X)={\iota }_{R_K}(X)=\{\begin{array}{ll}0& \text{if }\text{ rank}(X)\le K\\ +\infty & \text{otherwise}.\end{array}$$

$${\mathcal{Q}}_{\gamma }({\iota }_{R_K})(X)={\mathcal{Q}}_{\gamma }({\iota }_K)(\sigma (X))$$

$${\iota }_K^{+}(\mathbf{x})=\{\begin{array}{ll}0& \text{if }\Vert \mathbf{x}{\Vert }_0\le K\text{ and }\mathbf{x}\ge 0\\ +\infty & \text{otherwise}.\end{array},$$

$$\underset{X\in {\mathbb{H}}_N^{+}:\text{ rank}(X)\le K}{\mathrm{argmin}}\Vert \mathcal{A}(X)-\mathbf{b}{\Vert }^2$$

$$\underset{X\in {\mathbb{H}}_N}{\mathrm{argmin}}{\iota }_{R_K^{+}}(X)+\frac{1}{2}\Vert \mathcal{A}(X)-\mathbf{b}{\Vert }^2.$$

$$\underset{X\in {\mathbb{H}}_N}{\mathrm{argmin}}{\mathcal{Q}}_{\gamma }({\iota }_{R_K^{+}})(X)+\frac{1}{2}\Vert \mathcal{A}(X)-\mathbf{b}{\Vert }^2$$

$${\text{prox}}_g(w)=\underset{v\in \mathcal{V}}{\mathrm{argmin}}\left(g(v)+\frac{1}{2}\Vert v-w{\Vert }_{\mathcal{V}}^2\right).$$

$${\text{prox}}_{{\mathcal{Q}}_{\gamma }(F)/\rho }(X)=U\text{diag}({\text{prox}}_{{\mathcal{Q}}_{\gamma }(f)/\rho }(\lambda (X)))U^{\ast },$$

$${\mathcal{S}}_{\gamma }(F)(Y)=\underset{X}{\sup }-F(X)-\frac{\gamma }{2}\Vert X-Y{\Vert }_F^2=\underset{X}{\sup }-f(\lambda (X))-\frac{\gamma }{2}\left(\Vert X{\Vert }_F^2+\Vert Y{\Vert }_F^2-2\mathsf{Re}⟨X,Y{⟩}_F\right)$$

$${\mathcal{S}}_{\gamma }(F)(Y)=\underset{\xi \in {\mathbb{R}}^N}{\sup }-f(\xi )-\frac{\gamma }{2}\Vert \xi -\lambda (Y){\Vert }^2={\mathcal{S}}_{\gamma }(f)(\lambda (Y)).$$

$$\lambda (Z)-\lambda (X)=\lambda (Z-X)\in \partial {\mathcal{Q}}_{\gamma }(f)(\lambda (Z))/\rho .$$

$$|⟨x,a_{\mathbf{k}}⟩{|}^2=b_{\mathbf{k}},\mathbf{k}\in S.$$

$$⟨X,a_{\mathbf{k}}\otimes {\overline{a}}_{\mathbf{k}}⟩=b_{\mathbf{k}},\mathbf{k}\in S.$$

$$f_{\mathbf{k}}(\mathbf{n})=e^{-2\pi i\frac{\mathbf{k}\cdot \mathbf{n}}{m}}$$

$$f_{\mathbf{k}}(\mathbf{n})=e^{i{\zeta }_{\mathbf{k}}\cdot \mathbf{n}}$$

Full text

Keywords: Fourier phase retrieval, low rank matrices, non-convex optimization.

MSC2010: 49M20, 49N45, 65K10, 90C26.

1 Introduction

The (discretized) formulation of the phase retrieval problem consists in finding a complex vector $\mathbf{x}$, usually a discretized signal or image, given a certain amount of measurements in form of modulus of scalar products (see [SEC15]), i.e.

[TABLE]

(where $\mathbf{a}_{k}\in\mathbb{C}^{N}$ are known), and possibly additional geometric constraints. The aim is thus to reconstruct the discrete vector $\mathbf{x}\in\mathbb{C}^{N}$ representing the object. In [BCE06] is shown that $M\geq 4N-2$ generic complex measurements are needed in order to be able to distinguish two different signals (up to a phase); in [CEH15] the lower bound was improved to $M\geq 4N-4$. This number is believed to be close to optimal.

A possible combinatorial approach to the problem has been described in [SCa91], together with the proof that, as formulated, the problem is NP-hard. An optimization-based standard solution technique would be to solve the least-squares minimization problem

[TABLE]

but the quadratic terms in the objective function makes the problem highly non-convex and therefore hard to solve. A convex optimization approach was suggested in [CSV13], called PhaseLift, and this note considers an improvement of this approach. We recommend [FWd15] for a pleasant overview of recent advances regarding PhaseLift and related techniques.

A lot of work has been devoted to deal with cases in which $\mathbf{a}_{k}$ are Gaussian measurements, see e.g. [CLS15], but these algorithms seems to work poorly with Fourier data. The main objective of the present article is to provide improvements of the PhaseLift approach, with a particular focus on Fourier measurements, based on recent advances in non-convex low rank approximation algorithms. The theory relies on the quadratic envelope applied to the indicator functional of matrices with a fixed predetermined rank, and can be used also for any other PSD-estimation problem of fixed low rank (Section 2.1).

1.1 Physical background

The phase retrieval problem appears frequently for instance in elecrodynamics; the single scalar complex field $\Psi(x,y,z,t)$, often called wave-function, solution of the d’Alembert equation

[TABLE]

is enough to describe the electromagnetic disturbance in the free space. Making use of Fourier integral, the wave-function can be spectrally decomposed as superposition of monochromatic fields

[TABLE]

where the spatial wave-function $\psi_{\omega}(x,y,z)$ associated with a given monochromatic component of the spectral decomposition of $\Psi$ solves the Helmholz equation

[TABLE]

with suitable boundary conditions. In Coherent Diffractive Imaging (CDI) an object (sample) is illuminated by a coherent almost monochromatic wavefield and subsequentially the diffracted far-field intensity pattern is measured by detectors. The current detectors measure intensities but they are not able to measure the phase of the diffracted pattern; the phase retrieval problem consists then in retrieving the object knowing the aforementioned intensity plus, in general, some additional physical constraints. In the Fraunhofer regime, which occurs whenever the Fresnel number $N_{F}:=b^{2}/\lambda\Delta\ll 1$ ($b$ being the diameter of the region occupied by the sample and $\Delta$ the distance between the sample and the detector) it can be shown (see for instance [Pag06]) that the $2$-dimensional propagated disturbance $\psi_{\omega}$ measured at $z=\Delta\geq 0$ (where the optical axis is considered to be coincident with the $z$ axis) is

[TABLE]

where $k$ is the wavenumber and $\widehat{\cdot}$ the $2$-dimensional Fourier transform. In words this says that "the propagated disturbance has the form of a paraxial modulated spherical wave with the modulation being proportional to (a transversely scaled form of) the two-dimensional Fourier transform of the unpropagated disturbance" ([Pag06]).

A family of algorithms developed for solving the phase retrieval problem in combination with a priori knowledge of the support of the object, relies on the simple idea of alternatingly adjust the support of the image and the modulus of its Fourier transform. This goes back to Gerchberg and Saxton ([GSa72]) and to Fienup ([Fie82]), even though the pool of mathematical ideas, borrowed from convex analysis, goes back at least to Bregman ([Bre65]). The most famous is the error-reduction algorithm which is essentially a non-convex adaptation of the idea of Projection On Convex Sets (cfr. POCS algorithm, chatper 5 of [BCo10]): if $\psi_{O}$ is the object we want to reconstruct, supported in $C\subseteq\mathbb{R}^{2}$, error-reduction performs successive projections between the two sets

[TABLE]

The drawback is that the first set is not convex and therefore there are no theoretical guarantees of the convergence of this scheme. Further algorithms have of course been developed; in [Mar07] the projections are replaced with reflections, in [MBK16] a Newton-type scheme with Tikhonov regularization is proposed, to mention two recent contributions. These methods have been successfully applied with real data, but the non-convexity issue persists.

1.2 PhaseLift approach

In the present note we focus our attention to a lifting approach popularized by [CSV13] and [CES13]. For $A_{k}=\mathbf{a}_{k}\mathbf{a}_{k}^{*}$, where ∗ is the conjugate transpose and the latter product is the usual matrix multiplication, we define a linear operator $\mathcal{A}:{\mathbb{H}}_{N}(\mathbb{C})\to\mathbb{C}^{M}$ via

[TABLE]

being $\langle\cdot,\cdot\rangle_{F}$ the Frobenius inner product and ${\mathbb{H}}_{N}$ the space of $N\times N$ Hermitian matrices. Noticing that $|\langle\mathbf{a}_{k},\mathbf{x}\rangle|^{2}=\langle\mathbf{a}_{k}\mathbf{a}_{k}^{*},\mathbf{x}\mathbf{x}^{*}\rangle_{F}$ and recalling that every positive semidefinite matrix $X\in{\mathbb{H}}_{N}(\mathbb{C})$ with $\text{rank}(X)=1$ admits a factorization of the type $\mathbf{x}\mathbf{x}^{*}=X$ with $\mathbf{x}\in\mathbb{C}^{N}$, the quadratic phase retrieval problem is lifted onto a linear one with dimension squared; the problem (1) to be solved is now to find a suitable rank $1$ PSD-matrix $X$ such that

[TABLE]

In real applications, the “perfect measurement” $b_{k}=|\langle\mathbf{a}_{k},\mathbf{x}\rangle|^{2}$ is contaminated by noise, so the problem can be re-casted as an optimization one

[TABLE]

The constraint $\text{rank}(X)=1$ is however highly non-convex, so a convex relaxation has been proposed in [CSV13] based on the nuclear norm, known for being “low-rank inducing” (cfr. [RFP10]):

[TABLE]

The parameter $\lambda$ gives a tradeoff between low rank of the sought minimum $X$ and data fit, and this term also induces a bias which scales with $\lambda$. For this reason, in practice one does not in general find a rank 1 matrix, but a low rank one from which a suitable vector candidate $\mathbf{x}$ can be extracted, we refer to [CSV13] for the details.

1.3 Innovations and discussion of the lifting trick

While recognizing the PhaseLift as an important and celebrated contribution, we have found that this method suffers from a number of drawbacks making its practical use limited. The main issue is the size of the involved objects due to the lifting process. A typical (CDI) measurement generates an image of at least $100\times 100$ which leads to a vector $\mathbf{x}$ of size $10^{4}$, which in turn yields a $10^{4}\times 10^{4}$-matrix which is then fed to an iterative algorithm which need to compute eigenvalue factorization. The time complexity of the latter, without any implementation tricks, is $O((10^{4})^{3})=O(10^{12})$.

The second drawback stems from the fact that one usually need to run the algorithm several times in order to find a suitable value of $\lambda$, and a third drawback is that one typically need lots of linearly independent measurements $(A_{1},\ldots,A_{M})$ for stable recovery. As mentioned initially, one needs at least $M=4N$ measurements for having an (essentially) unique solution in the noiseless case. [CSV13] suggests working with $M\geq cN\log(N)$ where $c$ is some (unknown) constant, whereas in their numerical section they use $M=6N$ where $N=128$. When dealing with Fourier data, one commonly used method for boosting $M$ is simply oversampling in the Fourier domain [MSC98], but it was observed numerically in [CES13] that this gives an ill-posed inverse problem (see Sections 2.1 and 4.4.3), which led them to instead suggest the use of e.g. masks to increase the amount of measurements.

In this article we will prove that oversampling indeed does not give rise to more linearly independent data beyond $M\approx 2^{d}N$, where $d$ is the dimension of the problem. In other words it is pointless to oversample more than a factor 2 in each dimension. This also shows that any method based on the lifting principle (3) is bound to fail unless combined with masks or similar tricks, at least for $d\leq 2$. Furthermore, we will rewrite (4) into an optimization problem where the constrains are built into the functional to be minimized, and then use the quadratic envelope to regularize this functional. Ideally one would like to work with the convex envelope of the functional, but this is not doable in practice. The quadratic envelope yields a non-convex continuous surrogate which has recently been shown to have the desirable property of yielding a regularizer which, just like the convex envelope, does not move the global minima [Car19]. While not completely removing other local minimizers, it reduces the amount and in practice it seems that the corresponding algorithm can find the global minimizer in fairly difficult problems. In other words we suggest a new framework which seems to be able to solve the original problem (4). We refrain from proving this claim mathematically, but remark that in [CGO19] an analogous statement is rigorously shown for the similar problem of low rank recovery without PSD constraints. Here we satisfy with developing the algorithm, with a particular focus on efficient implementation, as well as numerically demonstrating several desirable features (see Section 5);

• •

the algorithm can solve any problem of the type (4) but with a predetermined rank $K$, not necessarily $K=1$.

• •

The nuclear norm approach (5) yields matrices with rank larger than 1, hence using more degrees of freedom. Despite this, the proposed method gives a better data fit while at the same time fulfilling the rank constraint.

• •

In [CES13] a reweighted version of (5) is used, which is known for better accuracy, (but this algorithm is non-convex as well). In terms of reconstruction accuracy the two algorithms now perform similarly. Benefits of the proposed algorithm is that it finds the correct rank and does so within fewer iterations, plus that it relies on a more developed theoretical framework.

The paper is structured as follows; Section 2 presents the mathematical details behind the new general fixed rank PSD-matrix estimation algorithm, Section 3 presents our theoretical findings regarding the application to phase retrieval and sampling issues, Section 4 covers implementation details and finally Section 5 concludes with some numerical illustrations.

1.4 Notation

${\mathbb{M}}_{N}$ will denote the set of $N\times N$ complex matrices and ${\mathbb{H}}_{N}$ the subset of Hermitian (i.e. self-adjoint) matrices. $N$ is the size of the unknown vector before lifting to a problem with $N^{2}$ unknowns, and $M$ is the amount of measurements, which we assume is larger than $N$. Images and higher dimensional objects are represented as tensors $\otimes_{j=1}^{d}{\mathbb{C}}^{n}$ and the measurements then take place in the Fourier domain which we discretize as $\otimes_{j=1}^{d}{\mathbb{C}}^{m}$, where $m\geq n$. Hence $N=n^{d}$ and $M$ is a constant multiple $N_{m}+1$ of $m^{d}$, where $N_{m}$ is the number of masks. ${\mathcal{A}}$ denotes the operator (3) and its “matrixification” (as described in Section 4.1) is denoted by A, except in the case of pure Fourier data, in which case we denote these objects by ${\mathcal{F}}$ and F respectively. The quadratic envelope is denoted by ${\mathcal{Q}}_{\gamma}(f)$. $\circ$ denotes composition of two functions, in particular if $f$ acts on ${\mathbb{R}}^{N}$ then $f\circ\lambda$ acts on ${\mathbb{H}}_{N}$, where $\lambda$ denotes the eigenvalues of any given $X\in{\mathbb{H}}_{N}$ ordered non-increasingly.

2 Quadratic envelope approach to low rank recovery

Let $\mathcal{H}$ be a Hilbert space, $f:\mathcal{H}\to\mathbb{R}\cup\{\infty\}$ a non-convex penalty functional and set

[TABLE]

where ${\mathcal{A}}$ is a linear operator from ${\mathcal{H}}$ to some other Hilbert space in which the “measurement” ${\bf b}$ lives. In [Car19] the quadratic envelope

[TABLE]

is studied as a means of regularizing functionals of this type. It is shown that the regularized functional

[TABLE]

has some desirable properties when $\|\mathcal{A}\|^{2}<\gamma$. Most notably, $\mathcal{R}_{reg}$ lies between $\mathcal{R}$ and its lower semi-continuous convex envelope, any local minimizer of $\mathcal{R}_{reg}$ is a local minimizer of $\mathcal{R}$ and the global minimizers of $\mathcal{R}$ and $\mathcal{R}_{reg}$ coincide.

In [CGO18] the quadratic envelope has been studied further for the case of the famous $\ell^{0}-\ell^{2}$ problem in which $\mathcal{H}=\mathbb{R}^{N}$ and $f({\bf x})=\|{\bf x}\|_{0}$ or, in case the sought degree of sparsity $K$ is known,

[TABLE]

Here $\|x\|_{0}$ denotes the number of non-zero entries in ${\bf x}$. In these cases $\mathcal{R}_{reg}$ has been numerically compared to the more classical $\ell^{1}$-penalty in solving the problem of retrieving a sparse signal. In the companion work [CGO19] the analysis has been lifted to the space of $N\times N$-matrices $\mathcal{H}=\mathbb{M}_{N}(\mathbb{C})$, and $f(X)={\text{ rank}}(X)$ or

[TABLE]

The vector problem and the matrix problem are closely related; note that $\iota_{R_{K}}(X)=\iota_{K}(\sigma(X))$ where $\sigma(X)$ denotes the singular values of $X$. It turns out that

[TABLE]

and, more importantly, if $X$ has SVD $X=U\Sigma V^{*}$, that ${\text{prox}}_{\iota_{R_{K}}}(X)=U{\text{diag}}({\text{prox}}_{\iota_{K}}(\sigma(X)))V^{*}$. In other words, if we can compute the proximal operator in the scalar case, the matrix case follows immediately. In [CGO19] the theory is much further developed for the case $f=\iota_{R_{K}}$; in particular we show that for noisy measurements the functional $\mathcal{R}_{reg}$ has a unique local minima which is also the solution of $\operatorname*{arg\,min}_{{\text{ rank}}(X)\leq K}\|{\mathcal{A}}(X)-{\bf b}\|$. The new ingredient in the present contribution is basically the PSD condition, and we satisfy with numerically testing the machinery on the Phase Retrieval problem as well as some generic low rank recovery problems with $K>1$.

2.1 Low-rank Positive Semi-Definite problems

We now come to the new material for this paper. In this section we set ${\mathcal{H}}$ to be the set of Hermitian matrices ${\mathbb{H}}_{N}\subset\mathbb{M}_{N}(\mathbb{C})$ and consider the problem of searching for a PSD matrix with fixed rank $K$. Of course, for the PhaseLift problem we will set $K=1$, but we develop the theory for general $K\in{\mathbb{N}}$. Motivated by the previous section we introduce the function $\iota_{K}^{+}:\mathbb{R}^{N}\to\{0,+\infty\}$ defined by

[TABLE]

where ${\bf x}\geq 0$ is interpreted elementwise. If $R_{K}^{+}$ denotes the set of positive semi-definite matrices with rank $\leq K$, then it is easy to see that $\iota_{R_{K}^{+}}(X)=\iota_{K}^{+}(\lambda(X))$ where $\lambda(X)$ denotes the eigenvalues of $X\in{\mathcal{H}}$. Given a transform ${\mathcal{A}}$ we then have that the fixed-rank minimization problem

[TABLE]

(where ${\mathbb{H}}_{N}^{+}$ denotes the subset of PSD matrices) is equivalent with

[TABLE]

This in turn has the same global minimizer as the regularized problem

[TABLE]

as long as $\gamma>\|{\mathcal{A}}\|^{2}$, as argued earlier. The latter is a continuous (in its domain ${\mathbb{H}}_{N}^{+}$, cfr. Proposition 3.2 of [Car19]) functional whose stationary points can be found for instance via Forward-Backward Splitting scheme (also know as Proximal Gradient Method: see for instance [Bec17])111This statement has a theoretical foundation: for a lower semicontinuous and semialgebric $f:\mathbb{R}^{N}\to\mathbb{R}\cup\{\infty\}$, the function ${\mathcal{Q}}_{\gamma}(f):\mathbb{R}^{N}\to\mathbb{R}\cup\{\infty\}$ is lower semicontinuous and semialgebric ([Car19]), and therefore by a theorem due to Bolte-Daniilidis-Lewis ([BDL07]) it has the local Lojasiewicz-Kurdyka property. In [ABS13] it is showed that the sequence generated by FBS, a particular case of their Algorithm 3, converges to a stationary point of the functional, or else is unbounded. , as long as the proximal operator of ${\mathcal{Q}}_{\gamma}(\iota_{R_{K}^{+}})$ is computable. We now describe how this can be done. We recall that a function $f:\mathbb{R}^{N}\to\mathbb{R}$ is said to be symmetric if $f({\bf x})=f(\Pi{\bf x})$ for every permutation $\Pi$ and for every ${\bf x}\in\mathbb{R}^{N}$. We also recall that for a lower semi-continuous function $g:\mathcal{V}\to[-\infty,\infty]$, being $\mathcal{V}$ an Hilbert space, the (possibly empty or set-valued) proximal operator is defined by

[TABLE]

Proposition 2.1**.**

Let $f:{\mathbb{R}}^{N}\rightarrow[0,\infty]$ be a symmetric function and consider $F:{\mathbb{H}}_{N}\rightarrow[0,\infty]$ defined by $F(X)=f(\lambda(X))$. Then ${\mathcal{Q}}_{\gamma}(F)={\mathcal{Q}}_{\gamma}(f)\circ\lambda$ and, for $\rho>\gamma$, we have

[TABLE]

where $U\emph{diag}(\lambda(X))U^{*}$ is a spectral decomposition of $X$.

Proof.

By Proposition 3.1 in [Car19], we have ${\mathcal{Q}}_{\gamma}(f)={\mathcal{S}}_{\gamma}({\mathcal{S}}_{\gamma}(f))$, where ${\mathcal{S}}_{\gamma}(f)(y)=\sup_{x}-f(x)-\frac{\gamma}{2}\|x-y\|^{2}$. Note that

[TABLE]

and that ${\mathsf{Re}}\langle X,Y\rangle_{F}\leq\sum_{j=1}^{N}{\lambda}_{j}(X){\lambda}_{j}(Y)$, by a variation of the von Neumann trace inequality (cfr. [Lew96]). Moreover the equality holds if and only if there exists a unitary matrix $V$ such that $V^{*}XV=\text{diag}{\lambda}(X)$ and $V^{*}YV=\text{diag}{\lambda}(Y)$. We may therefore assume without loss of generality that $X=V\text{diag}(\xi)V^{*}$ for some $\xi\in\mathbb{R}^{N}$, being $V\text{diag}\lambda(Y)V^{*}$ a spectral decomposition of $Y$. For such $X$ we have that $\langle X,Y\rangle_{F}=\sum\xi_{j}\lambda_{j}(Y)$, (where we do not require that $\xi$ is non-decreasing). By the symmetry of $f$ we have that $f(\lambda(X))=f(\xi)$, and therefore

[TABLE]

The corresponding claim for ${\mathcal{Q}}_{\gamma}$ follows immediately by ${\mathcal{Q}}_{\gamma}(f)={\mathcal{S}}_{\gamma}({\mathcal{S}}_{\gamma}(f))$.

Let us prove the second part of the statement. We have that $Z\in\text{prox}_{{\mathcal{Q}}_{\gamma}(F)/\rho}(X)$ if and only if $Z$ minimizes the functional $Y\mapsto{\mathcal{Q}}_{\gamma}(F)(Y)+\frac{\rho}{2}\|X-Y\|_{F}^{2}$ which is convex, since ${\mathcal{Q}}_{\gamma}(F)(Y)+\frac{\gamma}{2}\|X-Y\|_{F}^{2}$ is the l.s.c. convex envelope of $F(Y)+\frac{\gamma}{2}\|X-Y\|_{F}^{2}$ (by Theorem 3.1 [Car19]) and ${\mathcal{Q}}_{\gamma}(F)(Y)+\frac{\rho}{2}\|X-Y\|_{F}^{2}={\mathcal{Q}}_{\gamma}(F)(Y)+\frac{\gamma}{2}\|X-Y\|_{F}^{2}+\frac{\rho-\gamma}{2}\|X-Y\|_{F}^{2}$. This in turn happens if and only if $Z-X\in\partial{\mathcal{Q}}_{\gamma}(F)(Z)/\rho=\partial({\mathcal{Q}}_{\gamma}(f)\circ\lambda)(Z)/\rho$. From Section 5.2 of [BoL05] follows that the latter holds if and only if $Z-X$ and $Z$ share the same eigenvectors and

[TABLE]

This is equivalent to $\lambda(Z)\in\text{prox}_{{\mathcal{Q}}_{\gamma}(f)/\rho}(\lambda(X))$ and by the symmetry of $f$ we have that $\lambda_{i}(Z)=\lambda_{j}(Z)$ whenever $\lambda_{i}(X)=\lambda_{j}(X)$. Hence any unitary matrix $U$ of eigenvectors to $X$ are also eigenvectors to $Z$, and since $Z=U\text{diag}(\lambda(Z))U^{*}$, the desired conclusion follows. ∎

We remark that the maximum negative quadrature of ${\mathcal{Q}}_{\gamma}(F)$ is $\gamma$, as shown in [Car19], so the condition $\gamma<\rho$ ensures that the proximal operator is single valued (since the corresponding minimization problem is strongly convex).

Of course we are interested in the concrete case of $f=\iota_{K}^{+}$, but it turns out that ${\mathcal{Q}}_{\gamma}(\iota_{K}^{+})$ has a rather complicated expression. Luckily, the related transform ${\mathcal{S}}_{\gamma}(\iota_{K}^{+})$ has a simpler expression, and it follows that we can still compute the ${\text{prox}}_{{\mathcal{Q}}_{\gamma}(\iota_{K}^{+})}$ by duality; we postpone the details of this to Section 4.2. In the coming section we investigate the structure of the operator ${\mathcal{A}}$ for the phase retrieval problem with multidimensional data and Fourier measurements.

3 Images and Fourier data

Let us return to the phase retrieval problem. We consider “images” in $d$ dimensions, which can be realized as the tensor product vector space $\otimes_{j=1}^{d}{\mathbb{C}}^{n}$, where we use the same number of points in every dimension for simplicity only. This linear vector space can of course be identified with ${\mathbb{C}}^{N}={\mathbb{C}}^{n^{d}}$ by introducing some basis, but we will see that it is often convenient to actually skip this step and work directly in the more abstract setting. In this case a rank 1 matrix corresponds to a linear operator on $\otimes_{j=1}^{d}{\mathbb{C}}^{n}$ of the form $x\otimes y$, and the matrix is PSD if and only if the operator can be written $x\otimes\overline{x}$, where the bar indicates complex conjugation.

Given elements $a_{\mathbf{k}}\in\otimes_{j=1}^{d}{\mathbb{C}}^{n}$, where $\mathbf{k}$ runs over some set of (usually multidimensional) subindices $S$ (where $M=\#S$), we seek one element $x\in\otimes_{j=1}^{d}{\mathbb{C}}^{n}$ such that

[TABLE]

In this new framework, PhaseLift corresponds to the equivalent problem of finding a rank 1 PSD linear operator $X$ on $\otimes_{j=1}^{d}{\mathbb{C}}^{n}$ such that

[TABLE]

3.1 Fourier data and limitations to oversampling

In the common case of Fourier measurements, the $a_{\mathbf{k}}$’s are discretizations of pure oscillatory exponential functions, and in this case we denote them by $f_{\mathbf{k}}$ and the corresponding operator by ${\mathcal{F}}$ in place of ${\mathcal{A}}$. We denote by $\ell^{2}(S)$ the linear vector space of all functions on $S$, so that $b$ naturally identifies with an element of $\ell^{2}(S)$. Often, we measure on an $m^{d}$ grid where $(m/n)^{d}$ is the oversampling factor, i.e. $S=\{0,\ldots,m-1\}^{d}$ in which case $\ell^{2}(S)$ is readily identified with $\otimes_{j=1}^{d}{\mathbb{C}}^{m}$. To be more explicit, in this case we have

[TABLE]

for $\mathbf{k}\in\{0,\ldots,m-1\}^{d}$ and ${\bf n}\in\{0,\ldots,n-1\}^{d}$. However, to allow for unequally spaced sampling we stick to the more general setting where

[TABLE]

and $\{\zeta_{\mathbf{k}}\in{\mathbb{R}}^{d}\}_{\mathbf{k}\in S}$ are some frequencies and $S$ some set. The equation $\eqref{m1}$ can now be written

[TABLE]

where ${\mathcal{F}}:(\otimes_{j=1}^{d}{\mathbb{C}}^{n})\otimes(\otimes_{j=1}^{d}{\mathbb{C}}^{n})\longrightarrow\ell^{2}(S)$ is the linear operator defined by ${\mathcal{F}}(X)(\mathbf{k})=\left\langle X,f_{\mathbf{k}}\otimes\overline{f_{\mathbf{k}}}\right\rangle$.

Since we are working with an ill-posed inverse problem, it is crucial to get as much linearly independent equations as possible. In other words we want to choose $\{\zeta_{\mathbf{k}}\in{\mathbb{R}}^{n}\}_{\mathbf{k}\in S}$ so that ${\mathcal{F}}$ has maximal rank. The next result basically states that it is pointless to oversample beyond a factor of two.

Lemma 3.1**.**

Independent of how $\{\zeta_{\mathbf{k}}\in{\mathbb{R}}^{d}\}_{\mathbf{k}\in S}$ is chosen, the maximal rank of ${\mathcal{F}}$ is $(2n-1)^{d}$.

Proof.

We let $\{e_{\bf n}\}_{{\bf n}\in\{1,\ldots,n\}^{2d}}$ denote the canonical basis in $\otimes_{j=1}^{2d}{\mathbb{C}}^{n}$. The range of ${\mathcal{F}}$ is then spanned by the functions ${\mathcal{F}}(e_{\bf n})$ which, if we write ${\bf n}=({\bf n}^{1},{\bf n}^{2})$ with ${\bf n}^{1},{\bf n}^{2}\in\{1,\ldots,n\}^{d}$, takes the form

[TABLE]

The amount of tuples of the form ${\bf n}^{1}-{\bf n}^{2}$ equals $(2n-1)^{d}$, which gives the desired upper bound. ∎

We now prove that, given sufficient oversampling and a vise choice of $\{\zeta_{\mathbf{k}}\}_{\mathbf{k}\in S}$, the rank of ${\mathcal{F}}$ actually equals the maximal rank $(2n-1)^{d}$. More precisely, we will need that $\{\zeta_{\mathbf{k}}\}_{\mathbf{k}\in S}$ is such that ${\bf n}\mapsto e^{i{\zeta_{\mathbf{k}}\cdot{\bf n}}}$ is a linearly independent set on $\otimes_{j=1}^{d}{\mathbb{C}}^{2n-1}$ or, if $|S|>(2n-1)^{d}$, that these functions span the full space $\otimes_{j=1}^{d}{\mathbb{C}}^{2n-1}$. Such a choice of $\{\zeta_{\mathbf{k}}\}_{\mathbf{k}\in S}$ will be called non-degenerate. By considering the DFT it is clear that non-degenerate sets of frequencies exist for all cardinalities of $S$. More generally, say that we pick our $\zeta_{\mathbf{k}}$’s on an irregular product set, i.e. suppose that for each $j=1\ldots d$ there are “coordinate-frequencies” $\{\zeta_{k}^{j}\}_{k=1}^{2n-1}$ and the multi-frequencies are formed as $\zeta_{\mathbf{k}}=(\zeta^{1}_{k_{1}},\ldots,\zeta^{d}_{k_{d}})$ where $\mathbf{k}\in S=\{1,\ldots,2n-1\}^{d}$. Then, in order for the multidimensional set $\{\zeta_{\mathbf{k}}\}_{\mathbf{k}\in\{1,\ldots,2n-1\}^{d}}$ to be non-degenerate it suffices for each coordinate set $\{\zeta_{k}^{j}\}_{k=1}^{2n-1}$ is non-degenerate (see e.g. Theorem 2 in [Hay82] or Proposition 4.1 in [AnC17]).

Proposition 3.2**.**

Given a non-degenerate set of frequencies $\{\zeta_{\mathbf{k}}\}_{\mathbf{k}\in S}$, the rank of ${\mathcal{F}}$ equals

[TABLE]

Proof.

First assume that $|S|=(2n-1)^{d}$ and denote ${\mathcal{F}}$ in this case by ${\mathcal{F}}_{0}$. As in the previous proof we have that the range of ${\mathcal{F}}_{0}$ is spanned by $\mathbf{k}\mapsto e^{i{\zeta_{\mathbf{k}}\cdot{\bf n}}}$ where ${\bf n}\in\{-n+1,\ldots,n-1\}$. In this case, there are as many different ${\bf n}$’s as there are $\mathbf{k}$’s, and it follows by basic linear algebra that the set of functions of the form $\mathbf{k}\mapsto e^{i{\zeta_{\mathbf{k}}\cdot{\bf n}}}$ is linearly independent if and only if the set of functions ${\bf n}\mapsto e^{i{\zeta_{\mathbf{k}}\cdot{\bf n}}}$ is, which is true by assumption. This finishes the proof under the assumption that $|S|=(2n-1)^{d}$.

Given any particular basis for the domain and the codomain, the operator ${\mathcal{F}}_{0}$ has a matrix representation $\overrightarrow{{\mathcal{F}}_{0}}$ (of dimension $(2n-1)^{d}\times n^{2d}$). We know from what we have already shown that we have $(2n-1)^{d}$ linearly independent columns, i.e. $\overrightarrow{{\mathcal{F}}_{0}}$ has full rank. If $|S|>(2n-1)^{d}$ we can think of $\overrightarrow{{\mathcal{F}}}$ as adding rows to the matrix $\overrightarrow{{\mathcal{F}}_{0}}$, which then impossibly can result in a lower rank. Combined with the previous lemma, this shows that the rank in this case is $(2n-1)^{d}$. Similarly, if $|S|<(2n-1)^{d}$ we can think of $\overrightarrow{{\mathcal{F}}}$ as removing rows from the full rank matrix $\overrightarrow{{\mathcal{F}}_{0}}$, which then clearly results in a new full rank matrix. Hence the rank will be $|S|$, and the proof is complete. ∎

The above proposition is somehow bad news for the "oversampling approach", since it shows that the maximum oversampling factor one could hope for without adding more linearly dependent equations is roughly $2^{d}$. This conclusion is also backed by the numerical experiments in [CES13]: Section 4.4.3 is devoted to numerically demonstrate that oversampling alone does not yield a well-posed inverse problem. Hence Proposition 2.1 can be seen as a theoretical justification of numerical observations in [CES13]. We also validate this conclusion experimentally in Section 5.2. An interesting remark is that two is also a suitable oversampling parameter for other Phase Retrieval methods, see e.g. [MSC98] which backs up this conclusion experimentally. This can be understood since $|\hat{x}(\theta)|^{2}$ is a trigonometric polynomial with $(2n-1)^{d}$ undetermined coefficients, and it is well known that one needs precisely $(2n-1)^{d}$ (non-degenerate) measurements to uniquely determine these coefficients. The following interesting papers [BeP15, BBE17, BrS79, Hay82, HES16, Osh12] contain more information about uniqueness of the Phase Retrieval problem, (as well as other interesting algorithms and plenty of other information). A curious remark is that the factor $2^{d}$ does get better with the dimension $d$, and it is also well known in the experimental community that 3d reconstructions are more stable; see e.g. [CBM06, CZL14]. For the moment, the case $d=3$ is too computationally expensive with the present method, but we hope to address this shortcoming in future work. The proposition also shows that one does not increase stability, or degree of linear independence to be more precise, by picking $\{\zeta_{\mathbf{k}}\}$ from some irregularly sampled grid. Hence we find no motivation to deviate from the simplest possible choice, i.e. picking an $m$ in the range $n\leq m<2n$ and use (10) or, which is the same, setting $S=\{0,\ldots,m-1\}^{d}$ and $\zeta_{\mathbf{k}}=\mathbf{k}/m$ in (11).

3.2 Stabilizing PhaseLift by adding new equations

Section 2.1 of [CES13] provide a list of experimental methods which can be employed in order to add further linearly independent measurements to those given by the operator ${\mathcal{F}}$ from the previous section. In particular it is argued that one can use masks, that is, you create new measurement vectors $f_{\mathbf{k}}^{C}$ by introducing a mask which only lets light through in a region $C$. Mathematically, this amounts to multiplying the image ${\bf x}$ with the characteristic function $\chi_{C}$. If we use regularly sampled measurements as in (10), this gives new measurement-functions via the formula $f_{\mathbf{k}}^{C}({\bf n})=e^{-2\pi i\frac{\mathbf{k}\cdot{\bf n}}{m}}\chi_{C}({\bf n})$ where $C$ is realized as a subset of $\{0,\ldots,n-1\}^{d}$, and $\mathbf{k}$ runs over the index set $\{0,\ldots,m-1\}^{d}$. Also ptychography, optical grating and oblique illuminations are considered. However additional linearly independent measurements are created, we can form an operator ${\mathcal{A}}$ by extending ${\mathcal{F}}$ so that the problems (6)-(8) has an essentially unique solution.

While the methods mentioned above are great from a mathematical perspective, they are often not feasible or slow or expensive from a physical perspective in a concrete experiment. Due to this, a priori estimates on the support of the object measured is still by far the most common method used in practice. This was pioneered by Fienup in [Fie82] and presently the most popular methods to solve the phase retrieval problem in this context are Hybrid Input Output ([Fie82]), Difference Map ([Els03]) or Relaxed Averaged Alternating Reflection ([Luk05]). A simple way to incorporate support constraints into the scheme (6)-(8) is to construct ${\mathcal{A}}$ by extending ${\mathcal{F}}$ by simply adding linear equations $X({\bf m},{\bf n})=0$ for all pairs $({\bf m},{\bf n})\in(\{0,\ldots,n-1\}^{d})^{2}$ such that either ${\bf m}$ or ${\bf n}$ is outside of $C$ (recall that $X$ is a tensor in $(\otimes_{j=1}^{d}{\mathbb{C}}^{n})\otimes(\otimes_{j=1}^{d}{\mathbb{C}}^{n})$). However, this will yield an algorithm which balances meeting the support constraint versus the low rank inducing functional ${\mathcal{Q}}_{\gamma}(\iota_{R^{+}_{K}})$, and hence the output may still be non-zero off $C$, which may be suitable depending on how certain the support estimate is. An alternative is to set ${\mathcal{A}}={\mathcal{F}}$ (i.e. only “pure” Fourier data) and use ADMM on the problem

[TABLE]

which will force the solution to obey the support constraint. We leave it for further research to investigate these options from a practical perspective.

3.3 Estimating $\|{\mathcal{A}}\|$ for masked Fourier data

As explained in section 2, the parameter $\gamma$ in the quadratic envelope needs to be chosen considering the size of $\|{\mathcal{A}}\|^{2}$, and since ${\mathcal{A}}$ in practice is a huge matrix it is not desirable to have to compute its norm. We therefore provide some rough estimates here depending only on the dimension of the images, in the simplest case when $m=n$. We thus assume that ${\mathcal{A}}$ is formed as explained in the first paragraph of the previous section, using $N_{m}$ number of masks (plus pure Fourier data). Let ${\mathcal{A}}_{j}$, $j=0,\ldots,N_{m}$ be the suboperator of ${\mathcal{A}}$ connected to a particular mask (so ${\mathcal{A}}_{0}={\mathcal{F}}$ are those measurement where no mask is used.) In other words ${\mathcal{A}}_{0}(X)(\mathbf{k})=\langle X,f_{\mathbf{k}}\otimes\overline{f_{\mathbf{k}}}\rangle$. The tensors $\{f_{\mathbf{k}}\otimes\overline{f_{\mathbf{k}}}:~{}\mathbf{k}\in\{0,\ldots,n-1\}^{d}\}$ are mutually orthogonal (since the $f_{\mathbf{k}}$’s are). Moreover

[TABLE]

so it follows that ${\mathcal{A}}_{0}$ is $N$ times a unitary operator, which has operator norm 1. Now, if we represent ${\mathcal{A}}_{0}$ as a matrix $\mathbf{A}_{0}$, then it is clear that $\mathbf{A}_{j}$ (for each $j\geq 1$), is obtained from $\mathbf{A}_{0}$ by replacing entire rows and columns by 0. This means that $\|{\mathcal{A}}_{j}\|\leq\|{\mathcal{A}}_{0}\|=N$ by basic estimates, and hence the triangle inequality implies that $\|{\mathcal{A}}\|\leq\sum_{j}\|{\mathcal{A}}_{j}\|\leq(N_{m}+1)N=M$. Summing up we have shown that

[TABLE]

4 Implementation aspects

In this section we show how to efficiently minimize (8) without additional constraints, (or at least how to compute a stationary point). Since (8) is of the type $f+g$ where $f=\mathcal{Q}_{\gamma}(\iota_{R_{1}^{+}})$ is non-convex but $g=\frac{1}{2}\|{\mathcal{A}}(\cdot)-b\|^{2}_{\ell^{2}(S)}$ is smooth, we use the Forward-Backward Splitting scheme, as this has been established to converge to a stationary point under assumptions which are applicable in our setting. Moreover we suggest to use the accelerated version FISTA, since we have observed that this significantly speeds up convergence. The algorithm alternates between a gradient update step and a proximal operator step. The computation of the gradient of $g=\frac{1}{2}\|{\mathcal{A}}(\cdot)-b\|^{2}_{\ell^{2}(S)}$ is mathematically straightforward but very time consuming, and we will therefore introduce certain tricks based on FFT and Toeplitz matrices for efficient evaluation of this step when working with Fourier data, see Section 4.1. The computation of the proximal operator of $\mathcal{Q}_{\gamma}(\iota_{R_{K}^{+}})(X)=\mathcal{Q}_{\gamma}(\iota_{K}^{+})(\lambda(X))$ is rather tricky, the details are given in Section 4.2. We give here a general overview of all the functions involved. For concreteness we present the code when working with Fourier data and a number of masks, and leave it up to the reader to work out details for e.g. ADMM routines with support constraints.

In FISTA, short for Fast Iterative Shrinkage-Thresholding Algorithm, the proximal-gradient steps are taken at the interpolation

[TABLE]

where $\{\theta_{k}\}_{k\geq 1}$ is a sequence of positive real numbers with some properties that ensure convergence in the convex case (see for instance [Bec17]). We used the sequence $\theta_{k}=(k+1)/2$, for $k\geq 1$ as suggested in Remark 10.35 of [Bec17]. Both $X^{k}$ and $X^{k+1}$ are initialized at zero. Upon choosing a step-size $t$ (which we discuss a little further on), the FISTA algorithm now alternates between the update steps

Compute $X_{\text{int}}^{k}$.

• 2.

$X=X^{k}-t\cdot\text{grad}$, where grad is the gradient of $\|{{\mathcal{A}}}(\cdot)-b\|_{\ell^{2}(S)}^{2}$, evaluated at $X^{k}_{\text{int}}$.

• 3.

$X_{k+1}=\text{prox}_{{\mathcal{Q}}_{t\gamma}(\iota_{R_{1}}^{+})}(X)$.

For the latter to be single valued, we need to have $\frac{1}{\gamma}>t$, which puts an upper bound to the choice of step-size $t$.

We now discuss suitable choices of $\gamma$ and $t$. Based on (15) and the theory for the quadratic envelope, it would seem that $\gamma=M^{2}$ would be a natural choice, since it would guarantee that (6) and (8) have the same global minimizer. However, a large $\gamma$ also means that ${\mathcal{Q}}_{\gamma}(\iota_{R_{K}^{+}})$ is, informally speaking, more non-convex, and we have found that the choice $\gamma=N^{2}$ gives a better performance. We have also observed that performance is rather stable with respect to changes in $\gamma$ around this value. With this choice, we get the bound $t<1/N^{2}$ for the step-size. However, the convergence of FISTA is guaranteed (in both convex and non-convex cases) if the stepsize $t$ is $<1/L(\nabla g)$ (see [BCo10] and [ABS13]), where $g$ is the differentiable function and $L(\nabla g)$ the Lipschitz constant of its gradient, which in our case leads to $t\leq 1/\|{\mathcal{A}}\|^{2}$. Based on (15), we therefore used $t=1/(M^{2}+1)$ in our experiments. In our experience, larger step-size leads to faster convergence, but it is important to not exceed the bound.

4.1 Efficient computation of the gradient step

We shall here only treat Fourier measurements with masks, in which case the gradient can be computed very efficiently by using the Fast Fourier Transform (FFT in short), but the details are a bit intricate. A technical problem is the representation of tensors as vectors in a computer. Given any enumeration of the index set $\{0,\ldots,n-1\}^{d}$ and a vector $x\in\otimes_{j=1}^{d}{\mathbb{C}}^{n}$, we denote the corresponding vector by ${\bf x}\in{\mathbb{C}}^{N}$. More precisely, we let $\beta_{n}:\{0,\ldots,n-1\}^{d}\rightarrow\{1,\ldots,N\}$ be a bijection and set ${\bf x}(\beta_{n}({\bf j}))=x({\bf j})$. If for example $d=2$ (so we are dealing with images) and we vectorize by column-stacking, then $\beta_{n}({\bf j})=\beta_{n}(j_{1},j_{2})=j_{2}n+j_{1}+1.$

All vectors are realized as column-vectors, so that e.g. ${\bf x}{\bf x}^{*}$ becomes a matrix. Similarly, if $M$ denotes the amount of elements in $S$ the function $b\in\ell^{2}(S)$ can be identified with a vector ${\bf b}\in{\mathbb{C}}^{M}$ by ordering the elements of $S$. Once both the domain $({\mathbb{C}}^{n})^{d}$ and codomain $\ell^{2}(S)$ have been vectorized, the operators ${\mathcal{A}}$ and ${\mathcal{F}}$ get matrix representations that we denote by A and F. If $\mathbf{a}_{k}$ relates to ${\mathcal{A}}$ (and $\mathbf{A}$) as described shortly before (3), it is now easy to see that

[TABLE]

by basic multivariable calculus.

We first consider the case of no masks and follow the oversampling recommendations from Section 3.1; we pick a parameter $m$ with $n\leq m<2n$ and set $S=\{0,\ldots,m-1\}^{d}$. This set is then in bijective correspondence with $\{1,\ldots,M\}$ through $\beta_{m}$, and we will implicitly identify $k$ in the latter with $\mathbf{k}=\beta_{m}(k)$ in the former. Now $\mathbf{A}=\mathbf{F}$ and $\mathbf{a}_{k}$ derive from pure oscillatory exponentials as in (10). More precisely $\mathbf{a}_{k}$ corresponds to $a_{\mathbf{k}}={f}_{\mathbf{k}}({\bf j})=e^{-2\pi i\frac{{\bf j}\mathbf{k}}{m}}$ as in (10), where $a_{\mathbf{k}}$ is the tensor counterpart of $\mathbf{a}_{k}$.

Similarly, since ${\bf X}$ represents an arbitrary matrix in ${\mathbb{M}}_{N}$ it implicitly defines a tensor $X$ on the tensor product of $\otimes_{j=1}^{d}{\mathbb{C}}^{n}$ with itself, by the formula $X({\bf j},{\bf l})={\bf X}(\beta_{n}({\bf j}),\beta_{n}({\bf l}))$. Summing up we have that $\langle\mathbf{a}_{k}\mathbf{a}_{k}^{*},{\bf X}\rangle_{F}$ equals $\langle f_{\mathbf{k}}\otimes\overline{f_{\mathbf{k}}},X\rangle$ in the space \big{(}\otimes_{j=1}^{d}{\mathbb{C}}^{n}\big{)}\otimes\big{(}\otimes_{j=1}^{d}{\mathbb{C}}^{n}\big{)} or, more explicitly

[TABLE]

where the sums over ${\bf j},{\bf l}$ run over all ${\{0,\ldots,n-1\}^{d}}$, ${\bf p}$ runs through $\{-n+1,\ldots,n-1\}^{d}$ and each coordinate of ${\bf q}$, say $q_{1}$, satisfies $\max(0,-p_{1})\leq q_{1}\leq\min(n-1,n-1-p_{1})$ to ensure that both ${\bf q}$ and ${\bf p}+{\bf q}$ are in ${\{0,\ldots,n-1\}^{d}}$. We now introduce $Y({\bf p})=\sum_{{\bf q}}{X({\bf q},{\bf p}+{\bf q})}$ (where summation bounds for ${\bf q}$ are as before) for ${\bf p}\in\{-n+1,\ldots,n-1\}^{d}$, and define $Y({\bf p})=0$ outside of this grid. In the case $d=1$ the operation above corresponds to summing over lines in $X$ parallel to the diagonal, as in a Toeplitz matrix. For the multivariable case the above is easily computed using ${\bf X}$, by simply noting that $Y({\bf p})=\sum_{{\bf q}}{{\bf X}(\beta_{n}({\bf q}),\beta_{n}({\bf p}+{\bf q}))}$.

Clearly (17) is a sort of inverse Fourier transform on $\overline{Y}$, and moreover the computation of $Y$ is fast; $O(N^{2})$. To make use of FFT for fast evaluation of this expression jointly for all $\mathbf{k}\in\{0,\ldots,m-1\}^{d}$, we introduce $Y_{\textsc{mod}}({\bf p})=\sum_{{\bf n}\in\mathbb{Z}^{d}}f({\bf p}+m{\bf n})$ for all ${\bf p}\in\{0,\ldots,m-1\}^{d}$. Since usually $n\leq m\leq 2m$, note that for each ${\bf p}$ the sum contains at most 2 non-zero entries, so this operation is efficiently evaluated. Summing up we have that

[TABLE]

where $DFT$ denotes the discrete Fourier transform on the grid $\{0,\ldots,m-1\}^{d}$. This can be efficiently evaluated (i.e. $O(N^{2}\log(N)))$ in most programming languages using some preset implementation of the FFT. For example, when $d=3$ and when using MATLAB, we can represent $Y_{\textsc{MOD}}$ as a multidimensional array and evaluate (18) using the command fftn.

Let us now return to the formula (16). We denote the element computed in (18) by $c_{\mathbf{k}}$, and similarly we introduce $b_{\mathbf{k}}=b_{\beta_{m}(k)}$. Then (16) takes the form

[TABLE]

where $\mathbf{f}_{\mathbf{k}}$ is the vector representation of the tensor $f_{\mathbf{k}}$ via $\beta_{n}$, as before. At some fixed index pair $(p,q)\in\{1,\ldots,N\}^{2}$, this is evaluated as

[TABLE]

where ${\bf p}=\beta_{n}(p)$ and ${\bf q}=\beta_{n}(q)$. Clearly all these $N^{2}=n^{2d}$ values reduce to (at most) $M=m^{d}$ different values, and can be computed by an inverse Fourier transform on the tensor $(c_{\mathbf{k}}-b_{\mathbf{k}})$, which has time complexity $O(M^{2}\log(M))$.

This concludes the explanation of how to efficiently compute the gradient step (16) in the case of no masks. Luckily, the incorporation of masks poses very little additional difficulty. Say we have $N_{m}$ masks $C_{1},\ldots,C_{N_{m}}$ as described in Section 3.2. Clearly then the computation of (16) can be separated in $N_{m}+1$ independent pieces, one for each mask. We just describe how to compute one of these contributions for a fixed mask $C$. In this case we again have that $k=1,\ldots,m^{d}$ and we may as well index the vectors using $\mathbf{k}\in\{0,\ldots,m-1\}^{d}$ instead. Now ${a}_{\mathbf{k}}$ corresponds to the tensor $\chi_{C}({\bf j})e^{-2\pi i\frac{{\bf j}\mathbf{k}}{m}}$ and $\mathbf{a}_{\mathbf{k}}$ to its vectorization via $\beta_{n}$. If $\mathbf{w}$ is the vectorization of $\chi_{C}$ via $\beta_{n}$ and $I_{\mathbf{w}}$ denotes the diagonal operator on ${\mathbb{C}}^{N}$ with $\mathbf{w}$ as diagonal, we thus have that $\mathbf{a}_{\mathbf{k}}=I_{\mathbf{w}}\mathbf{f}_{\mathbf{k}}$ with $\mathbf{f}_{\mathbf{k}}$ as before. Returning to the coefficients in (16) we have that

[TABLE]

which means that we can use formula (18) upon replacing ${\bf X}$ with $I_{\mathbf{w}}{\bf X}I_{\mathbf{w}}.$ Finally, since $\mathbf{a}_{\mathbf{k}}\mathbf{a}_{\mathbf{k}}^{*}=I_{\mathbf{w}}\mathbf{f}_{\mathbf{k}}\mathbf{f}_{\mathbf{k}}^{*}I_{\mathbf{w}}$ the summation in (16) can be evaluated as in (19) with the difference that we insert a 0 in the final matrix on positions $(p,q)$ where either of the two coordinates ${\bf p}$ and ${\bf q}$ lie outside of $C$.

4.2 Computation of the proximal operator

The expression of the proximal operator of $\mathcal{Q}_{\gamma}(\iota_{K}^{+})$ is quite involved. We first recall that $\mathcal{Q}_{\gamma}(f)=\mathcal{S}_{\gamma}(\mathcal{S}_{\gamma}(f))$ where $\mathcal{S}_{\gamma}(h)(x)$ is computed by first taking the Legendre transform of $h+\frac{\gamma}{2}\|\cdot\|^{2}$ and then subtracting $\frac{\gamma}{2}\|\cdot\|^{2}$. By an alteration of the Moreau decomposition (cfr. Chapter 14 of [BCo10]) it is possible to show (cfr. [Car16], Proposition 3.3) that if $\rho>\gamma$ we have

[TABLE]

where (see again [Car16])

[TABLE]

being $\hat{{\bf x}}$ the decreasing rearrangement of ${\bf x}$. We now sketch the routine for computing $\text{prox}_{\frac{\rho-\gamma}{\rho\gamma}\mathcal{S}_{\gamma}(\iota_{K}^{+})}({\bf y})$. After some basic algebra it turns out that

[TABLE]

By the rearrangement inequality we may assume that ${\bf y}$ is already ordered non increasingly and therefore the latter turns into

[TABLE]

For a given ${\bf x}$ let $\tilde{k}({\bf x})$ be the minimum between $K$ and the last index $j$ for which $x_{j}$ is non-negative. Then (21) becomes

[TABLE]

Now it is clear that if the vector

[TABLE]

is non-increasing then it is a global minimizer for (21). In particular this happens whenever ${\bf y}$ switches sign before $K$ (i.e. $\tilde{k}({\bf y})<K$) or whenever $y_{K}\geq\rho y_{K+1}/\gamma$. We exclude these two scenarios from the further analysis.

If $\tilde{k}({\bf x})<K$ then $x_{\tilde{k}({\bf x})}\geq 0$ and $x_{\tilde{k}({\bf x})+1}<0$, but looking at (22) we see that setting $x_{\tilde{k}({\bf x})+1}=0$ would diminish the quantity (since $y_{\tilde{k}({\bf x})+1}>0$), a contradiction. Thus we have that $\tilde{k}({\bf x})=K$, and so (22) becomes

[TABLE]

Let ${\bf x}$ be some candidate for the global minimum and set $s:=x_{K}$. It is easy to see that ${\bf x}$ must have the following structure:

[TABLE]

Now consider (25) inserted into (24); except from an additive constant this function takes the form

[TABLE]

By inspection it is clear that this function is convex and has a unique minimum on the interval $V=[y_{K},\rho y_{K+1}/\gamma]$; indeed the first sum is constant on $(-\infty,y_{K}]$ and strictly convex (increasing) on $(y_{K},\infty)$, whereas the second sum is strictly convex (decreasing) $(-\infty,\rho y_{K+1}/\gamma)$ and constant on $[\rho y_{K+1}/\gamma,\infty)$. We seek for the unique solution of $\frac{d}{ds}F(s)=0$ in $V$. Let $j^{*}$ be the smallest index such that $y_{j^{*}}\leq\rho y_{K+1}/\gamma$ and let $l^{*}$ be the biggest index such that $\rho y_{l^{*}}/\gamma\geq y_{K}$; now the set $\{y_{i}\}_{i=j^{*}}^{K}\cup\{\rho y_{i}/\gamma\}_{i=K+1}^{l^{*}}$ provides a partition of $V$. Let $I$ be one of the subintervals. For all values of $s$ in $I$ let $j$ be the first index such that $x_{j}=s$ and let $l$ be the last, where $x_{j}$ is given by (25). The formula for the solution of $\frac{d}{ds}F(s)=0$ becomes

[TABLE]

If this value is outside $I$ then the optimal $s$ is to be found in another interval. By strict convexity $s_{I}\in I$ will hold in at least one interval, and at most two intervals (when $s_{I}$ lies on the boundary).

In conclusion we can summarize the previous observations in an algorithm for computing the proximal operator:

If $y_{K}<0$ or $\rho y_{K+1}/\gamma<y_{K}$ return (23), else

• 2.

compute $j^{*}$ and $l^{*}$;

• 3.

sort $\{y_{i}\}_{i=j^{*}}^{K}\cup\{\rho y_{i}/\gamma\}_{i=K+1}^{l^{*}}$ decreasingly and call it ${\bf z}$;

• 4.

for $m\in\{1,\dots,l^{*}-j^{*}-1\}$ set $s=(z_{m}+z_{m+1})/2$ and compute the indices $j$ and $l$ as described above. Notice that the indices $j$ and $l$ are the same for all $t\in[z_{m},z_{m+1}]$, and this is why it is enough to consider the midpoints;

• 5.

compute the new value of $s$ according to (26);

• 6.

if $z_{m}\geq s\geq z_{m+1}$ stop and return (25). Otherwise increase $m$ with $1$ and repeat the steps 4-5.

The “for-loop” 4-6 must stop since $F$ is strictly convex and has a unique minimizer.

5 Numerical examples

We will illustrate the various results by conducting four experiments. In the first experiment we will compare the proposed method with PhaseLift and reweighted PhaseLift, in the second we test the oversampling conjectures from Section 3, in the third we discuss the capacity of the proposed method of solving low rank PSD problems of any given rank as predicted in Section 2, and finally we end with the reconstruction of a 2D image using the tricks of Section 4.

5.1 1D synthetic "masked" Fourier measurements

In this subsection we consider 1-dimensional signals ${\bf x}_{0}\in\mathbb{C}^{100}$ with $\|{\bf x}_{0}\|=1$ randomly generated with a Gaussian distribution; they represent our "ground truth", i.e. the signals that we are interested in retrieving. The method we proposed is compared to the reweighted nuclear norm technique, following the numerical section of [CES13]. The idea of reweighting the nuclear norm was heuristically introduced by [FHB01]. Despite of the lack of a systematic mathematical theory behind, this trick seemed to work quite well for our problem too, proving to be able to significantly enhance the capacity of the nuclear norm of finding low rank matrices.

We consider binary masks, as introduced in Section 3.2. Following [CES13] we use $m=n$ (oversampling is considered separately in the next section) and three masks, so $n=m=N=100$ and $M=400$. Note that this means that we use 400 equations plus a rank-PSD constraint to determine roughly 5000 variables in the lifted problem. We focus our attention in comparing the approximated solutions to the two problems (8) and

[TABLE]

where $W^{(l)}$ are weight (diagonal) matrices and where ${\mathcal{A}}$ is as described in section 3.2 and ${{\bf b}}={\mathcal{A}}({\bf x}_{0}{\bf x}_{0}^{*})+\boldsymbol{\epsilon}$ with $\epsilon\in\mathbb{C}^{400}$ additive random gaussian noise.

The proposed method only contains one user parameter $\gamma$ which is fixed as described in Section 4, independent of noise level. The minimization of (27) on the other hand relies on $\lambda$ which must scale with the noise for optimal performance, and moreover the update of the weights requires yet another parameter $\delta$ to be chosen; following [CWB08] we set

[TABLE]

being $X^{(l)}$ the outcome of the algorithm after some fixed number of iterations with $W^{(l)}$. $W^{(1)}$ is just the identity matrix.

We have not found any clear guidance in the literature for how to pick $\lambda$ and $\delta$; in [CES13] the algorithm is run several times for different $\lambda$’s picked according to a golden section search or cross validation. This is very time consuming, especially taking into account that the algorithm is slow to begin with due to the lifting trick. We have found that $\lambda=0.01+0.75\|\epsilon\|$ and $\delta=0.01+0.05\|\epsilon\|$ seemed to work fairly well in our experiments, and these are the values used in the graphs presented. The bisection search only gives a marginal improvement to the results shown, but it would be unfair to do this in comparison with 8, since this is run only once with one fixed choice of $\gamma$. For practical purposes, the fact that (8) does not rely on intricate parameter choices is clearly a strength.

To compute the numerical approximations we used the FISTA algorithm (already previously described) with $10000$ iterations for (8) and with $10000$ iterations for each weights update for (27); we found out that $l>2$ does not give any significant contribution to reconstruction accuracy. The stepsize was fixed to $t=1/(16N^{2}+1)=1/1601$ (see section 4) and 5 different trials for each level of noise were carried out, in the sense that for each $d\in\{0,0.5,1,1.5,2,2.5,3,3.5,4\}$ we run five different experiments, i.e. we randomly generated five different ground truths $X_{0}$ (and consequently five different noise vectors $\epsilon$ with $d=\|\epsilon\|$), and we averaged the various outcomes over the number of trials. All the masks were randomly generated and the norm of the measured data vector is averagely $\approx 20$; in the graphs below we cover cases where the noise has up to $20\%$ of the signal magnitude. According to [CES13] the approximated signal reconstructed using (27) is chosen as the eigenvector of the greatest eigenvalue of the output $\tilde{X}$ of the algorithm. Since in (8) only rank 1 matrices are involved, the approximated signal will be the only nonzero eigenvector of $\tilde{X}$. The underlying signal is unique only up to a global phase, therefore the distance between ${\bf x}_{0}$ and $\tilde{{\bf x}}$ is computed as

[TABLE]

In Figure 3 we plot $D$, $\|\tilde{X}-X_{0}\|_{F}$ and in Figure 4 the average rank of $\tilde{X}$. The graphs show that the performances of the three methods. It is clear that nuclear norm without reweighting is inferior in all aspects. We focus therefore on discussing the proposed method versus reweighted nuclear norm. In terms of “norms precision”, these are essentially comparable; nevertheless the reweighted nuclear norm fails in finding the correct rank for high levels of noise. Thus, when this method is equivalent with the proposed method, it uses more degrees of freedom. On the other hand the method that we propose is always able to retrieve rank 1 solutions, in perfect compliance with the theory developed.

Since determining the numerical rank is somewhat fishy, we include the following tables which show the first ten eigenvalues (decreasingly ordered) of the approximation $\tilde{X}$, computed respectively with (8) and (27), in each of the five different experiments run for the noise level $\|\epsilon\|=3$. As is plain to see, the solutions obtained by the reweighted nuclear norm are far from rank 1.

5.2 Stabilizing by oversampling

In broad terms, the conclusion of Section 3 is that it is desirable to oversample a factor 2 (in each dimension), but beyond this point no further gain is expected. We test this for the one dimensional case on synthetic data of size $N=25$ with an identical setup to the one above. As measurement we use ${\bf b}={\mathcal{A}}({\bf X}_{0})+\boldsymbol{\epsilon}$ where $\boldsymbol{\epsilon}$ is noise which we vary between 0 and $20\%$ of the magnitude of ${\bf b}$. Each data point in the below graphs is the average of five distinct trials.

A numerical observation in [CES13] is that one needs 3 masks to have stable inversion. The authors also point out that a low residual error $\|{\mathcal{A}}(\tilde{\bf X})-{\bf b}\|$ in combination with a poor actual error $\|\tilde{\bf X}-{\bf X}_{0}\|$ is an indicator of having an ill-posed inverse problem, i.e. one where several different points minimize the functional in question. Below we plot both graphs as a function of the Noise to Signal Ratio $\|\boldsymbol{\epsilon}\|/\|{\bf b}\|$.

To the left we see that the residual errors are almost identical, and moreover close to the line $y=x$, as expected, since it is clearly unlikely to beat the noise level by any notable margin. Indeed, if we were to get $\tilde{\bf X}={\bf X}_{0}$ then clearly $\|{\mathcal{A}}(\tilde{\bf X})-{\bf b}\|=\|\boldsymbol{\epsilon}\|$, so the above graph actually tells us that we most of the time find a solution better than ground truth.

To the right we see the normalized reconstruction error, and we can confirm that the expectations from Section 3 is in perfect compliance with the outcome. Oversampling by a factor two has a significant effect in improving the actual error, whereas oversampling by a factor 3 only has a marginal effect, as predicted by Proposition 3.2. Noteworthy is that also the latter two curves are close to “$y=x$”.

Encouraged by this result we now try the same setup with 2 masks. As is plain to see, the pattern repeats itself with the major difference that $m=n$ is completely unreliable except for 0 noise, whereas $m=2n$ and $m=3n$ do a similar job which, although not extremely convincing, clearly outperforms $m=n$. In conclusion, it seems that the recommendation of oversampling with a factor of two has both theoretical and numerical support, at least in the case $d=1$.

5.3 Finding the correct rank for general PSD estimation problems

We have ran extensive tests for minimizing (8) with different choices of ${\mathcal{A}}$, ${\bf b}$ and ranks $K$, and never found a single instance where the algorithm converges to a point with the wrong rank, i.e. different from $K$, as long as $\gamma>\|A\|^{2}$.222In practice it may be beneficial to violate this restriction, as explained in Section 4 This is in accordance with the theory for quadratic envelopes as developed in [Car19] and summarized in Section 2, stating that local minimizers of (8) form a subset of those for (6), except for the fact that FBS is only guaranteed to find stationary points,333At least relying on the theory summarized in Section 4 not necessarily local minima. This indicates that either there are no saddle point type stationary points, or that FBS actually has the capacity to avoid them. We do not underline these observations with any specific graph, since it seems impossible to design one experiment that would cover all different types of potential applications.

5.4 2D image reconstruction

We complement this section by displaying a 2D reconstruction example. The example chosen is a $27\times 27$ pixels grayscale real image; the initial measurements were contaminated with $1\%$ of additive gaussian random noise. The algorithm used was again FISTA with eight binary masks and the parameters choice indicated in Section 3.3 where we employed the tricks of Section 4.1 for efficient evaluation via FFT of the gradient-step. The reconstruction was made on a standard laptop and the bottleneck timewise for the algorithm is the eigendecomposition, which prohibits larger images to be processed. We plan to address this shortcoming in future works, which clearly needs to be solved for the algorithm (as well as PhaseLift) to be of practical use for 2D or 3D imaging.

Our experience suggested that the 2D problem needs a higher number of masks (with respect to the 1D case) in order to be stabilized; this observation is present in [CES13] as well where, in the noisefree case, eight binary masks are used. We used $m=n$ for simplicity, it is of course possible that oversampling could lead to fewer masks, but we have not yet tested this.

On the side, the original image is the second while the first is our reconstruction.