A Partial Derandomization of PhaseLift using Spherical Designs

Introduction

It is by no means clear how many such amplitude measurements are necessary to allow for recovery. Thus from the very beginning, there have been a number of works regarding injectivity conditions for this problem in the context of the specific applications .

Balan et al. consider the scenario of $\mathcal{O}(d^{2})$ measurements, which form a complex projective $2$-design (cf. Def. 3 below). They derive an explicit reconstruction formula for this setup based on the following observation well known in conic programming. Namely, the quadratic constraints on $x$ are linear in the outer product $xx^{*}$:

This “lifts” the problem to matrix space of dimension $d^{2}$, where it becomes linear and can be explicitly solved to find the unique solution.

As we will show in Theorem 2, it is, without making additional assumptions on the $2$-design, not possible to use as measurements a random subset of this $2$-design which is of size $o(d^{2})$. In other words, for the measurement scenario described in , the quadratic scaling in $d$ is basically unavoidable.

To contrast these two extreme approaches, ref. works with a number of measurements close to the absolute minimum, but there are no tractable reconstruction schemes provided, the question of numerical stability is not considered, and it is unclear whether non-generic measurements – i.e., vectors with additional structural properties – can be employed. On the other hand, the number of measurements in is much larger, while the measurements are highly structured and there is an explicit reconstruction method. A number of recent works including this paper aim to balance between these two approaches, working with a number of measurements only slightly larger while having at least some of the desired properties mentioned above.

Ref. introduces a reconstruction method called polarization that works for $\mathcal{O}(d\log d)$ measurements and can handle structured measurement vectors, including the masked illumination setup that appears in diffraction imaging , where the measurements are generated by the discrete Fourier transform preceded by a random diagonal matrix. For Gaussian measurements, the polarization approach has also shown to be stable with respect to measurement noise . While simulations seem to suggest stability also for the derandomized masked illumination setup, a proof of stability is – to our knowledge – not available yet.

An alternative approach, which we will also follow in this paper, is the PhaseLift algorithm, which is based on the lifted formulation (1). The algorithm was introduced in and reconstruction guarantees have been provided in . The central observation is that the matrix $xx^{*}$, while unknown, is certainly of rank one. This connects the phase retrievel problem with the young but already extensive field of low-rank matrix recovery . Over the past years, this research program has rigorously identified many instances in which low-rank matrices can be efficiently reconstructed from few linear measurements. The existing results on low-rank matrix recovery were not directly applicable to phase retrieval, because the measurement matrices $a_{i}a_{i}^{*}$ failed to be sufficiently incoherent in the sense of (the incoherence parameter captures the well-posedness of a low-rank recovery problem). For the case of Gaussian measurement vectors $a_{i}$, Candès, Strohmer, Voroninski and Li were able to circumvent this problem, providing problem-specific stable recovery guarantees for a number of measurements of optimal order $\mathcal{O}(d)$. For recovery, they use a convex relaxation of the rank minimization problem, which makes the reconstruction algorithm tractable.

It should be noted, however, that because of the significantly increased problem dimensions, PhaseLift is not as efficient as many phase retrieval algorithms developed over the last decades in the physics literature (such as ) and the optimization literature (for example ). Recently there have been attempts to provide recovery guarantees for alternating minimization algorithms , which are somewhat closer to the algorithms used in practice, but this direction of research is only at its beginnings.

While the above mentioned recovery guarantees for PhaseLift address the issues of tractable reconstruction and stability with respect to noise, these results leave open the question of whether measurement systems with additional structure and less randomness still allow for guaranteed recovery. There are both practical and theoretical motivations for pursuing such generalizations: A practitioner may be constrained in the choice of measurements by the application at hand or reduce the amount of randomness required for implementation purposes. The most prominent example are again masked Fourier measurements, which appear as a natural model in diffraction imaging, but a lot of different scenarios imposing different structure are conceivable. From a theoretical point of view, the use of Gaussian vectors obscures the specific properties that make phase retrieval possible. As discussed in the following subsection, it is a common thread in randomized signal processing that results are first established for Gaussian measurements and later generalized to structured ensembles.

In this paper, we focus on the theoretical aspect: which properties of a measurements are sufficient for PhaseLift to succeed? We prove recovery guarantees for ensembles of measurement vectors drawn at random from a finite set whose first $2t$ moments agree with those of Haar-random vectors (or, essentially, Gaussian vectors). A configuration of finite vectors which gives rise to such an ensemble is known as a complex projective $t$-design The definition of a $t$-design varies between authors. In particular, what is called a $t$-design here (and in most of the physics literature), would sometimes be referred to as a $2t$ or even a $(2t+1)$-design. See Section 3.3 for our precise definition. . Designs were introduced by Delsarte, Goethals and Seidel in a seminal paper and have been studied in algebraic combinatorics , coding theory , and recently in quantum information theory . Furthermore, complex projective $2$-designs were the key ingredient for the reconstruction formula for phase retrieval proposed in .

One may see a more general philosophy behind this approach. In the field of sparse and low-rank reconstruction, a number of recovery results had first been established for Gaussian measurements. In subsequent works, it has then been proven that measurements drawn at random from certain fixed orthonormal bases are actually sufficient. Examples include uniform recovery guarantees for compressed sensing ( vs. ) and low-rank matrix recovery ( vs. ), respectively. Typically, the de-randomized proofs require much higher technical efforts and deliver slightly weaker results. For a recent survey on structured random measurements in signal processing see .

As the number of measurements needed for phase retrieval is larger than the signal space dimension, one cannot expect these results to exactly carry over to the phase retrieval setting. Nevertheless, the question remains whether there is a larger, but preferably not too large, set such that measurements drawn from it uniformly at random allow for phase retrieval reconstruction guarantees. In some sense, the sampling scenario we seek can be interpreted as an interpolation between the maximally random setup of Gaussian measurement with an optimal order of measurements and the construction in , which is completely deterministic, but suboptimal in terms of the embedding dimension. While in this paper, we will focus on the phase retrieval problem, we remark that such an interpolating approach between measurements drawn from a basis and maximally random measurements may also be of interest in other situations where constructions from bases are known, but lead to somewhat suboptimal embedding dimensions.

The concept of $t$-designs, as defined in Section 3.3, provides such an interpolation. The intuition behind that definition is that with growing $t$, more and more moments of the random vector corresponding to a random selection from the $t$-design agree with the Haar measure on the unit sphere. In that sense, as $t$ scales up further, $t$-designs give better and better approximations to Haar-random vectors.

From a practical point of view, the usefulness of these concepts hinges on the availability of constructions for designs. Explicit constructions for any order $t$ and any dimension $d$ are known – however, they are typically “inefficient” in the sense that they require a vector set of exponential size. For example, the construction in uses $\mathcal{O}(t)^{d}$ vectors which is exponential in the dimension $d$.

Tighter analytic expressions for exact designs are notoriously difficult to find. Designs of degree 2 are widely known . A concrete example is used for the converse bound in Section 7 (as well as for the converse bounds for low-rank matrix recovery from Fourier-type bases in ). For degree 3, both real While stated only for dimensions that are a power of $2$, the results can be used for construtions in arbitrary dimensions . and complex designs are known. For higher $t$, there are numerical methods based on the notion of the frame potential , non-constructive existence proofs , and constructions in sporadic dimensions (c.f. and references thererin).

Importantly, almost-tight randomized constructions for approximate designs for arbitrary degrees and dimensions are known . The simplest results show that collections of Haar-random vectors form approximate $t$-designs. This indeed can reduce randomness: One only needs to expend a considerable amount of randomness once to generated a design – for subsequent applications it is sufficient to sample small subsets from it The situation is comparable to the use of random graphs as randomness expanders .. Going further, there have been recent deep results on designs obtained from certain structured ensembles . We do not describe the details here, as they are geared toward quantum problems and may have to be substantially modified to be applicable to the phase retrivial. The only connection to phase retrieval to date is the estimation of pure quantum states .

Finally we point out that the notion of the frame potential above is no coincidence. In a frame-theoretic approach to designs is provided, underlining their close connection.

2. Main results

In this paper, we show that spherical designs can indeed be used to partially derandomize recovery guarantees for underdetermined estimation problems; we generalize the recovery guarantee in to measurements drawn uniformly at random from complex projective designs, at the cost of a slightly higher number of measurements.

Here $\omega\geq 1$ is an arbitrary parameter and $C$ is a universal constant.

As the discussion of the previous subsection suggests, the bounds on the sampling rate decrease as the order of the design increases. For fixed $t$, and up to poly-log factors, it is proportional to $\mathcal{O}(d^{1+2/t})$. This is sub-quadratic for the regime $t\geq 3$ where our arguments apply. If the degree is allowed to grow logarithmialy with the dimension (as $t=2\log d$), we recover an optimal, linear scaling up to a polylog overhead, $m=\mathcal{O}(d\,\log^{3}d)$.

In light of the highly structured, analytical and exact designs known for degree 2 and 3, it is of great interest to ask whether a linear scaling can already be achieved for some small, fixed $t$. As shown by the following theorem, however, for $t=2$ not even a subquadratic scaling is possible if no additional assumptions are made, irrespective of the reconstruction algorithm used.

Suppose that $m$ measurement vectors $y_{1},\dots,y_{m}$ are sampled independently and uniformly at random from $D_{2}$. Then, for any $\omega\geq 0$, the number of measurements must obey

It is worthwhile to put this statement in perspective with other advances in the field. Throughout our work, we have only demanded that the set of all possible measurement vectors forms a $t$-design and have not made any further assumptions. Theorem 2 has to be interpreted in this regard: The 2-design property alone does not allow for a sub-quadratic scaling when a “reasonably small” probability of failure is required in the recovery process.

Note that this does not exclude the possibility that certain realizations of 2-designs can perform better, if additional structural properties can be exploited. A good example for such a measurement process is the multi-illumination setup provided in . In the authors of this paper verified that the set of all measurement vectors used in the framework of does constitute a 2-design (Lemma 6). Additional structural properties – most notably a certain correlated Fourier basis structure in the individual measurements – allowed for establishing recovery guarantees already for $m=\mathcal{O}\left(d\log^{4}d\right)$ measurements and $m=\mathcal{O}(d\log^{2}d)$ , respectively – which both clearly are sub-quadratic sampling rates.

3. Outlook

There are a number of problems left open by our analysis. First, recall that our results achieve linear scaling up to logarithmic factors only when samples are drawn from a set of superpolynomial size. Thus it would be very interesting to find out whether there are polynomial size sets such that sampling from them achieves such a scaling, in particular, if $t$-designs for some fixed $t$ can be used. The case of $t=3$ seems particularly important in that regard, since the converse bound (Theorem 2) shows that a design order of at least 3 is necessary. Also, highly structued 3-designs are known to exist (see above).

Another important follow-up problem concerns approximate $t$-designs. While our main result is phrased for exact $t$-designs, certain scenarios will only exhibit approximate design properties. We expect that our proofs can be generalized to such a setup, but also leave this problem for future work. Lastly, the reconstruction quality for noisy measurements is also an important issue yet to be investigated.

Numerical Experiments

In this section we complement our theoretical results with numerical experiments, which we have implemented in Matlab using CVX . As may have been expected, these experiments suggest that PhaseLift from designs actually works much better than our main theorem suggests. To be concrete, we use stabilizer states – a highly structured vector set which is very prominent in quantum information theory . Stabilizer states exist in any dimension, though their properties are somewhat better-behaved in prime power dimensions. In this case, there exists $\mathcal{O}(d^{\log d})$ stabilizer state vectors. Due to their rich combinatorial structure, these vectors can be constructed efficiently. For dimensions $d=2^{n}$ that are a power of two, it is known that the set of stabilizer states forms a 3-design. This statement is false for other prime power dimensions ($d\neq 2^{n}$ for some $n$), where they only form an exact 2-design. However, weighted 3-designs can be constructed for arbitrary dimensions $d$ by projecting down stabilizer states from the next largest power-of-2-dimension $2^{n}$ obeying $2^{n-1}~{}<d<2^{n}$ . For further clarification of the concept of exact and weighted $t$-designs we defer the reader to and references therin.

We obtain the picture of a relatively sharp phase transition along a line that scales linearly in the problem dimension. In fact, the transtion seems to occur in the vicinity of the line $m=4d-4$ – drawn in red in Figure 1. This seems to agree with the conjecture that $4d-4$ measurements are required for injectivity (see e.g. ). However, there are a few differences in the problem setup: Firstly, the conjecture only asks whether there is a unique solution, while the numerical simulations study whether the PhaseLift algorithm can find it. Secondly, the conjecture concerns unique solutions for all possible inputs, while numerically, we estimate the success probability. And thirdly, the conjecture states that generic measurements work, while our simulations use a specific random procedure (drawn uniformly from a 3-design) to generate them.

Technical Background and Notation

In this work we require three different objects of linear algebra: vectors, matrices and operators acting on matrices.

We will work with vectors in a $d$-dimensional complex Hilbert space $V^{d}$ equipped with an inner product $\langle\cdot,\cdot\rangle$. We refer to the associated induced norm by

We will denote such vectors by latin characters. For $z\in V^{d}$, we define the dual vector $z^{*}\in(V^{d})^{*}$ via

On the level of matrices we will exclusively consider $d\times d$ dimensional hermitian matrices, which we denote by capital latin characters. Endowed with the Hilbert-Schmitt (or Frobenius) scalar product

the space $H^{d}$ becomes a Hilbert space. In addition to that, we will require the 3 different Schatten-norms

where the second one is induced by the scalar product (3). These three norms are related via the inequalities

We call a hermitian matrix $Z$ positive-semidefinite ($Z\geq 0$), if $\langle y,Zy\rangle\geq 0$ for all $y\in V^{d}$. Positive semidefinite matrices form a cone (Chapter II,12), which induces a partial ordering of matrices. Concretely, for $Z,Y\in H^{d}$ we write $Y\geq Z$ if $Y-Z$ is positive-semidefinite ($Y-Z\geq 0$).

Finally, we will frequently encouter matrix-valued operators acting on the space $H^{d}$. We label such objects with capital caligraphic letters and introduce the operator norm

induced by the Frobenius norm on $H^{d}$. It turns out that only very few matrix-valued operators will appear below. These are: the identity map

and (scalar multiples of) projectors onto some matrix $Y\in H^{d}$. The latter corresponds to

The notion of positive-semidefiniteness directly translates to matrix valued operators. Concretely, we call $\mathcal{M}$ positive-semidefinite ($\mathcal{M}\geq 0$) if $(Z,\mathcal{M}Z)\geq 0$ for all $Z\in H^{d}$. Again, this induces a partial ordering. Like in the matrix case, we write $\mathcal{N}\geq\mathcal{M}$, if $\mathcal{N}-\mathcal{M}\geq 0$. It is easy to check that all the operators introduced so far are positive semidefinite and in particular we obtain the ordering

2. Multilinear Algebra

The properties of $t$-designs are most naturally stated in the framework of ($t$-fold) tensor product spaces. This motivates recapitulating some basic concepts of multilinear algebra that are going to greatly simplify our analysis later on. The concepts presented here are standard and can be found in any textbook on multilinear algebra. Our presentation has been influenced in particular by .

Let $V_{1},\ldots,V_{k}$ be (finite dimensional, complex) vector spaces, and let $V_{1}^{*},\ldots,V_{k}^{*}$ be their dual spaces. A function

is multilinear, if it is linear in each $V_{i}$, $i=1,\ldots,k$. We denote the space of such functions by $V_{1}^{*}\otimes\cdots\otimes V_{k}^{*}$ and call it the tensor product of $V_{1}^{*},\ldots,V_{k}^{*}$. Consequently, the tensor product $\left(V^{d}\right)^{\otimes k}=\bigotimes_{i=1}^{k}V^{d}$ is the space of all multilinear functions

and we call the elementary elements $z_{1}\otimes\cdots\otimes z_{k}$ the tensor product of the vectors $z_{1},\ldots,z_{k}\in V^{d}$. Such an element can alternatively be defined more concretely via the Kronecker product of the individual vectors. However, such a construction requires an explicit choice of basis in $V^{d}$ which is not the case in (6).

With this notation, the space of linear maps $V^{d}\to V^{d}$ ($d\times d$-matrices) corresponds to the tensor product $M^{d}:=V^{d}\otimes\left(V^{d}\right)^{*}$ which is spanned by $\left\{y\otimes z^{*}:\;y,z\in V^{d}\right\}$ – the set of all rank-1 matrices. For this generating set of $M^{d}$, we define the trace to be the natural bilinear map

for all $y,z\in V^{d}$. The familiar notion of trace is obtained by extending this definition linearly to $M^{d}$.

Using $M^{d}=V^{d}\otimes\left(V^{d}\right)^{*}$ allows us to define the (matrix) tensor product $\left(M^{d}\right)^{\otimes k}$ to be the space of all multilinear functions

in complete analogy to the above. We call the elements $Z_{1}\otimes\cdots\otimes Z_{k}$ the tensor product of the matrices $Z_{1},\cdots,Z_{k}\in M^{d}$.

On this tensor space, we define the partial trace (over the $i$-th system) to be

Note that with the identification $M^{d}=V^{d}\otimes\left(V^{d}\right)^{*}$, $\operatorname{tr}_{i}$ corresponds to the natural contraction at position $i$. The partial trace over more than one system can be obtained by concatenating individual traces of this form, e.g. for $1\leq i<j\leq k$

In particular, the full trace then corresponds to

Let us now return to the tensor space $\left(V^{d}\right)^{\otimes k}$ of vectors. We define the (symmetrizer) map $P_{\operatorname{Sym}^{k}}:\left(V^{d}\right)^{\otimes k}\to\left(V^{d}\right)^{\otimes k}$ via their action on elementary elements:

where $S_{k}$ denotes the group of permutations of $k$ elements. This map projects $\left(V^{d}\right)^{\otimes k}$ onto the totally symmetric subspace $\operatorname{Sym}^{k}$ of $\left(V^{d}\right)^{\otimes k}$ whose dimension is

3. Complex projective designs

The idea of (real) spherical designs originates in coding theory and has been extended to more general spaces in . We refer the interested reader to Levenshtein for a unified treatment of designs in general metric spaces and from now on focus on designs in the complex vector space $V^{d}$.

Roughly speaking, a complex projective $t$-design is a finite subset of the complex unit sphere in $V^{d}$ with the property that the discrete average of any polynomial of degree $t$ or less equals its uniform average. Many equivalent definitions – see e.g. – capture this essence. However, there is a more explicit definition of a $t$-design that is much more suitable for our purpose:

A finite set $\{w_{1},\ldots,w_{N}\}\subset V^{d}$ of normalized vectors is called a $t$-design of dimension $d$ if and only if

where $P_{\operatorname{Sym}^{t}}$ denotes the projector onto the totally symmetric subspace (7) of $(V^{d})^{\otimes t}$ and consequently $\dim(\operatorname{Sym}^{t})=\binom{d+t-1}{t}$.

where the right hand side is integrated with respect to the Haar measure. This form makes the statement that $t$-designs mimic the first $2t$ moments of Haar measure more explicit.

P. Seymor and T. Zaslavsky proved in that $t$-designs on $V^{d}$ exist for every $t,d\geq 1$, provided that $N$ is large enough ($N\geq N(d,t)$), but they do not give an explicit construction. A necessary criterion – cf. – for the $t$-design property is that the number of vectors $N$ obeys

However, the proof in is non-constructive and known constructions are “inneficient” in the sense that the number of vectors required greatly exceeds (10). Hayashi et al. proposed a construction requiring $\mathcal{O}(t)^{d}$ vectors. For real spherical designs other “inefficient” constructions have been proposed ($N=t^{\mathcal{O}(d^{2})}$) which can be used to obtain complex projective designs.

Adressing this apparant lack of efficient constructions, Ambainis and Emerson proposed the notion of approximate desings. These vector sets only fulfill property (9) only up to an $\epsilon$-precision, but their great advantage is that they can be constructed efficiently. Concretely, they show that for every $d\geq 2t$, there exists an $\epsilon=\mathcal{O}(d^{-1/3})$ approximate $t$-design consisting of $\mathcal{O}(d^{3t})$ vectors only.

The great value of $t$-designs is due to the following fact: If we sample $m$ vectors $a_{i},\ldots,a_{m}$ iid from a $t$-design $D_{t}=\left\{w_{1},\ldots,w_{N}\right\}$, the design property guarantees (with $A_{i}=a_{i}a_{i}^{*}$ and $W_{i}=w_{i}w_{i}^{*}$)

for all $1\leq k\leq t$. This knowledge about the first $t$ moments of the sampling procedure is the key ingredient for our partial derandomization of Gaussian PhaseLift .

4. Large Deviation Bounds

This approach makes heavy use of operator-valued large deviation bounds. They have been established first in the field of quantum information by Ahlswede and Winter . Later the first author of this paper and his coworkers successfully applied these methods to the problem of low rank matrix recovery . By now these methods are widely used and we borrow them in their most recent (and convenient) form from Tropp .

5. Wiring Diagrams

The defining property (9) of $t$-designs is phrased in terms of tensor spaces. To work with these notions practically, we need tools for efficiently computing contractions between high-order tensors. The concept of wiring diagrams provides such a method – see for an introduction and also (however, they use a slightly different notation). Here, we give a brief description that should suffice for our calculations.

Roughly, the calculus of wiring diagrams associates with every tensor a box, and with every index of that tensor a line emanating from the box. Two connected lines represent contracted indices. (More precisely, we place contravariant indices of a tensor on top of the associated box and covariant ones at the bottom. However, one should be able to digest our calculations without reference to this detail). A matrix $A:V^{d}\to V^{d}$ can be seen as a two-indexed tensor ${A^{i}}_{j}$. It will thus be represented by a node $A$ with the upper line corresponding to the index $i$ and the lower one to $j$. Two matrices $A,B$ are multiplied by contracting $B$’s “contravariant” index with $A$’s “covariant” one:

corresponds to a contraction of the two indices of a matrix:

Tensor products are arranged in parallel:

Hence, a partial trace takes the following form:

The last ingredient we need are the transpositions $\sigma_{(i,j)}$ on $(V^{d})^{\otimes t}$ which act by interchanging the $i$th and the $j$th tensor factor. For example

with $x,y\in V^{d}$ arbitrary. Transpositions suffice, because they generate the full group of permutations. For $\left(V^{d}\right)^{\otimes 2}$ we only have

but for higher tensor systems more permutations can occur. Consequently, permutations act by interchanging different input and output lines and the wiring diagram representation allows one to keep track of this pictorially. In fact, only the input and output position of a line matters. We can use diagrams to simplify expressions by disentangling the corresponding lines. Take $\sigma_{(1,2)}$ on $\left(V^{d}\right)^{\otimes 2}$ as an example. Using wiring diagrams we can derive the standard result

pictorially. We are now ready to prove some important auxiliary results.

Let $A,B\in H^{d}$ be arbitrary. Then it holds that

which is, in our experience, a common misconception.

The basic formula (7) for $P_{\operatorname{Sym}^{2}}$ is given by

and the concepts from above allow us to translate this into the following wiring diagram:

(Note that this operator acts on the full tensor space $\left(V^{d}\right)^{\otimes 2}$, hence in the wiring diagram it is represented by a two-indexed box.) Applying the graphical calculus yields

Obviously, it is also possible to obtain (11) by direct calculation. We have included such a calculation in the appendix (Section 9.1) to demonstrate the complexity of direct calculations as compared to graphical ones.

We conclude this section with the following slightly more involved result.

Let $A,B,C\in H^{d}$ be arbitrary. Then it holds that

The proof can in principle be obtained by evaluating all permutations of 3 tensor systems algebraically and taking the partial trace afterwards. However, a pictorial calculation using wiring diagrams is much faster and more elegant.

For permutations of three elements, formula (7) implies

where. $\sigma_{2,1,3}(u\otimes v\otimes w)=(v\otimes u\otimes w)$, etc. This in turn allows us to write

Problem Setup

In the sampling process, we start by measuring the intensity of the signal:

which is a renormalized version of $\mathcal{A}^{*}\mathcal{A}:H^{d}\to H^{d}$. Concretely

The scaling is going to greatly simplify our analysis, because it guarantees that $\mathcal{R}$ is “near-isotropic”, as the following result shows.

The operator $\mathcal{R}$ defined in (16) is near-isotropic in the sense that

Let us start with deriving (17). For $Z\in H^{d}$ arbitrary we have

Here, (18) follows from the fact that the $a_{i}$’s are chosen iid from a $t$-design, (4.1) uses the fact that $\dim(\operatorname{Sym}^{2})=\binom{d+1}{2}^{-1}$ together with Definition 3, and the final line is an application of Lemma 6. ∎

Let now $x\in V^{d}$ be the signal we want to recover. As in we consider the space

(which is the tangent space of the manifold of all hermitian matrices at the point $X=xx^{*}$). This space is of crucial importance for our analysis. The orthogonal projection onto this space can be given explicitly:

We denote the projection onto its orthogonal complement with respect to the Frobenius inner product by $\mathcal{P}_{T}^{\perp}$. Then for any matrix $Z\in H^{d}$ the decomposition

is valid. We point out that in particular

holds. We will frequently use this fact. For a proof, consider $Z\in H^{d}$ arbitrary and insert the relevant definitions:

2. Convex Relaxation

Following the measurements (13) and (14) can be translated into matrix form by applying the following “lifts”:

By doing so the measurements assume the a linear form:

Hence, the phase retrivial problem becomes a matrix recovery problem. The solution to this is guaranteed to have rank 1 and encodes (up to a global phase) the unknown vector $x$ via $X=xx^{*}$. Relaxing the rank minimization problem (which would output the correct solution) to a trace norm minimization yields the now-familiar convex optimization problem

While this convex program is formally equivalent to the previously studied general-purpose matrix recovery algorithms , there are two important differences:

The measurement matrices $A_{i}$ are rank-1 projectors: $A_{i}=a_{i}a_{i}^{*}$.

The unknown signal is known to be proportional to a rank-1 projector ($X=xx^{*}$) as well.

3. Well-posedness / Injectivity

In this section, we follow to establish a certain injectivity property of the measurement operator $\mathcal{A}$. Compared to , our injectivity properties are somewhat weaker. Their proof used the independence of the components of the Gaussian measurement operator, which is not available in this setting, where individual vector components might be strongly correlated. We will pay the price for these weaker bounds in Section 6. There, we construct an “approximate dual certificate” that proves that the sought-for signal indeed minimizes the nuclear norm. Owing to the weaker bounds found here, the construction is more complicated than in . In the language of , we will have to carry out the full “golfing scheme”, as opposed to the “single leg” that proved sufficient in .

With probability of failure smaller than $d^{2}\exp(-\frac{3m}{384d})$ the inequality

is valid for all matrices $Z\in T$ simultaneously.

We aim to show the more general statement

Note that these summands have mean zero by construction. Furthermore observe that the auxiliary result (23) implies

follows. For the variance we use the standard identity

and focus on the last expression. Writing it out explicitly yields

where we have used the basic definition of $\mathcal{P}_{T}$ and $0\leq\operatorname{tr}(A_{i}X)=|\langle a_{i},x\rangle|^{2}\leq 1$. Consequently, for $Z\in T$ arbitrary

Here we have applied $\dim\operatorname{Sym}^{3}=\binom{d+2}{3}^{-1}$ and Lemma 7 in lines 3 and 4, respectively. Furthermore we used $Z\in T$ – hence $\mathcal{P}_{T}Z=Z$ and $\operatorname{tr}(Z)=\operatorname{tr}(XZ)$ – as well as the basic definition (22) of $\mathcal{P}_{T}$ to simplify the terms occuring in the fourth line. Putting everything together yields

and we can safely set $\sigma^{2}:=\frac{12d}{m}$. Now Theorem 5 tells us

for all $0\leq\delta\leq 1\leq 6d=\sigma^{2}/\underline{R}$. This gives the desired bound on the event

occuring. If this is not the case, (26) implies

for all matrices $Z\in T$ simultaneously. This is the general statement at the beginning of the proof and setting $\delta=1/2$ yields Proposition 9. ∎

Let $\mathcal{A}$ be as above with vectors sampled from a $t$-design ($t\geq 1$). Then the statement

holds with probability one for all matrices $Z\in H^{d}$ simultaneously.

where we have used $0\leq\Pi_{A_{i}}\leq\mathcal{I}$. ∎

Note that equation (27) can be improved. Indeed, a standard application of the Operator Bernstein inequality (Theorem 4) gives

for all matrices $Z\in T$ with probability of failure smaller than $d^{2}\exp\left(-Cm/d\right)$ for some $0<C\leq 1$. However, we actually do not require this tighter bound.

Proof of the Main Theorem / Convex Geometry

In this section, we will follow to prove that the convex program (24) indeed recovers the sought for signal $x$, provided that a certain geometric object – an approximate dual certificate – exists.

Consequently $X$ is guaranteed to be the unique minimum of (24), if

is true for every feasible $\Delta$. In order to show this we combine feasibility of $\Delta$ with inequalities (25) and (27) to obtain

Feasibility of $\Delta$ also implies $(Y,\Delta)=0$, because by defnition $Y$ is in the range of $\mathcal{A}^{*}$. Combining this insight with the defining property (28) of $Y$ and (30) yields

which is just the desired optimality criterion (29). ∎

Constructing the Dual Certificate

A straightforward approach to construct an approximate dual certificate would be to set

The key observation here is that the $t$-design property provides one with useful information about the first $t$ moments of the random variable $|\langle x,a_{i}\rangle|^{2}$. This knowledge allows us to explicitly bound the probability of “dangerously large overlaps” or “coherent measurement vectors” occurring.

Let $x\in V^{d}$ be an arbitrary vector of unit length. If $a$ is chosen uniformly at random from a $t$-design ($t\geq 1$) $D_{t}\subset V^{d}$, then the following is true for every $\gamma\leq 1$:

We aim to prove the slightly more general statement

which is valid for any $\delta\geq 1$. Setting $\delta=4$ then yields (32). The $t$-design property provides us with useful information about the first $t$ moments of the non-negative random variable $\xi=|\langle a,x\rangle|^{2}$. Indeed, with $A=aa^{*}$ it holds for every $k\leq t$ that

These inequalities are tight for the mean $\mu=\tau_{1}$ of $\xi$ and hence

Now we aim to use the well-known $t$-th moment bound

which is a straightforward generalization of Chebyshev’s inequality. Applying it, yields the desired result. Indeed,

The previous lemma bounds the probability of the undesired events

where $0\leq\gamma\leq 1$ is a fixed parameter which we refer to as the truncation rate. It turns out that a single truncation of this kind does not quite suffice yet for our purpose. We need to introduce a second truncation step.

Fix $Z\in T$ arbitrary and decompose it as

and define the two-fold truncated operator

where $1_{E_{i}}$ and $1_{G_{i}}$ denote the indicator functions associated with the events $E_{i}$ and $G_{i}$, respectively.

The following result shows that due to Lemma 13 this truncated operator is in expectation close to the original $\mathcal{R}$.

Fix $Z\in T$ arbitrary and let $\mathcal{R}_{Z}$ be as in (34). Then

We start by introducing the auxiliar (singly truncated) operator

Now use Lemma 13 to bound the first term:

and inserting these bounds into (36) yields the desired statement. ∎

We now establish a technical result which will allow us to find a suitable approximate dual certificate using the “golfing scheme” construction .

Fix $Z\in T$ arbitrary, let $\mathcal{R}_{Z}$ be as in (34). Assume that that the design order $t$ is at least 3 and the truncation rate $\gamma$ satisfies

Then for $1/4\leq b\leq 1$ and $c\geq\sqrt{2}b$ with probability at least $1-d\exp(-\frac{9mb}{640td^{2-\gamma}})$ one has

The statement is invariant under rescaling of $Z$. Therefore it suffices to treat the case $\|Z\|_{2}=1$. In this case we can decompose

which follows from $\gamma\leq 1-2/t$, $t\geq 3$ and $b\geq 1/4$. To obtain (38) we use a similar reasoning:

where we have used the fact that $\mathcal{P}_{T}$ projects onto a subspace of at most rank-2 matrices in the third line and (43) in the fourth. This motivates to define the event

which guarantees both (37) and (38) due to the assumption on $c$ and $\|Z\|_{2}=1$. So everything boils down to bounding the probability of $E^{c}$. We decompose

We will estimate this sum using the Operator Bernstein inequality (Theorem 4). Thus we need an a priori bound for the summands

as well as a bound for the variance. First observe that

With this ingredient we can now construct a suitable approximate dual certificate $Y$, closely following .

holds and the total number of measurements fulfills

The randomzied construction of $Y$ is summarized in Algorithm 1. If this algorithm succeeds, it outputs three lists

The recursive construction yields the following expressions (c.f. [70, Lemma 14]):

Then, in case of success, the validity of properties (37) and (38) for $c=1/2$ and $b=1/8$ in each step ($Q_{i}\to Q_{i+1}$ and $Y_{i}\to Y_{i+1}$, respectively) guarantee

Thus, $Y_{r}$ constitutes an approximate dual certificate in the sense of Def. 11.

Recall that the $\xi_{i}$’s are Bernoulli random variables which indicate whether the $i$-th iteration of the algorithm has been succesful ($\xi_{i}=1$), or failed ($\xi_{i}=0$). Our aim is to bound the probability of the event in (42) by a similar expression involving independent It was pointed out to us by A. Hansen that in some previous papers which involve a similar construction to the one presented here, it was tacitly assumed that the $\xi_{i}$ are independent. This will of course not be true in general. Fortunately, a more careful argument shows that all conclusions remain valid . Our treatment here is similar to the one presented in . Bernoulli variables $\xi_{i}^{\prime}$. To this end, write

is valid if $\xi_{l}^{\prime}$ is an independent $p^{\prime}$-Bernoulli distributed with

Proposition 16 provides a uniform lower bound on the success probability $p(\xi_{l-1},\dots,\xi_{1})$. Indeed, there is a universal constant $C_{1}$ such that invoking Prop. 16 with

and $Z=Q$ gives a probability of success of at least $9/10$ for any $Q$ (in particular, independently of the $\xi_{l-1},\dots,\xi_{1}$). Thus, choosing $p^{\prime}=9/10$ and $m_{i}=m$ for all $i$, we can then iterate the estimate (44) to arrive at

where the $\xi_{i}^{\prime}$ are independent Bernoulli variables with parameter $9/10$. A standard Chernoff bound (e.g. [72, Section Concentration: Theorem 2.1]) gives

Setting the number of iterations generously to

where we have used $\omega\geq 1\geq\log 2$ in the last inequality. Together with (42), (45) and (46) this gives the desired bound

on our construction of $Y$ failing. The total number of measurement vectors sampled is

Finally we are ready to put all pieces together and show or main result – Theorem 1.

Converse Bound

In this paper, we require designs of order at least three. Here we prove that this criterion is fundamental in the sense that sampling from 2-designs in general cannot guarantee a sub-quadtratic sampling rate. In order to do so, we will use a particular sort of 2-design, called a maximal set of mutually unbiased bases (MUBs) . Two orthonormal bases $\left\{u_{i}\right\}_{i=1}^{d}$ and $\left\{v_{i}\right\}_{i=1}^{d}$ are called mutually unbiased if their overlap is uniformly minimal. Concretely, this means that

must hold for all $i,j=1,\ldots,d$. Note that this is just a generalization of the incoherence property between standard and Fourier basis. In prime power dimensions, a maximal set of $(d+1)$ such MUBs is known to exist (and can be constructed) . Such a set is maximal in the sense that it is not possible to find more than $(d+1)$ MUBs in any Hilbert space. Among other interesting properties – cf. for a detailed survey – maximal sets of MUBs are known to form 2-designs .

The defining properties of a maximal set of MUBs allow us to derive the converse bound – Theorem 2.

Suppose that $m$ measurement vectors $a_{1},\dots,a_{m}$ are sampled independently and uniformly at random from $D_{2}$. Then, for any $\omega\geq 0$, the number of measurements must obey

Consequently a scaling of $\mathcal{O}(d^{2})$ in general cannot be avoided when demanding only the property of being a 2-design and simultaneously requiring a “reasonably small” probability of failure in the recovery process.

Suppose that $\left\{u_{i}\right\}_{i=1}^{d}$ is one orthonormal basis contained in the maximal set of MUBs $D_{2}$ and set $x:=u_{1}$ as well as $z:=u_{2}$. Note that by definition these vectors are orthogonal and normalized. Due to the particular structure of MUBs, $x$ and $z$ can only be distinguished if either $u_{1}$ or $u_{2}$ is contained in $\left\{a_{1},\ldots,a_{m}\right\}$. Since each $a_{i}$ is chosen iid at random from $D_{2}$ containing $(d+1)d$ elements, the probability of obtaining either $u_{1}$ or $u_{2}$ is $p=\frac{2}{(d+1)d}$. As a result, the problem reduces to the following standard stopping time problem (cf. for example Example (2) in Chapter 6.2 in ):

To answer this question, we have to find the smallest integer $m$ such that

for any $p\in[0,1/2]$ implies that (48) is a necessary criterion for (49) and we are done. ∎

Conclusion

In this paper we have derived a partly derandomized version of Gaussian PhaseLift . Instead of Gaussian random measurements, our method guarantees recovery for sampling iid from certain finite vector configurations, dubbed $t$-designs. The required sampling rate depends on the design order $t$:

For small $t$ this rate is worse than the Gaussian analogue – but still non-trivial. However, as soon as $t$ exceeds $2\log d$, we obtain linear scaling up to a polylogarithmic overhead.

In any case, we feel that the main purpose of this paper is not to present yet another efficient solution heuristics, but to show that the phase retrieval problem can be derandomized using $t$-designs. These finite vector sets lie in the vast intermediate region between random Fourier vectors and Gaussian random vectors (the Fourier basis is a $1$-design, whereas normalized Gaussian random vectors correspond to an $\infty$-design). Therefore the design order $t$ allows us to gradually transcend between these two extremal cases.

Acknowledgements DG and RK are grateful to the organizers and participants of the Workshop on Phaseless Reconstruction, held as part of the 2013 February Fourier Talks at the University of Maryland, where they were introduced to the details of the problem. This extends, in particular, to Thomas Strohmer.

The work of DG and RK is supported by the Excellence Initiative of the German Federal and State Governments (Grant ZUK 43), by scholarship funds from the State Graduate Funding Program of Baden-Württemberg, and by the US Army Research Office under contracts W911NF-14-1-0098 and W911NF-14-1-0133 (Quantum Characterization, Verification, and Validation), and the DFG. FK acknowledges support from the German Federal Ministry of Education and Reseach (BMBF) through the cooperative research project ZeMat.

References

Appendix

Here we briefly state an elementary proof of Lemma 6. In the main text we proved this result using wiring diagrams. The purpose of this is to underline the relative simplicity of wiring diagram calculations. Indeed, the elementary proof below is considerably more cumbersome than its pictorial counterpart.

Let us choose an arbitrary orthonormal basis $b_{1},\ldots,b_{d}$ of $V^{d}$. In the induced basis $\left\{b_{i}\otimes b_{j}\right\}_{i,j=1}^{d}$ of $V^{d}\otimes V^{d}$ the transpositions then correspond to

This choice of basis furthermore allows us to write down $\operatorname{tr}_{2}(A)$ for $A\in M^{d}\otimes M^{d}$ explicity:

Consequently we get for $A,B\in H^{d}$ arbitrary

The latter term can be evaluated explicitly: