A Generalization of Wirtinger Flow for Exact Interferometric Inversion

TL;DR

This paper extends Wirtinger Flow to interferometric inversion, enabling exact signal recovery from cross-correlations under broader measurement models with theoretical guarantees and practical validation.

Contribution

It introduces Generalized Wirtinger Flow (GWF), unifying low rank matrix recovery and non-convex optimization, and provides less restrictive conditions for exact recovery.

Findings

GWF covers a larger class of measurement models.

The sufficient condition for exact recovery is less stringent than LRMR.

Numerical experiments validate GWF in radar imaging scenarios.

Abstract

Interferometric inversion involves recovery of a signal from cross-correlations of its linear transformations. A close relative of interferometric inversion is the generalized phase retrieval problem, which consists of recovering a signal from the auto-correlations of its linear transformations. Recently, significant advancements have been made in phase retrieval methods despite the ill-posed, and non-convex nature of the problem. One such method is Wirtinger Flow (WF), a non-convex optimization framework that provides high probability guarantees of exact recovery under certain measurement models, such as coded diffraction patterns, and Gaussian sampling vectors. In this paper, we develop a generalization of WF for interferometric inversion, which we refer to as Generalized Wirtinger Flow (GWF). GWF theory extends the probabilistic exact recovery results of WF to arbitrary measurement models characterized in the equivalent lifted problem, hence covers a larger class of measurement models. Our framework unifies the theory of low rank matrix recovery (LRMR) and the non-convex optimization approach of WF, thereby establishes theoretical advantages of the non-convex approach over LRMR. We identify a new sufficient condition on the lifted forward model that directly implies exact recovery conditions of standard WF. This condition is less stringent than those of LRMR, which is the state of the art approach for exact interferometric inversion. We establish our sufficient condition for the cross-correlations of linear measurements collected by complex Gaussian sampling vectors, and show that the exact recovery conditions of standard WF imply our sufficient condition. As a result, the regularity condition of WF becomes redundant in solving the interferometric inversion problem. Finally, we demonstrate the effectiveness of GWF numerically in a deterministic multi-static radar imaging scenario.

Equations

$$f_i^m=⟨{\mathbf{L}}_i^m,{\mathbf{\rho }}_t⟩,f_j^m=⟨{\mathbf{L}}_j^m,{\mathbf{\rho }}_t⟩,m=1,\cdots ,M,$$

$$d_{ij}^m=f_i^m\overline{f_j^m}=({\mathbf{L}}_i^m{)}^H{\mathbf{\rho }}_t{\mathbf{\rho }}_t^H{\mathbf{L}}_j^{m}m=1,\cdots M,$$

$$\mathcal{J}(\mathbf{\rho }):=\frac{1}{2M}\sum _{m=1}^M|{\left({\mathbf{L}}_i^m\right)}^H\mathbf{\rho }{\mathbf{\rho }}^H{\mathbf{L}}_j^m-d_{ij}^m{|}^2,$$

$$\hat{\mathbf{\rho }}=\underset{\mathbf{\rho }}{\mathrm{argmin}}\mathcal{J}(\mathbf{\rho }).$$

$$\mathcal{L}(\mathbf{\rho })=\mathbf{y},\text{where}\mathbf{\rho }\in {\mathbb{C}}^N,y_{ij}^m={\left({\mathbf{L}}_i^m\right)}^H\mathbf{\rho }{\mathbf{\rho }}^H{\mathbf{L}}_j^m.$$

$$P:=\{e^{\mathrm{i}\varphi }{\mathbf{\rho }}_t:\varphi \in [0,2\pi )\},$$

$$P_{\mathbf{\rho }}:=\{e^{\mathrm{i}\varphi }\mathbf{\rho },\varphi \in [0,2\pi )\},$$

$$d_{ij}^m=⟨{\mathbf{F}}^m,{\tilde{\mathbf{\rho }}}_t{⟩}_{F}m=1,...,M$$

$$\mathbf{d}=\mathcal{F}({\tilde{\mathbf{\rho }}}_t)$$

$$\text{find:}\mathbf{X}\text{s.t. }\mathbf{d}=\mathcal{F}(\mathbf{X}),\mathbf{X}⪰0,\text{rank}(\mathbf{X})=1.$$

$$\text{minimize:}\text{rank}(\mathbf{X})\text{s.t. }\mathbf{d}=\mathcal{F}(\mathbf{X}),\mathbf{X}⪰0,$$

$$\text{minimize:}\text{tr}(\mathbf{X})\text{s.t. }\mathbf{d}=\mathcal{F}(\mathbf{X}),\mathbf{X}⪰0,$$

$$\text{minimize:}\text{tr}(\mathbf{X}),\text{s.t. }\frac{1}{2}\Vert \mathcal{F}(\mathbf{X})-\mathbf{d}{\Vert }_2^2\le \mathcal{O}({\sigma }^2),\mathbf{X}⪰0,$$

$$\underset{\mathbf{X}⪰0}{\text{minimize:}}\frac{1}{2}\Vert \mathcal{F}(\mathbf{X})-\mathbf{d}{\Vert }_2^2+\lambda \text{tr}(\mathbf{X}).$$

$${\mathbf{X}}_k$$

$${\mathbf{\nu }}_k$$

$${\mathbf{\rho }}_{k+1}={\mathbf{\rho }}_k-\frac{{\mu }_{k+1}}{\Vert {\mathbf{\rho }}_0{\Vert }^2}\nabla \mathcal{J}({\mathbf{\rho }}_k).$$

$$\nabla \mathcal{J}$$

$$\frac{\partial }{\partial \mathbf{\rho }}$$

$$\nabla \mathcal{J}({\mathbf{\rho }}_k)=\frac{1}{2M}\sum _{m=1}^M\left[\overline{e_{ij}^m}({\mathbf{L}}_j^m{\left({\mathbf{L}}_i^m\right)}^H{\mathbf{\rho }}_k)+e_{ij}^m({\mathbf{L}}_i^m{\left({\mathbf{L}}_j^m\right)}^H{\mathbf{\rho }}_k)\right],$$

$$\nabla \mathcal{J}({\mathbf{\rho }}_k)=\frac{1}{M}\left[\frac{1}{2}\left(\sum _{m=1}^M\overline{e_{ij}^m}({\mathbf{L}}_j^m{\left({\mathbf{L}}_i^m\right)}^H)+\sum _{m=1}^M{e}_{ij}^m({\mathbf{L}}_i^m{\left({\mathbf{L}}_j^m\right)}^H)\right)\right]{\mathbf{\rho }}_k.$$

$${\mathcal{F}}^H(\mathbf{y})=\sum _{m=1}^M{y}^m({\mathbf{L}}_i^m{\left({\mathbf{L}}_j^m\right)}^H),$$

$$\nabla \mathcal{J}({\mathbf{\rho }}_k)=\frac{1}{M}{\mathcal{P}}_S\left({\mathcal{F}}^H(\mathbf{e})\right){\mathbf{\rho }}_k,$$

$$\text{minimize: }\frac{1}{2M}\Vert \mathcal{F}(\mathbf{X})-\mathbf{d}{\Vert }_2^2\text{s.t. }\text{rank}(\mathbf{X})=1\text{and}\mathbf{X}⪰0.$$

$$\text{minimize: }\frac{1}{2M}\Vert \mathcal{F}(\mathbf{X})-\mathbf{d}{\Vert }_2^2\text{s.t.}\mathbf{X}=\mathbf{\rho }{\mathbf{\rho }}^H.$$

$${\mathbf{\rho }}_{k+1}={\mathbf{\rho }}_k-{\mu }_{k+1}{\left(\frac{\partial \overline{\mathbf{X}}}{\partial \bar{\mathbf{\rho }}}(\frac{\partial \mathcal{J}}{\partial \overline{\mathbf{X}}}+\frac{\partial \mathcal{J}}{\partial \mathbf{X}})\right)}^T,\text{where}$$

$${\left(\frac{\partial \overline{\mathbf{X}}}{\partial \overline{\mathbf{\rho }}}\frac{\partial \mathcal{J}}{\partial \overline{\mathbf{X}}}\right)}^T=\frac{1}{2M}({\mathcal{F}}^H\mathcal{F}(\mathbf{X})-{\mathcal{F}}^H\mathbf{d})\mathbf{\rho },{\left(\frac{\partial \overline{\mathbf{X}}}{\partial \overline{\mathbf{\rho }}}\frac{\partial \mathcal{J}}{\partial \mathbf{X}}\right)}^T=\frac{1}{2M}({\mathcal{F}}^H\mathcal{F}(\mathbf{X})-{\mathcal{F}}^H\mathbf{d}{)}^H\mathbf{\rho }.$$

$$\mathrm{dist}(\mathbf{\rho },{\mathbf{\rho }}_t)=\Vert \mathbf{\rho }-e^{\mathrm{i}\Phi (\mathbf{\rho })}{\mathbf{\rho }}_t\Vert ,\text{where}\Phi (\mathbf{\rho }):=\underset{\varphi \in [0,2\pi ]}{\mathrm{argmin}}\Vert \mathbf{\rho }-{\mathbf{\rho }}_t{e}^{\mathrm{i}\varphi }\Vert .$$

$${\mathrm{dist}}^2(\mathbf{\rho },{\mathbf{\rho }}_t)=\Vert \mathbf{\rho }{\Vert }^2+\Vert {\mathbf{\rho }}_t{\Vert }^2-2|⟨\mathbf{\rho },{\mathbf{\rho }}_t⟩|,$$

$$\cos \theta =\frac{|⟨\mathbf{\rho },{\mathbf{\rho }}_t⟩|}{\Vert \mathbf{\rho }\Vert \Vert {\mathbf{\rho }}_t\Vert },$$

$$\Vert \mathbf{\rho }{\mathbf{\rho }}^H-{\mathbf{\rho }}_t{\mathbf{\rho }}_t^H{\Vert }_F^2=\Vert \mathbf{\rho }{\mathbf{\rho }}^H{\Vert }_F^2+\Vert {\mathbf{\rho }}_t{\mathbf{\rho }}_t^H{\Vert }_F^2-2\mathrm{Re}(⟨\mathbf{\rho }{\mathbf{\rho }}^H,{\mathbf{\rho }}_t{\mathbf{\rho }}_t^H{⟩}_F).$$

$${\mathcal{P}}_{\mathcal{X}}={\mathcal{P}}_{r=1}\circ {\mathcal{P}}_{+},$$

Full text

Key words— interferometric inversion, Wirtinger Flow, phase retrieval, wave-based imaging, interferometric imaging, low rank matrix recovery, PhaseLift

1 Introduction

Interferometric inversion involves the recovery of a signal of interest from the cross-correlations of its linear measurements, each collected by a different sensing process. Let $\mathbf{L}^{m}_{i},\mathbf{L}^{m}_{j}\in\mathbb{C}^{N}$ denote the $m^{th}$ sampling vectors of the $i^{th}$ and $j^{th}$ sensing processes and $\bm{\rho}_{t}\in\mathbb{C}^{N}$ be the ground truth/signal of interest. We define

[TABLE]

as the linear measurements and describe the cross-correlated measurements as

[TABLE]

where $\overline{(\cdot)}$ denotes complex conjugation. Thus, interferometric inversion involves recovery of $\bm{\rho}_{t}\in\mathbb{C}^{N}$ from $d^{m}_{ij}\in\mathbb{C},\ m=1,...,M$ using the model in (2).

Interferometric inversion problem arises in many applications in different disciplines. These include radar and sonar interferometry [2, 3, 4], passive imaging in acoustic, electromagnetic and geophysical applications [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], interferometric microscopy [21], beamforming and sensor localization in large area networks [22], [23, 24, 25] among others. Additionally, cross-correlations were shown to provide robustness to statistical fluctuations in scattering media or incoherent sources in wave-based imaging [26, 27], and with respect to phase errors in the correlated linear transformations [28, 29, 30, 31]. Therefore, in applications such as passive imaging [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20] and interferometry [2, 3, 4], cross-correlations are formed as a part of the inference process after acquiring linear measurements by sensors that are configured differently in space, time or frequency. Additionally, the cross-correlated measurement model arises naturally from the underlying physical sensing processes in certain applications such as optical and radio astronomy [32, 33], or quantum optical imaging [34].

A special case of the interferometric inversion problem is when $i=j$ in (2), in which the model becomes the auto-correlations of linear measurements collected by a single sensing process. In this case, the interferometric inversion problem reduces to the well-known phase retrieval problem. Notably, both problems are non-convex due to the quadratic equality constraints enforced by the correlated measurement model. In recent years, several phase retrieval methods with exact recovery guarantees have been developed despite the non-convex nature of the problem. These methods are characterized by either one or both of the following two principles: convexification of the solution set, which includes lifting based approaches [35, 36, 37], or a provably accurate initialization, followed by an algorithmic map that refines the initial estimate which is most prominently established by Wirtinger Flow (WF) [1] and its variants [38, 39, 40, 41, 42, 43, 44, 45].

In methods that deploy lifting, such as PhaseLift [35, 36], signal recovery from quadratic measurements is reformulated as a low rank matrix recovery (LRMR) problem. While the LRMR approach offers convergence guarantees via convexification, it has limited practical applicability in typical sensing and imaging problems since lifting increases the dimension of the inverse problem by an order of magnitude and requires demanding memory storage in implementation. The WF framework, on the other hand, avoids lifting hence offers considerable advantages in computational complexity and memory requirements. Despite solving the non-convex problem directly, convergence to a global solution at a geometric rate is guaranteed by WF for coded diffraction patterns and Gaussian sampling vectors [1], and more recently for short time Fourier transforms [44]. These advantages promote WF as a theoretical framework suitable for large scale imaging problems.

Conventionally, interferometric inversion in imaging applications has been approached by Fourier based techniques, such as time or frequency difference of arrival (TDOA/FDOA) backprojection [5, 9, 6, 16, 13, 46, 18, 19]. While these methods are practical and computationally efficient, their applicability is limited to scenes composed of well-separated point targets due to underlying assumptions. As an alternative, LRMR theory has been explored for interferometric inversion [11, 47]. In particular, [47] assesses the robustness of solutions to several convex semi-definite programs (SDPs) in a deterministic framework after lifting the interferometric measurements. Notably, these SDP solvers are inspired by the PhaseLift method [35, 36, 48], hence suffer from the same drawbacks in practice. In [11], an iterative optimization approach to LRMR was developed for interferometric passive imaging to circumvent the poor scaling properties of SDP approaches. The methodology in [11] is based on Uzawa’s method for matrix completion [49, 50]. While this method is computationally more efficient than the SDP solvers, it still operates on the lifted domain, hence requires significant memory and computational resources. Additionally, these convexified lifting based solvers require stringent theoretical conditions on the measurement model, which poses a major theoretical barrier for interferometric inversion problems with deterministic forward models.

In this paper, motivated by its advantages over lifting based methods, we develop a generalization of WF applicable to interferometric inversion problems with exact recovery guarantees. We refer to this method as the Generalized Wirtinger Flow (GWF). Beyond the immediate extension of algorithmic principles of WF to the interferometric inversion, our mathematical framework differs from that of WF in two significant ways. i) The GWF theory is established in the lifted domain. This lifting-based perspective allows us to develop a novel approach that bridges the theory between LRMR and the non-convex methodology of WF, and culminates into derivation of a new sufficient condition for exact recovery in the context of interferometric inversion. Namely, GWF guarantees exact recovery to a general class of problems that are characterized over the equivalent lifted domain by the restricted isometry property (RIP) on the set of rank-1, positive semi-definite (PSD) matrices. Specifically, we provide an upper bound for the restricted isometry constant (RIC) over the set of rank-1, PSD matrices so that RIP directly implies convergence to a global solution. ii) Given the sufficient condition, we derive exact recovery guarantees for an arbitrary mapping over the lifted domain using solely geometric arguments in a deterministic framework, rather than model specific probabilistic arguments of standard WF. Notably, our analysis proves that the regularity condition, which forms the basis of convergence in standard WF, is redundant for the case of interferometric inversion.

In addition, developing the GWF framework through the equivalent lifted problem allows us to identify the key theoretical advantages of the non-convex optimization approach over lifting based convex solvers beyond the immediate gains in computation and applicability. Our exact recovery condition for GWF proves to be a less stringent restricted isometry property than that of the sufficient conditions of LRMR methods. As a result of its fundamental differences from standard WF, and less stringent exact recovery requirements than those of LRMR methods, GWF offers a promising step towards exact recovery guarantees for interferometric inversion problems which are governed by deterministic measurement models. One such application is array imaging, for which RIP was shown to hold over rank-1 matrices for sufficiently high central frequencies, albeit asymptotically in the number of receivers [51]. Accordingly, we consider GWF as an exact interferometric imaging method alternative to LRMR for multi-static radar imaging, and design imaging system parameters for the lifted forward model to satisfy of our sufficient condition in the non-asymptotic regime [52].

In our key results, we first prove that the RIP over rank-1, PSD matrices on the lifted forward model implies an accurate initialization by the spectral method, and ensures that the regularity condition is satisfied in the $\epsilon$-neighborhood defined by the initialization if RIC is less than $0.214$. Next, although not applicable to the case of auto-correlations in [1], we show that this sufficient condition is satisfied for the cross-correlation of linear measurements collected by i.i.d. Gaussian sampling vectors. Essentially, our results establish that the probabilistic arguments used for the Gaussian sampling model in [1] are fully captured by the single event that the RIP is satisfied over rank-1, PSD matrices in the case $i\neq j$ in (2). To validate our analysis we conduct numerical experiments for the Gaussian sampling model by counting empirical probability of exact recovery. We then demonstrate the effectiveness of GWF in a realistic passive multi-static radar imaging scenario. Our preliminary numerical simulations confirm our theoretical results and show that the interferometric inversion is solved in an exact manner by GWF.

The rest of our paper is organized as follows. We first introduce the problem formulation for interferometric inversion on the signal domain, and discuss lifting based prior art in Section 2. We then formulate the GWF algorithm in Section 3, present key definitions and terminology followed by the main theorem statements in Section 4. The proofs of the main theorems are provided in Section 5. The numerical simulations for Gaussian sampling model and interferometric multi-static radar imaging are presented in Section 6. Section 7 concludes our paper. Appendices A, B, and C include proofs of lemmas used in Section 4.

2 Problem Formulation and Prior Art

In this section, we introduce the non-convex formulation of the interferometric inversion problem and its key challenges. Next we discuss lifting based, convex formulations which address these key challenges in solving quadratic systems of equations via low rank matrix recovery theory.

2.1 The Non-Convex Objective Function

To address the interferometric inversion problem, we define the following objective function, and set up the corresponding optimization problem:

[TABLE]

Let $\mathbf{y}=[y_{ij}^{1},y_{ij}^{2},\cdots y_{ij}^{M}]^{T}$ and $\mathcal{L}:\mathbb{C}^{N}\rightarrow\mathbb{C}^{M}$ be the cross-correlated measurement map defined as

[TABLE]

The objective function in (3) is the $\ell_{2}$ mismatch in the range of $\mathcal{L}$, i.e., the space of cross-correlated measurements, which is solved over the signal domain $\mathbb{C}^{N}$ in (4). The difficulty in (4) is that the objective function $\mathcal{J}$ is non-convex over the variable $\bm{\rho}$ due to the invariance of cross-correlated measurements to global phase factors. Essentially (4) has a non-convex solution set with infinitely many elements, which casts interferometric inversion a challenging, ill-posed problem.

Definition 2.1** (Global Solution Set).**

We say that the points

[TABLE]

form the global solution set for the interferometric inversion from the cross-correlated measurements (2).

More generally, for any $\bm{\rho}\in\mathbb{C}^{N}$, let $\mathcal{E}_{\bm{\rho}}=\{\mathbf{z}\in\mathbb{C}^{N}:\mathcal{L}(\mathbf{z})=\mathcal{L}(\bm{\rho})\}$ be the equivalence class of $\bm{\rho}$ under $\mathcal{L}$. We then define the following collection of signals as the equivalence set of $\bm{\rho}$.

Definition 2.2** (Equivalence Set).**

Let $\bm{\rho}\in\mathbb{C}^{N}$ and

[TABLE]

We refer to $\mathit{P}_{\bm{\rho}}$ as the equivalence set of $\bm{\rho}$.

Remark**.**

Note that $\mathit{P}_{\bm{\rho}}\subset\mathcal{E}_{\bm{\rho}}$; and $\mathit{P}_{\bm{\rho}_{t}}$ is identical to the global solution set in Definition (2.1).

Alleviating the non-injectivity of the measurement map is a key step in formulating methods that guarantee exact recovery in phase retrieval literature [53], and offers us a blueprint in addressing (4). A key observation is that one can consider (5) as a mapping from a rank-1, positive semi-definite matrix $\bm{\rho}\bm{\rho}^{H}\in\mathbb{C}^{N\times N}$ instead of a quadratic map from the signal domain in $\mathbb{C}^{N}$, and attempt to recover $\bm{\rho}_{t}\bm{\rho}_{t}^{H}$. This approach is known as the lifting technique which is the main premise of LRMR based phase retrieval [35, 36, 37, 54] and interferometric inversion methods [11, 47, 48].

2.2 Low Rank Matrix Recovery via Lifting

We adopt the concepts of the LRMR approach to the interferometric inversion problem based on PhaseLift [11], [36], and introduce the following definitions.

Definition 2.3**.**

Lifting. Each correlated measurement in (2) can be written in the form of an inner product of two rank-1 operators, $\tilde{\bm{\rho}}_{t}=\bm{\rho}_{t}\bm{\rho}_{t}^{H}$ and $\overline{\mathbf{F}}^{m}=\mathbf{L}_{j}^{m}(\mathbf{L}_{i}^{m})^{H}$ such that111$(\bar{\cdot})$ denotes element-wise complex conjugation

[TABLE]

where $\langle\cdot,\cdot\rangle_{F}$ is the Frobenius inner product. We refer to the procedure of transforming interferometric inversion over $\mathbb{C}^{N}$ to the recovery of the rank-1 unknown $\tilde{\bm{\rho}}_{t}$ in $\mathbb{C}^{N\times N}$ as lifting.

The lifting technique introduces a new linear measurement map which we define as follows:

Definition 2.4**.**

Lifted Forward Model. Let $\mathbf{d}=[d_{ij}^{1},d_{ij}^{2},\cdots d_{ij}^{M}]\in\mathbb{C}^{M}$ denote the vector obtained by stacking the cross-correlated measurements in (2). Then using (7), we define $\mathcal{F}:\mathbb{C}^{N\times N}\rightarrow\mathbb{C}^{M}$ as follows

[TABLE]

and refer to $\mathcal{F}$ as the lifted forward model/map.

Remark**.**

The map $\mathcal{F}$ can be interpreted as a $M\times N^{2}$ matrix with $\bar{\mathbf{F}}^{m}$ as its rows in the vectorized problem in which $N\times N$ variable $\tilde{\bm{\rho}}_{t}$ is concatenated into a vector.

We refer to the problem of recovering of $\tilde{\bm{\rho}}_{t}$ from $\mathbf{d}$ using the model (8) as the lifted problem, or interferometric inversion in the lifted domain.

The main advantage of lifting is that for all $\bm{\rho}\in\mathbb{C}^{N}$, each non-convex equivalence set $\mathit{P}_{\bm{\rho}}$ is now mapped to a set with a single element $\tilde{\bm{\rho}}=\bm{\rho}\bm{\rho}^{H}$. Using the definition of the lifted forward model, quadratic equality constraints reduce to affine equality constraints to define a convex manifold in $\mathbb{C}^{N\times N}$. Using the rank-1 positive semi-definite (PSD) structure of the unknown $\tilde{\bm{\rho}}_{t}$, interferometric inversion in the lifted domain can be formulated as the following optimization problem:

[TABLE]

Here, we refer to $\mathbf{X}\in\mathbb{C}^{N\times N}$ as the lifted variable. Due to the fact that there surely exists a rank-1 solution from (8), (9) is equivalent to:

[TABLE]

which is known to be an NP-hard problem [49, 50]. Given its rank-minimization form, (10) is approached by low rank matrix recovery theory analogous to compressive sensing [55, 56]. Most prominently, the non-convex rank term of the objective function is relaxed by a convex surrogate, which under the positive semi-definite constraint, corresponds to the trace norm. This results in the following formulation [51, 35, 36]:

[TABLE]

which can be solved in polynomial time via semi-definite programming. To obtain a robust program, (11) is commonly perturbed using the least squares criterion, which under the additive i.i.d. noise assumption, can be written as

[TABLE]

where $\sigma^{2}$ corresponds to the variance in the case of Gaussian noise, and the order in the threshold can be tuned to obtain a desired lower bound on the log-likelihood of observing $\mathbf{d}$. Alternatively, one can equivalently formulate the Lagrangian of (12) with a proper choice of the regularization parameter $\lambda$ in relation to the threshold, as follows:

[TABLE]

(13) can be solved by Uzawa’s method [49, 50] which is analogous to the singular value thresholding algorithm with a PSD constraint [11] with the following iterations:

[TABLE]

In (14) and (15), $\bm{\nu}=[\nu^{1},\cdots\nu^{M}]^{T}$ denote the Lagrange multipliers initialized as $\bm{\nu}_{0}=0$, $\mu_{k}$ is the step size, subscript $k$ denotes the iteration number. $\mathcal{P}_{\tau}$ is the shrinkage operator acting on the singular values of its argument with threshold $\tau_{k}=\mu_{k}\lambda$, which through the parameter $\lambda$ of trace regularization, enforces the low-rank constraint. $\mathcal{P}_{+}$ is the projection operator onto the PSD cone.

As $M\ll N^{2}$ in typical estimation problems, the lifted forward model has a non-trivial null space. LRMR theory encapsulates identifying necessary and sufficient conditions on $\mathcal{F}$ in order to guarantee exact recovery despite having an under-determined system of equations in (8). The key conditions on $\mathcal{F}$ are primarily characterized by its null space [57, 58, 59] or restricted isometry properties [50, 60, 61] on low rank matrices. Methods such as PhaseLift [35, 36] and PhaseCut [37] assert conditions on the mapping $\mathcal{F}$ such that there exists no feasible element in the PSD cone with a smaller trace than the true solution $\tilde{\bm{\rho}_{t}}$, from which exact recovery results to a unique minimizer are directly implied by the standard arguments of semi-definite programming [36]. Iterative optimization via Uzawa’s method on the other hand requires RIP over rank-5 matrices with a sufficiently small RIC ($\leq{1}/{10}$) for the convex problem in (13) to have the identical solution to the original non-convex problem in (10). Notably, in this scenario the PSD constraint can simply be dropped due to the guaranteed uniqueness of the solution to the original problem, and the convexity of the problem in (13). Furthermore, under the properties on $\mathcal{F}$ asserted by PhaseLift, it is observed in [48] that the lifted problem can be robustly solved as a convex feasibility problem by Douglas-Rachford splitting by eliminating the trace minimization step completely. These indicate a level of redundancy between trace regularization and the PSD constraint, which arises from the convexification of the problem.

Altogether, lifting based approaches provide a profound perspective to the interferometric inversion problem. Our observation is that the key principles of lifting based methods in establishing exact recovery guarantees are reciprocated in the non-convex framework of WF. In fact, WF corresponds to solving a perturbed non-convex feasibility problem over the lifted domain, and in this sense, is reminiscent of the optimizationless PhaseLift method of [48], and Uzawa’s iterations in [11]. To observe this, we introduce the GWF iterations for interferometric inversion, and develop the method as a solver in the lifted problem framework. The basis of our extension from WF to GWF framework is the identification of conditions on the lifted forward model for the exactness of a solution in the lifted domain.

It is worth noting that low rank models and LRMR theory cover a rich area of research with wide range of applications in fields such as control theory and machine learning. Several of the LRMR methods discussed in this subsection were first proposed for matrix completion, and relate strongly to the problem of low rank matrix factorization. Non-convex optimization theory for low rank matrix factorization has undergone notable developments in recent years, with methods that offer exact recovery results for quadratic or bilinear equations of rank-$r$ matrices in the random Gaussian model [62], or if the measurement map satisfies RIP over rank-$6r$ with a RIC $\leq 1/10$ [63]. For further discussion on advances in non-convex low rank matrix recovery, we refer the reader to [64, 65].

3 Generalized Wirtinger Flow Framework for Exact Interferometric Inversion

In this section, we present the algorithmic principles of the GWF framework for exact interferometric inversion. We specifically identify and present the theoretical advantages obtained from viewing the algorithm in the lifted domain, despite operating on the signal domain. We first present the GWF iterations, then proceed with generalizing the spectral method for initialization using our lifted formulation of GWF. We finally provide an algorithm summary with specifications of computational complexity of each step.

3.1 GWF Iterations

In presenting the GWF iterations, we begin by extending the iterative scheme of WF to the case of interferometric measurement model in (2), with $i\neq j$. We next introduce an equivalent, novel interpretation of these iterations in the lifted domain, which yields the basis of our GWF framework for exact recovery with arbitrary lifted forward models.

3.1.1 Extending WF

In contrast to the lifting based approaches, the non-convex form of the problem in (4) is preserved in WF. Given an accurate initial estimate $\bm{\rho}_{0}$, WF involves using the following updates to refine the current estimate $\bm{\rho}_{k}$:

[TABLE]

Notably, $\mathcal{J}$ is a real-valued function of a complex variable $\bm{\rho}\in\mathbb{C}^{N}$, and therefore, non-holomorphic. Hence, the gradient over $\mathcal{J}$ is defined by the means of Wirtinger derivatives

[TABLE]

and $\bm{\rho}=\bm{\rho}_{R}+\mathrm{i}\bm{\rho}_{I}$, with $\bm{\rho}_{R},\bm{\rho}_{I}\in\mathbb{R}^{N}$. Thus, the iterations in (16) correspond to that of the steepest descent method [1], where $\mu_{k+1}$ is the step size. For interferometric inversion by solving (4), $\nabla\mathcal{J}$ evaluated at $\bm{\rho}_{k}$ is given by

[TABLE]

where $e_{ij}^{m}=(\left(\mathbf{L}_{i}^{m}\right)^{H}\bm{\rho}_{k}\bm{\rho}_{k}^{H}\mathbf{L}_{j}^{m}-d_{ij}^{m})$ is the mismatch between the synthesized and cross-correlated measurements. Note that in the case of phase retrieval we have $i=j$. Letting $\mathbf{L}^{m}=\mathbf{L}_{i}^{m}=\mathbf{L}_{j}^{m}$ and $d^{m}=|\langle\mathbf{L}^{m},\bm{\rho}_{t}\rangle|^{2}$, for $m=1\cdots M$, (19) reduces to the standard WF iterations of [1]. See Appendix A.1 for the derivation of $\nabla\mathcal{J}$ for $i\neq j$.

3.1.2 Interpretation of GWF Updates in the Lifted Domain

Our formulation of GWF framework resides in the lifted domain, and reveals an illuminating interpretation of the updates in (19). Moving the common term of the current iterate $\bm{\rho}_{k}$ outside the summation, (19) can be expressed as

[TABLE]

From the definition of the lifted forward model in (8), the second term inside the brackets in (20) becomes the backprojection of the measurement error onto the adjoint space of $\mathcal{F}$.

Definition 3.1**.**

Backprojection. Let $\mathbf{y}=[y^{1},y^{2},\cdots y^{M}]\in\mathbb{C}^{M}$. For the lifted forward model $\mathcal{F}:\mathbb{C}^{N\times N}\rightarrow\mathbb{C}^{M}$ in Definition 8, we define the adjoint operator $\mathcal{F}^{H}:\mathbb{C}^{M}\rightarrow\mathbb{C}^{N\times N}$ as,

[TABLE]

and refer to (21) as the backprojection of $\mathbf{y}$.

Since (20) consists of the average of terms that are Hermitian transposes of each other, the update term in (16) corresponds to backprojecting the mismatch between synthesized and true measurements, and then projecting it onto the set of symmetric matrices $S$, i.e.,

[TABLE]

where $\mathbf{e}=[e_{ij}^{1},e_{ij}^{2},\cdots e_{ij}^{M}]\in\mathbb{C}^{M}$ is the measurement mismatch vector and $\mathcal{P}_{S}(\cdot)$ denotes the projection operator onto $S$.

The representation in (22) provides a novel perspective in interpreting GWF as a solver of the lifted problem. Consider the structure enforced by (16) for the non-relaxed form of the low rank recovery problem in (10). In the GWF updates, the rank of the lifted variable $\mathbf{X}$ in (10) is merely fixed at one, and the rank minimization problem is converted to its original non-convex feasibility problem in (9). Knowing that (3) corresponds to the $\ell_{2}$ mismatch in the space of cross-correlated measurements, the GWF solver of (3) is equivalently the solver of the following lifted problem:

[TABLE]

Observe that the rank-1, PSD constraint precisely corresponds to minimization over the set of elements of $\mathcal{X}=\{\bm{\rho}\bm{\rho}^{H},\bm{\rho}\in\mathbb{C}^{N}\}$ and that (23) can be equivalently cast as:

[TABLE]

The updates to solve (24) then can be performed on the leading eigenspace of the lifted variable $\mathbf{X}$ directly by means of the Jacobian $\frac{\partial\bar{\mathbf{X}}}{\partial\bar{\bm{\rho}}}$. Since $\mathbf{X}=\mathbf{X}^{H}$, this yields222We use the property of Wirtinger derivatives in writing the update (25), such that the derivative of a real valued function of a complex variable has the property that $\overline{(\frac{\partial\mathcal{J}}{\partial\mathbf{X}})}={(\frac{\partial\mathcal{J}}{\partial\overline{\mathbf{X}}})}$

[TABLE]

[TABLE]

Initializing the algorithm with a proper estimate in the constraint set $\mathcal{X}$, i.e. $\mathbf{X}_{0}=\bm{\rho}_{0}\bm{\rho}_{0}^{H}$, and substituting $\mathbf{X}=\bm{\rho}\bm{\rho}^{H}$, we precisely obtain $\nabla\mathcal{J}(\bm{\rho}_{k})$ derived in (22) when (25) is evaluated at $\bm{\rho}=\bm{\rho}_{k}$.

3.2 The Distance Metric and Equivalence of Convergence in Signal and Lifted Domains

While the view of GWF in the lifted problem and the formulation in (24) are illuminating, the algorithmic map of GWF operates exclusively on the signal domain in $\mathbb{C}^{N}$. The duality between the lifted domain and the signal domain is established by the distance metric of WF framework which is defined as follows [1]:

Definition 3.2**.**

Let $\bm{\rho}_{t}\in\mathbb{C}^{N}$ be an element of the global solution set in (2.1). The distance of an element $\bm{\rho}\in\mathbb{C}^{N}$ to $\bm{\rho}_{t}$ is defined as [1]

[TABLE]

For the purpose of convergence, the metric implies that we are primarily interested in the distance of an estimate $\bm{\rho}$ to any of the elements in the non-convex solution set $\mathit{P}$. In technical terms, invoking Definition 2.2, (27) is a measure of distance between the equivalence sets $\mathit{P}_{\bm{\rho}}$ and $\mathit{P}$. By Definition 3.2, the ambiguity due to the invariance of the cross-correlation map, $\mathcal{L}$, to the global phase factors is evaded on $\mathbb{C}^{N}$, without lifting the problem. Observe that the phase ambiguity is indeed removed analytically, since the $\ell_{2}$ norm is minimized when $\mathrm{Re}(\langle\bm{\rho},e^{\mathrm{i}\Phi(\bm{\rho})}\bm{\rho}_{t}\rangle)=|\langle\bm{\rho},e^{j\Phi(\bm{\rho})}\bm{\rho}_{t}\rangle|=|\langle\bm{\rho},\bm{\rho}_{t}\rangle|$, which is achieved at $e^{\mathrm{i}\Phi(\bm{\rho})}=\frac{\langle\bm{\rho}_{t},\bm{\rho}\rangle}{|\langle\bm{\rho},\bm{\rho}_{t}\rangle|}$. Hence, the squared distance becomes

[TABLE]

and is independent of any global phase factor on $\bm{\rho}$ or $\bm{\rho}_{t}$. Geometrically, the metric suggests that the cosine of the angle $\theta$ between the elements $\bm{\rho}$ and $\bm{\rho}_{t}$ is given by

[TABLE]

which can be viewed as the angle between the subspaces spanned by the elements $\bm{\rho}\bm{\rho}^{H}$ and $\bm{\rho}_{t}\bm{\rho}_{t}^{H}$ in the lifted domain $\mathbb{C}^{N\times N}$. This can be seen by evaluating the error between the lifted terms as follows:

[TABLE]

For two rank-1 arguments, the Frobenius inner product in (30) reduces to $|\langle\bm{\rho},\bm{\rho}_{t}\rangle|^{2}$, and the cosine of the angle between the elements becomes equal to $\cos^{2}(\theta)$. Since (29) is non-negative, the relationship between the two angles is one-to-one. Therefore, WF distance metric can be interpreted as a measure of distance between the lifted variables, and the convergence with respect to the metric (27) in signal domain is equivalent to convergence with respect to (30) in the lifted space.

3.3 Theoretical Advantages of GWF via the Lifted Perspective

The lifted perspective of GWF reveals its connection to convex low rank matrix recovery methods, such as those proposed in [48] and [11], as the objective in (23) simply corresponds to solving the quadratic data-fit over a more exclusive search space. In the scope of lifting based approaches, we discussed Uzawa’s method as a solver for the trace regularized problem in (13). Consider solving (23) by Uzawa’s method described in (14)-(15), which essentially reduces to projected gradient descent: a gradient step over the smooth $\ell_{2}$ mismatch term is followed by a projection onto the intersection of the PSD cone and the set of rank-1 matrices. While in general projections to the intersection of two sets is an optimization problem on its own, and the rank-1 constraint constitutes a non-convex manifold, there exists a simple projection onto the set $\mathcal{X}=\{\mathbf{X}\in\mathbb{C}^{N\times N}:\text{rank}(\mathbf{X})=1\cap\ \mathbf{X}\succeq 0\}$ such that

[TABLE]

which is the successive operation of the projection onto the PSD cone, followed by a rank-1 approximation.

Despite that a solver can be formulated, recovery guarantees of Uzawa’s method do not cover this case due to the non-convexity of the projected set. Simply, a gradient step over the convex PSD cone is not guaranteed to improve the rank-1 approximation of an estimate, which is after all the original motivation behind pursuing the convex solvers. Essentially, GWF framework circumvents convexification and presents an alternative update scheme to Uzawa’s for minimizing the objective in (23). This alternative update form stems from the fact PSD rank-1 constraint of the lifted unknown can be enforced by a variable transformation in the update equation, rather than by projections given in (31).

Clearly, the immediate advantage of the GWF formulation in (24) over the convex relaxations is the dimensionality reduction of the search space, attributed to iterating in the signal domain. There is yet another significant advantage to the GWF formulation relating to exact recovery guarantees. By (24), an iterative scheme for the unrelaxed, non-convex form of the lifted problem is formulated, which enforces the rank-1, PSD structure on the iterates. This allows the constraint set to be considerably smaller than that of the trace relaxation or the convex feasibility problems. Formally, the problem (24) has a unique solution if the equivalence set of any $\bm{\rho}\in\mathbb{C}^{N}$ is its equivalence class under the correlation map $\mathcal{L}$ as defined in Definition 2.2.

Condition 3.1**.**

Uniqueness Condition*. There exists a unique solution $\bm{\rho}_{t}\bm{\rho}_{t}^{H}\in\mathcal{X}$ for the problem in (24) if*

[TABLE]

In other words, there should exist no element $H$ in the null space of $\mathcal{F}$, such that $\bm{\rho}_{t}\bm{\rho}_{t}^{H}+H$ is a rank-1, PSD matrix. Therefore, for exact interferometric inversion by GWF, the null space condition of the lifted forward model has to hold over a much less restrictive set than any of the approaches discussed in Section 2.2.

3.4 The Spectral Method for GWF Initialization

Having to solve a non-convex problem, exact recovery guarantees of the WF framework depend on the accuracy of the initial estimate $\bm{\rho}_{0}$. The initial estimate of the standard WF algorithm is computed by the spectral method which corresponds to the leading eigenvector of the following positive semi-definite matrix:

[TABLE]

where $\mathbf{L}^{m}=\mathbf{L}_{i}^{m}=\mathbf{L}_{j}^{m}$, and $d^{m}=|\langle\mathbf{L}^{m},\bm{\rho}_{t}\rangle|^{2}$ for $i=j$. The leading eigenvector is scaled by the square root of the corresponding largest eigenvalue $\lambda_{0}$ of $\mathbf{Y}$. In [1], the spectral method is described from a stochastic perspective. By the strong law of large numbers, under the assumption that we have $\mathbf{L}^{m}\sim\mathcal{N}(0,\frac{1}{2}\mathbf{I})+\mathrm{i}\mathcal{N}(0,\frac{1}{2}\mathbf{I}),m=1,\cdots,M$, the spectral matrix $\mathbf{Y}$ becomes equal to

[TABLE]

as $M\rightarrow\infty$, which has the true solution $\bm{\rho}_{t}$ as its leading eigenvector. The concentration of the spectral matrix around its expectation is used to show that the leading eigenvector of $\mathbf{Y}$ is sufficiently accurate, such that the sequence of iterates $\{\bm{\rho}_{k}\}$ of (16) converges to an element in the global solution set $P$.

In developing the GWF framework, we view the spectral method as a procedure in the lifted domain. In fact, we observe that the spectral matrix of phase retrieval in (33) is the backprojection estimate of the lifted unknown $\tilde{\bm{\rho}}_{t}={\bm{\rho}}_{t}{\bm{\rho}}_{t}^{H}$. Having different measurement vectors $\mathbf{L}_{i}^{m}$ and $\mathbf{L}_{j}^{m}$ in the cross-correlated measurement case, using the definition of the backprojection operator in (21), we extend (33) and redefine $\mathbf{Y}$ as follows:

[TABLE]

As noted, the true solution $\tilde{\bm{\rho}}_{t}=\bm{\rho}_{t}\bm{\rho}_{t}^{H}$ of the lifted problem lies in the positive semi-definite (PSD) cone. In standard WF for phase retrieval, the spectral matrix $\mathbf{Y}$ is formed by summation of positive semi-definite outer products that are scaled by $\mathbb{R}^{+}$ valued measurements $\{d^{m}\}_{m=1}^{M}$, as auto-correlations are by definition squared magnitudes. Hence, the WF spectral method generates an estimate of the lifted unknown within the constraint set by default. This obviously is not the case for the backprojection estimate (35) with the cross-correlated measurement model. Therefore, the extension of the spectral method to cross-correlations includes a projection step onto the PSD cone. Since the PSD cone is convex, its projection operator is non-expansive and yields a closer estimate to $\bm{\rho}_{t}\bm{\rho}_{t}^{H}$ than that of (35). The GWF spectral matrix then becomes

[TABLE]

We discard the positive semi-definitivity in (36), and only project onto the set of symmetric matrices, which is also convex. This is simply because only the leading eigenvector will be kept from the generated lifted estimate, which is unaffected by the projection onto the PSD cone333Unless the leading eigenvalue is negative, a scenario that is excluded due to the conditions for exact recovery in Section 4.. Letting $\lambda_{0},\mathbf{v}_{0}$ denote the leading eigenvalue-eigenvector pair of $\hat{\mathbf{X}}$, the GWF initial estimate $\bm{\rho}_{0}$ is given by444Note that (36) yields a symmetric matrix, hence have eigenvalues $\lambda_{i}\in\mathbb{R}$.

[TABLE]

Using the representation in the lifted problem in (35) and plugging in (8) for the measurements, the GWF spectral matrix $\hat{\mathbf{X}}$ can be written as

[TABLE]

Hence, the leading eigenvalue-eigenvector extraction corresponds to keeping the rank-1, PSD approximation, $\tilde{\bm{\rho}}_{0}:=\bm{\rho}_{0}\bm{\rho}_{0}^{H}$, of the backprojection estimate in the lifted domain. Furthermore, by definition in (38), the accuracy of the spectral estimate fully hinges on the properties of the normal operator of the lifted problem, i.e., $\mathcal{F}^{H}\mathcal{F}$ over the set of rank-1, PSD matrices, $\mathcal{X}$.

As a comparison, consider the LRMR approach for the interferometric inversion by the PSD constrained singular value thresholding (SVT) algorithm in [11]. It can be seen that, setting $\mu_{k}=1/M$ in (15), the first iterate generated by (14) is $\hat{\mathbf{X}}=\mathcal{P}_{\tau}(\frac{1}{M}\mathcal{F}^{H}\mathbf{d})$, which, prior to the singular value thresholding, is identical to the backprojection estimate computed by the spectral method. Hence, the spectral method simply differs from the first Uzawa iteration by keeping the rank-1 approximation via the projection operator $\mathcal{P}_{\mathcal{X}}$ defined in (31), instead of a low-rank approximation. This precisely corresponds to the first Uzawa iteration to solve the non-convex rank-1 constrained problem in (23).

3.5 Algorithm Summary and Computational Complexity

Compared to the lifting based LRMR approaches, GWF provides significant reductions in computational complexity, and memory requirements per iteration. As shown in Section 3.4, GWF uses the first iteration of the Uzawa’s method to compute an initial estimate $\bm{\rho}_{0}$, and replaces the following iterations over the lifted domain with iterations on the leading eigenspace. The GWF algorithm is summarized as follows:

• •

Input: Interferometric measurements $d_{ij}^{m}$ and measurement vectors $\mathbf{L}_{i}^{m},\mathbf{L}_{j}^{m}$, $\forall i\neq j$, $m=1,\cdots,M$.

• •

Initialization: Run Uzawa’s method in (14) initialized with $\mathbf{y}=0$, $\mu_{0}=1/M$, and $\lambda=0$, i.e., trace regularization free, for $1$ iteration, yielding

[TABLE]

Keep the rank-1 approximation $\lambda_{0}\bm{\rho}_{0}\bm{\rho}_{0}^{H}$. The initialization step consists of the outer product of the two measurement vectors for each of the $M$ samples, resulting in $\mathcal{O}(MN^{2})$ multiplications, followed by an eigenvalue decomposition with $\mathcal{O}(N^{3})$ complexity.

• •

Iterations: Perform gradient descent updates as $\bm{\rho}_{k+1}=\bm{\rho}_{k}-\frac{\mu_{k+1}}{\|\bm{\rho}_{0}\|^{2}}\nabla\mathcal{J}(\bm{\rho}_{k})$, with

[TABLE]

where $(\mathbf{e}_{k})^{m}=(\left(\mathbf{L}_{i}^{m}\right)^{H}\bm{\rho}_{k}\bm{\rho}_{k}^{H}\mathbf{L}_{j}^{m}-d_{ij}^{m})$.

Each iteration requires the following operations:

Computing and storing the linear terms $(\mathbf{L}_{i,j}^{m})^{H}\bm{\rho}_{k}$, requiring $M$ number of $N$ multiplications for each, resulting in $\mathcal{O}(MN)$ multiplications. 2. 2.

Computing the error by cross correlating linear terms, requiring $\mathcal{O}(M)$ multiplications. 3. 3.

Multiplication of the linear terms $(\mathbf{L}_{i,j}^{m})^{H}\bm{\rho}_{k}$ and the error $e_{ij}^{m}$ for each $m=1,\cdots M$, requiring $\mathcal{O}(M)$ multiplications. 4. 4.

Multiplication of the result in $3$ with vectors $\{\mathbf{L}_{i}^{m}\}_{m=1}^{M}$ and $\{\mathbf{L}_{j}^{m}\}_{m=1}^{M}$, requiring $\mathcal{O}(MN)$ multiplications.

These operations result in $\mathcal{O}(MN)$ multiplications for each iteration.

4 Theory of Exact Interferometric Inversion via GWF

In this section we present our exact recovery guarantees for interferometric inversion by GWF. Notably, we merge the exact recovery conditions of standard Wirtinger Flow that rely on statistical properties of the sampling vectors into a single condition on the lifted forward model in the context of interferometric inversion ($i\neq j$). In Theorem 4.1, we state that if $\mathcal{F}$ satisfies the restricted isometry property on the set of rank-1, positive semidefinite matrices with a sufficiently small restricted isometry constant, GWF is guaranteed to recover the true solution upto a global phase factor, for any $\bm{\rho}_{t}\in\mathbb{C}^{N}$. Following Theorem 4.1, we establish the validity of our condition in Theorem 4.2 for the case of $\mathcal{O}(N\log N)$ measurements that are cross-correlations of i.i.d. complex Gaussian sampling vectors.

We begin by introducing the definitions of some concepts that appear in our theorem statements.

4.1 The $\epsilon$-Neighborhood of $\mathit{P}$ and the Regularity Condition

Definition 4.1**.**

$\epsilon$-Neighborhood of $\mathit{P}$. We denote the $\epsilon$-neighborhood of the global solution set $\mathit{P}$ in (2.1) by $\mathit{E}(\epsilon)$ and define it as follows [1]:

[TABLE]

The set $E(\epsilon)$ is determined by the distance of the spectral initialization to the global solution set, i.e., $\epsilon=\text{dist}(\bm{\rho}_{0},\bm{\rho}_{t})$. The main result of the standard WF framework is that, for Gaussian and coded diffraction pattern measurement models [1], $\epsilon$ is sufficiently small so that the objective function $\mathcal{J}$ satisfies the following regularity condition:

Condition 4.1**.**

Regularity Condition*. The objective function $\mathcal{J}$ satisfies the regularity condition if, for all $\bm{\rho}\in\mathit{E}(\epsilon)$, the following holds*

[TABLE]

for fixed $\alpha>0$ and $\beta>0$ such that $\alpha\beta>4$.

The regularity condition guarantees that the iterations in (19) are contractions with respect to the distance metric (27), which ensures that all gradient descent iterates remain in $E(\epsilon)$, which is established by the following lemma from [1]:

Lemma 4.1**.**

Assume that $\mathcal{J}$ obeys the regularity condition in (40) for some fixed $\alpha,\beta$ for all $\bm{\rho}\in E(\epsilon)$. Having $\bm{\rho}_{0}\in E(\epsilon)$, and assuming $\mu\leq{2}/{\beta}$, consider the following update

[TABLE]

Then, for all $k$ we have $\bm{\rho}_{k}\in E(\epsilon)$ and

[TABLE]

Proof.

See [1], Lemma 7.10. ∎

Furthermore, from the definition of $\nabla\mathcal{J}$ in (19), the regularity condition implies that there exists no $\bm{\rho}\in E(\epsilon)$ that belongs to the equivalence class of $\bm{\rho}_{t}$ under the mapping $\mathcal{L}$. Hence, the uniqueness condition for exact recovery is satisfied locally, and 4.1 is a sufficient condition by Lemma 7.1 of [1] such that the algorithm iterates $\{\bm{\rho}_{k}\}$ converge to $\mathit{P}$ at a geometric rate. The spectral initialization is said to be sufficiently accurate, if (40) holds $\forall\bm{\rho}\in\mathit{E}(\epsilon)$.

4.2 Sufficient Conditions For Exact Recovery

The condition we assert on the lifted forward model $\mathcal{F}$ is the restricted isometry property on the set of rank-1, PSD matrices, $\mathcal{X}$.

Definition 4.2**.**

Restricted Isometry Property (RIP). Let $\mathcal{A}:\mathbb{C}^{K\times N}\rightarrow\mathbb{C}^{M}$ denote a linear operator. Without loss of generality assume $K\leq N$. For every $1\leq r\leq K$, the $r-\text{restricted}$ isometry constant (RIC) is defined as the smallest $\delta_{r}<1$ such that

[TABLE]

holds for all matrices $\mathbf{X}$ of rank at most $r$, where $\|\mathbf{X}\|_{F}=\sqrt{\text{Tr}(\mathbf{X}^{H}\mathbf{X})}$ denotes the Frobenius norm.

Suppose (42) holds for all $\mathbf{X}\in\mathcal{D}\subset\mathbb{C}^{K\times N}$ that have rank-r with some constant $0<\delta_{r}<1$, then $\mathcal{A}$ is said to satisfy the RIPD with RIC-$\delta_{r}$. For the interferometric inversion problem, having $K=N$, if there exists $\delta_{1}<1$ such that $\mathcal{F}$ in (8) satisfies the restricted isometry property on the PSD cone, we say that the lifted forward model satisfies RIP+ with RIC-$\delta_{1}$, i.e., RIP over the set of rank-1, PSD matrices, $\mathcal{X}$.

(42) quantifies how close $\mathcal{F}^{H}\mathcal{F}$ is to an identity over rank-1, PSD matrices through the following lemma:

Lemma 4.2**.**

Suppose $\mathcal{F}$ satisfies RIP on the set rank-1, PSD matrices of size $N\times N$, i.e., $\mathcal{X}=\{\bm{\rho}\bm{\rho}^{H}:\bm{\rho}\in\mathbb{C}^{N}\}$ with RIC-$\delta_{1}$. Then, for any $\mathbf{X}\in\mathcal{X}$ we have

[TABLE]

where $\bm{\delta}:\mathcal{X}\rightarrow\mathbb{C}^{M}$ is a bounded operator such that

[TABLE]

where $\|\cdot\|$ denotes the spectral norm.

Proof.

See Appendix A.2. ∎

In the GWF framework, we identify RIP over rank-1, PSD matrices with RIC-$\delta_{1}\leq 0.214$ as a sufficient condition for an arbitrary lifted forward model $\mathcal{F}$, which guarantees that the spectral initialization of GWF provides an initial estimate that is sufficiently accurate, i.e., GWF iterates are guaranteed to converge to a global solution in $\mathit{P}$ starting from $\bm{\rho}_{0}$.

Theorem 4.1**.**

Exact Recovery by GWF*. Assume the lifted forward model $\mathcal{F}$555Upto a normalization factor, such as $\frac{1}{\sqrt{M}}$. satisfies the RIP condition over rank-1, PSD matrices with the restricted isometry constant (RIC)-$\delta_{1}$. Then, the initial estimate $\bm{\rho}_{0}$ obtained from the spectral method by (36) and (37) satisfies*

[TABLE]

where $\epsilon$ is an $\mathcal{O}(1)$ constant with $\epsilon^{2}=(2+\delta_{1})(1-\sqrt{1-\frac{\delta_{1}}{1-\delta_{1}}})+\frac{\delta_{1}^{2}}{8}$, and the regularity condition in (40) surely holds for the objective function in (3) if $\delta_{1}\leq 0.214$, with any $\alpha$, $\beta$ satisfying

[TABLE]

where $\delta_{2}=\frac{\sqrt{2}(2+\epsilon)\delta_{1}}{\sqrt{(1-\epsilon)(2-\epsilon})}$ and $c(\delta_{1})=(2+\epsilon)(1+\epsilon)(1+\delta_{1})$. Thus, for the iterations of (41) with the update term in (19) and the fixed step size $\mu\leq{2}/{\beta}$, we have

[TABLE]

Proof.

See Section 5.1. ∎

We refer to $\delta_{2}$ as the restricted isometry constant of the local RIP-2 condition by Lemma 5.5 which is exclusively over elements of the form $\{\bm{\rho}\bm{\rho}^{H}-\bm{\rho}_{t}\bm{\rho}_{t}^{H}:\bm{\rho}\in E(\epsilon)\}$ , and $c$ as the local Lipschitz constant of $\nabla\mathcal{J}$ by Lemma 5.6, both stated in Section 5.

Remark**.**

We summarize the implications of Theorem 4.1 by the following remarks.

Theorem 4.1 establishes a regime in which the regularity condition holds by default as a result of the RIP over rank-1, PSD matrices, captured by (45). This regime also defines the range of values the RIC-$\delta_{1}$ can attain, through the quantity $\delta_{2}$ which must satisfy $\delta_{2}<1$, as shown in Figure 1. Notably, having $\delta_{2}<1$ guarantees the uniqueness of a solution locally in $E(\epsilon)$, and the restricted strong convexity of $\mathcal{J}$ since there exists a fixed $\alpha,\beta$ satisfying (45). Numerically, plugging in the $\epsilon$ constant defined by $\delta_{1}$, this is satisfied for $\delta_{1}\leq 0.214$. 2. 2.

Figure 1 demonstrates the values the constants $c$ and $h$ attain in the valid range for $\delta_{1}$. These $\mathcal{O}(1)$ constants directly impact the convergence rate of the algorithm, as $\alpha$ and $\beta$ are required to be sufficiently large values for (45) to hold. Observe that (45) implies setting $\alpha=\alpha^{\prime}/{\|\bm{\rho}_{t}\|^{2}}$, and $\beta=\beta^{\prime}\|\bm{\rho}_{t}\|^{2}$, where $\alpha^{\prime},\beta^{\prime}=\mathcal{O}(1)$, hence $\alpha\beta=\mathcal{O}(1)$. Since clearly $h<2$, and $c^{2}>4$, and the regularity condition can at best hold when $\alpha^{\prime}\beta^{\prime}>4$ by definition. Therefore, the RIP over rank-1, PSD matrices with RIC-$\delta_{1}\leq 0.214$ is a sufficient condition for the exact recovery via GWF. 3. 3.

Since $\|\bm{\rho}_{t}\|$ is unknown a priori, a suitable approximate scaling in setting constants $\alpha,\beta$ is $\|\bm{\rho}_{0}\|^{2}$. This scaling factor is merely the leading eigenvalue $\lambda_{0}$ of the spectral matrix in (36). From Lemma 5.2, for $\mathcal{F}$ satisfying RIP over rank-1, PSD matrices with the RIC-$\delta_{1}$, $\lambda_{0}$ is lower bounded by $1-\delta_{1}$, hence the condition in (45) can be enforced by setting $\alpha^{\prime},\beta^{\prime}=\mathcal{O}(1)$ such that

[TABLE]

Similarly, picking $\mu\leq\frac{2}{\beta}$ to yield convergence rate of $\frac{2\mu}{\alpha}\leq\frac{4}{\alpha^{\prime}\beta^{\prime}}$, the iterations in (41) must have a step size $\mu$ that is $\mathcal{O}({1}/{\|\bm{\rho}_{t}\|^{2}})$. This is essentially where the normalization term in (16) originates from, with $\|\bm{\rho}_{0}\|^{2}$ serving as an approximation to $\|\bm{\rho}_{t}\|^{2}$. As a result, Lemma 4.1 for (16) simply holds for $\mu_{k}\leq(1-\delta_{1})\frac{2}{\beta^{\prime}}$ due the effect of the mismatch in the scaling factors, in agreement to what is noted in [1]. 4. 4.

Overall, Theorem 4.1 establishes that the convergence speed of GWF algorithm is controlled by the restricted isometry constant $\delta_{1}$. As $\delta_{1}$ approaches to the critical limit of $0.214$, $c$ and $\delta_{2}$ values increase super-linearly as shown in Figure 1. This has strong implications on the convergence speed, as $\beta^{\prime}$ is inversely proportional to the choice of step sizes $\mu_{k}$, and quadratically related to the magnitude of $c$. 5. 5.

A consequence of our result is a universal upper bound on $\epsilon$ under our sufficient condition of RIP over rank-1, PSD matrices. As depicted in Figure 1, Theorem 4.1 determines what sufficiently close means numerically.

Remark**.**

Despite solving the identical perturbed problem, GWF iterations provably converge, whereas the convergence guarantees of Uzawa’s method vanish due to inclusion of a non-convex constraint. Similarly, the special structure of the constraint set and the GWF iterates suffice the RIP condition to be satisfied only over rank-1 matrices in the PSD cone, whereas the uniqueness condition of Uzawa’s method requires RIP over the set of rank-2 matrices in the non-convex rank minimization problem.

4.3 Restricted Isometry Property for Cross-Correlation of Gaussian Measurements

In standard WF for phase retrieval with the i.i.d. complex Gaussian model, i.e., $\mathbf{L}_{m}\sim\mathcal{N}(0,\mathbf{I}/2)+\mathrm{i}\mathcal{N}(0,\mathbf{I}/2)$, the accuracy of the spectral estimate (33) is established as

[TABLE]

with probability $1-10e^{\gamma N}-8/N^{2}$ where $\gamma$ is a fixed positive numerical constant. This result is derived from the concentration bound of the Hessian of the objective function around its expectation at a global minimizer $\bm{\rho}_{t}$ such that

[TABLE]

where $\delta$ is the concentration bound.

In the problem of phase retrieval, plugging in the definition of phaseless measurements into (33), the auto-correlations yield the $4^{th}$ moments of the elements of the Gaussian measurement vectors. This introduces a bias of $\|\bm{\rho}_{t}\|^{2}\mathbf{I}$ in the spectral estimate $\mathbf{Y}$ as can be seen in (47). Moving from the auto-correlations to cross-correlations removes this bias component from the spectral matrix of GWF. Hence, for the Gaussian model, cross-correlated measurement map, i.e., $\mathcal{F}$ when $i\neq j$, satisfies the RIP over the set of rank-1, PSD matrices.

Without loss of generality and following [1], we present our result for the case $\|\bm{\rho}_{t}\|=1$.

Theorem 4.2**.**

RIP over Rank-1, PSD Matrices for Cross-Correlated Gaussian Measurements.* Let the measurement vectors $\mathbf{L}_{i}^{m},\mathbf{L}_{j}^{m}$ in (2) follow the i.i.d. complex Gaussian model, i.e. $\mathbf{L}_{i}^{m},\mathbf{L}_{j}^{m}\sim\mathcal{N}(0,\mathbf{I}/2)+\mathrm{i}\mathcal{N}(0,\mathbf{I}/2)$. Then, the lifted forward model $\frac{1}{\sqrt{M}}\mathcal{F}$ in (8) for cross-correlated measurements (when $i\neq j$), satisfies RIP defined in (42) for $r=1$ over the PSD cone, with probability $1-8e^{-\gamma N}-5/N^{2}$ given $\mathcal{O}(N\log N)$ measurements, where $\gamma$ is a fixed positive numerical constant. Moreover, the spectral matrix $\hat{\mathbf{X}}$ defined in (36) satisfies*

[TABLE]

*where $\bm{\rho}_{t}$ is the ground truth signal with $\|\bm{\rho}_{t}\|=1$ and $\delta_{1}$ is the restricted isometry constant of $\frac{1}{\sqrt{M}}\mathcal{F}$. *

Proof.

See Section 5.2. ∎

Remark**.**

Outcomes of Theorem 4.2 are explained with the following remarks.

Theorem 4.2 establishes the relationship between the concentration bound of the spectral matrix, and the RIP-1 condition for interferometric inversion. This indicates that the regularity condition of WF framework is redundant for our problem if $\delta_{1}$ is picked properly, since by Theorem 4.1, RIP+ with RIC-$\delta_{1}\leq 0.214$ directly implies the regularity condition. 2. 2.

Note that the equivalent linear model in the lifted domain actually has $N^{2}$ unknowns. By our formulation of GWF in the lifted problem, having measurements of the order of $N\log N$, therefore, corresponds to an under-determined system of equations in which exact recovery guarantees of GWF hold. 3. 3.

Note that our measurement complexity is identical to that of standard WF. This is an expected result, as the cross-correlations only impact the removal of the diagonal bias in the spectral matrix $\mathbf{Y}$, not the concentration of $\mathbf{Y}$ around its expectation. 4. 4.

It should be re-iterated that when $i=j$, the backprojection estimate $\mathbf{Y}$ contains a diagonal bias of the form $\|\bm{\rho}_{t}\|^{2}\mathbf{I}$, which breaks our RIP condition over rank-1, PSD matrices. In such a case, the concentration bound then holds around the expectation term that includes this diagonal bias. Establishing a sufficient condition for the phaseless case under our framework is not in the scope of this paper, which we omit for future work.

Remark**.**

It should be noted that phaseless measurements in the Gaussian model is known to satisfy a restricted isometry property under the $\ell_{1}$ norm in the range of the lifted forward map. This is referred to as RIP-1 condition in literature, and is prominently featured in exact recovery theory of standard WF, and more recently in [45] for non-convex optimization via gradient descent, even in the presence of non-convex regularizers.

A critical difference of our work is that such a RIP type condition isn’t used as a tool to obtain local curvature and local smoothness conditions to assert the regularity condition with a high probability. Instead, under the RIP studied in this paper, a regime is derived in which the regularity condition is guaranteed to hold, deterministically. Notably, for the interferometric inversion problem with $i\neq j$, $\delta_{1}$ values tested in [1] are confidently within the range in which our GWF theory applies.

5 Proofs of Theorems 4.1 and 4.2

In this section, we present the proofs of Theorems 4.1 and 4.2. We first present key lemmas, and next prove the theorems using these lemmas. We provide the detailed proofs of the lemmas in Appendix B and C.

To prove Theorem 4.1, we begin by showing that the RIC-$\delta_{1}$ of our RIP condition determines the distance $\epsilon$ of the spectral initialization in a one-to-one manner. We then establish that the regularity condition is directly implied by RIP+ with RIC-$\delta_{1}\leq 0.214$. In achieving this result we first show that the structure of the rank-1, PSD set allows for restricted isometry property to hold locally for the difference of two rank-1 PSD matrices, i.e., a local RIP-2 condition similar to the one in [66], with a RIC-$\delta_{2}$. The upper bound on $\delta_{1}$ ensures that RIC-$\delta_{2}$ of the local RIP-2 satisfies $\delta_{2}<1$. The local RIP-2 condition, in turn, ensures that restricted strong convexity holds in the $\epsilon$-neighborhood of the global solution set, which leads to exact recovery conditions of GWF.

For Theorem 4.2, we first show that the bias term in (47) resulting from the $4^{th}$ moments of the random Gaussian entries disappear when we have cross-correlations instead of auto-correlations of measurements. We then establish that the spectral matrix is concentrated around its expectation, using the machinery in [1], adapted for cross-correlations. Finally, we use the definition of the spectral matrix to derive the RIP over rank-1, PSD matrices from the concentration bound, which yields the RIC-$\delta_{1}$.

5.1 Proof of Theorem 4.1

Without the loss of generality, we assume $\bm{\rho}_{t}$ is a solution with $\|\bm{\rho}_{t}\|=1$. In establishing the exact recovery guarantees for GWF, we take a two-step approach. For a lifted forward map $\mathcal{F}$ satisfying RIP+ with RIC-$\delta_{1}$ , we first show that the initialization by spectral method yields an estimate that is in the set $E(\epsilon)$. We then establish the regularity condition (40) for the objective function (3) in the $\epsilon$-neighborhood defined by the initialization. These two results culminate into convergence to a global solution at a geometric rate as stated in Theorem 4.1.

5.1.1 $\epsilon$-Neighborhood of Spectral Initialization

Rather than the law of large numbers approach in [1], we take the geometric point of view of [43] in establishing the $\epsilon$-neighborhood of the spectral initialization. We begin by evaluating the distance of the leading eigenvector $\mathbf{v}_{0}\in\mathbb{C}^{N}$ of the spectral matrix in (36) to the global solution set (2.1). Recall the definition of the distance metric

[TABLE]

which is essentially the Euclidean distance of $\mathbf{v}_{0}$ to the closest point in the solution set $\mathit{P}$ in (2.1). Without loss of generality, we fix $\Phi(\mathbf{v}_{0})=\Phi_{0}$ and incorporate it into $\bm{\rho}_{t}$ such that $\hat{\bm{\rho}}_{t}=e^{\mathrm{i}\Phi_{0}}\bm{\rho}_{t}$ represents the closest solution to $\mathbf{v}_{0}$ in $\mathit{P}$. We breakdown the key arguments of our proof into the following three lemmas.

Lemma 5.1**.**

Let $\bm{\rho}_{t}$ be a solution with $\|\bm{\rho}_{t}\|=1$ and $\hat{\bm{\rho}}_{t}$ is the closest solution in $\mathit{P}$ to $\mathbf{v}_{0}$. Then, $\mathrm{Re}\langle\hat{\bm{\rho}}_{t},\mathbf{v}_{0}\rangle=\langle\hat{\bm{\rho}}_{t},\mathbf{v}_{0}\rangle$, and

[TABLE]

where $\|\mathbf{v}_{0}\|=1$, $\cos(\theta)=\mathrm{Re}\langle\hat{\bm{\rho}}_{t},\mathbf{v}_{0}\rangle$, and $\mathbf{v}_{0}^{\perp}$ is a unit vector lying in a plane whose normal is $\mathbf{v}_{0}$. Similarly, there exists a perpendicular unit vector, $\hat{\bm{\rho}}_{t}^{\perp}$, to $\hat{\bm{\rho}}_{t}$ such that

[TABLE]

Proof.

See Appendix B.1. ∎

Lemma 5.2**.**

Consider the spectral matrix $\hat{\mathbf{X}}$ given by (38), and denote the spectral matrix projected onto the positive semi-definite cone as $\hat{\mathbf{X}}_{PSD}$. Then, for a lifted mapping $\mathcal{F}$ satisfying RIP+ with RIC-$\delta_{1}=\delta$, $\hat{\mathbf{X}}$ and $\hat{\mathbf{X}}_{PSD}$ have the identical leading eigenvalue-eigenvector pair $\lambda_{0}$, $\mathbf{v}_{0}$ such that

[TABLE]

Furthermore, $\hat{\mathbf{X}}$ and $\hat{\mathbf{X}}_{PSD}$ generate identical spectral initialization, $\bm{\rho}_{0}$.

Proof.

See Appendix B.2. ∎

Lemma 5.2 allows us to analyze the distance of the initial estimate $\bm{\rho}_{0}$ to the solution set by the convenience of either the positive semi-definite $\hat{\mathbf{X}}_{PSD}$ or the symmetric spectral estimate $\hat{\mathbf{X}}$, since they generate the same initial estimate $\bm{\rho}_{0}$.

Next, using the Lemmas 5.1 and 5.2 we reach the following key result.

Lemma 5.3**.**

In the setup of Lemmas 5.1 and 5.2, for the angle $\theta$ between the one-dimensional sub-spaces spanned by $\hat{\bm{\rho}}_{t}$ and $\mathbf{v}_{0}$ we have

[TABLE]

where $\delta$ is the RIC-$\delta_{1}$ of the lifted map $\mathcal{F}$ satisfying RIP+ .

Proof.

See Appendix B.3. ∎

From Lemma 5.3, we can now lower bound the inner product of $\hat{\bm{\rho}}_{t}$ and $\mathbf{v}_{0}$ such that

[TABLE]

where $\kappa=\frac{\delta}{1-\delta}$. Writing the distance of the spectral initialization $\bm{\rho}_{0}=\sqrt{\lambda_{0}}\mathbf{v}_{0}$ to the solution set, we have

[TABLE]

It is easy to see that $\text{Re}\langle e^{\mathrm{i}\Phi(\bm{\rho}_{0})}{\bm{\rho}}_{t},\sqrt{\lambda_{0}}\mathbf{v}_{0}\rangle$ is maximized when $\Phi(\bm{\rho}_{0})=\Phi(\mathbf{v}_{0})$, hence we get

[TABLE]

From Lemma 5.2, we know that $1-\delta\leq\lambda_{0}\leq 1+\delta$. Moreover, the upper bound on the right hand side of (54) is simply a quadratic term with respect to $\sqrt{\lambda_{0}}$ since $\sqrt{\lambda_{0}}(\sqrt{\lambda_{0}}-2\sqrt{1-\kappa})+1$, which is maximized at the boundary of the domain of values $\sqrt{\lambda_{0}}$ takes. Since the quadratic equation is minimized at $\sqrt{1-\kappa}$ and we have ${1-\kappa}\leq{1-\delta}$ for the domain of possible values of $0\leq\delta<1$, $\lambda_{0}=1+\delta$ is an upper bound for the right-hand-side of (54). Hence, we obtain

[TABLE]

Writing the Taylor series expansion of $\sqrt{1+\delta}$ around [math], and discarding the components of order $\mathcal{O}(\delta^{3})$ and higher, we have the final upper bound

[TABLE]

which sets the $\epsilon$-neighborhood as

[TABLE]

5.1.2 Proof of the Regularity Condition

Recall that we seek a solution to the interferometric inversion problem by minimizing the following loss function

[TABLE]

and address the optimization by forming the steepest descent iterates

[TABLE]

where $\mu$ is the learning rate and $\nabla\mathcal{J}$ is the complex gradient defined by the Wirtinger derivatives.

As shown in Section 3.1.2, the gradient evaluated at a point $\bm{\rho}$ can be expressed as

[TABLE]

where

[TABLE]

with $\tilde{\bm{\rho}}$ and $\tilde{\bm{\rho}}_{t}$ denoting the lifted variables $\tilde{\bm{\rho}}=\bm{\rho}\bm{\rho}^{H}$ and $\tilde{\bm{\rho}}_{t}=\bm{\rho}_{t}\bm{\rho}_{t}^{H}$, respectively. Invoking Lemma 4.2 and the linearity of $\mathcal{F}^{H}\mathcal{F}$ and $\bm{\delta}$ of (43), (60) can be represented as

[TABLE]

where $\tilde{\mathbf{e}}=\tilde{\bm{\rho}}-\tilde{\bm{\rho}_{t}}$ is the error in the lifted problem. Since the lifted variables are already symmetric we can take them out of the projection operator due its linearity. Hence, for the update term we obtain

[TABLE]

Reprising the regularity condition under consideration, we need to establish that there exists constants $\alpha$ and $\beta$, such that $\alpha\beta>4$ for all $\bm{\rho}\in E(\epsilon)$ and

[TABLE]

To show the existence of constants $\alpha$ and $\beta$ that satisfy (63), we upper bound the gradient term, which converts the regularity condition to a restricted strong convexity condition [67, 68]. We begin the proof by introducing the following key lemmas.

Lemma 5.4**.**

Let $\bm{\rho}_{t}$ be the ground truth signal with $\|\bm{\rho}_{t}\|=1$, and $\hat{\bm{\rho}}_{t}$ denote the global solution closest to $\bm{\rho}$ such that $\hat{\bm{\rho}}_{t}=e^{\mathrm{i}\Phi(\bm{\rho})}\bm{\rho}_{t}$. Then, for any $\bm{\rho}\in E(\epsilon)$, we have

[TABLE]

Proof.

See Appendix B.4. ∎

Lemma 5.5**.**

Let $\bm{\rho}_{t}$ be the ground truth signal with $\|\bm{\rho}_{t}\|=1$, and let the linear map $\mathcal{F}$ satisfy RIP+ with RIC-$\delta_{1}$. Then, for $\delta_{1}\leq 0.214$ and any $\bm{\rho}\in E(\epsilon)$, we have $\delta_{2}<1$ such that

[TABLE]

where $\delta_{2}=\frac{\sqrt{2}(2+\epsilon)}{\sqrt{(1-\epsilon)(2-\epsilon})}\delta_{1}$.

We refer to (65) as the local RIP-2 condition in the lifted domain with RIC-$\delta_{2}$ for the mapping $\mathcal{F}$. The two lemmas culminate into the local Lipschitz continuity of $\nabla\mathcal{J}$.

Proof.

See Appendix B.5. ∎

Lemma 5.6**.**

In the setup of Lemmas 5.4 and 5.5, for any $\bm{\rho}\in E(\epsilon)$, the objective function $\mathcal{J}$ in (57) is Lipschitz differentiable with

[TABLE]

where $c=(1+\epsilon)(2+\epsilon)(1+\delta_{1})$ is the Lipschitz constant. Furthermore, to establish the regularity condition for $\mathcal{J}$, it is sufficient to show that

[TABLE]

for any $\bm{\rho}\in E(\epsilon)$.

Proof.

See Appendix B.6. ∎

We finally utilize the following lemma to obtain an alternative form of the restricted strong convexity condition [69].

Lemma 5.7**.**

For the objective function $\mathcal{J}$ in (57), the condition in (67) is satisfied if

[TABLE]

where $\eta=\frac{1}{\alpha}+\frac{c^{2}}{\beta}$.

Proof.

See Appendix B.7. ∎

Writing the objective function explicitly in terms of the lifted terms, and applying the lower bound from the local RIP-2 condition of the lifted forward model, we can express the regularity condition simply as

[TABLE]

where $\delta_{2}$ is as defined in Lemma 5.5. From Lemma 5.4, the regularity condition is then satisfied by identifying $\alpha,\beta$ with $\alpha\beta>4$ such that

[TABLE]

where, from Lemma 5.6, $c=(1+\epsilon)(2+\epsilon)(1+\delta_{1})$.

5.2 Proof of Theorem 4.2

Lemma 5.8**.**

Expectation of Spectral Matrix. Let measurement vectors $\mathbf{L}_{i}^{m},\mathbf{L}_{j}^{m}$ be statistically independent and distributed according to the complex Gaussian model as $\mathbf{L}_{i},\mathbf{L}_{j}^{m}\sim\mathcal{N}(0,\frac{1}{2}\mathbf{I})+\mathrm{i}\mathcal{N}(0,\frac{1}{2}\mathbf{I})$. Let $\bm{\rho}_{t}$ be independent of the measurement vectors, and $\mathbf{Y}$ denote the backprojection estimate of the lifted signal generated by the spectral method given as

[TABLE]

where $\mathcal{F}$ is the lifted forward map in (8). Then,

[TABLE]

Proof.

See Appendix C.1. ∎

Lemma 5.9**.**

Concentration Around Expectation. In the setup of Lemma 5.8, assume that the number of measurements is $M=C(\delta)\cdot N\log N$, where $C$ is a constant that depends on $\delta$. Then,

[TABLE]

holds with probability at least $p=1-8e^{-\gamma N}-5N^{-2}$ where $\gamma$ is a fixed positive constant.

Proof.

See Appendix C.2. ∎

Plugging in the expectation from Lemma 5.8 into (69), and using the definition of $\bm{\delta}$ from Lemma 4.2, the Lemmas 5.8 and 5.9 culminate to

[TABLE]

for any $\bm{\rho}_{t}$ with $\|\bm{\rho}_{t}\|=1$ with probability at least $p$. From the definition of the spectral norm on $\mathbb{C}^{N\times N}$, we can write the left-hand-side of (70) equivalently as

[TABLE]

Since the spectral norm corresponds to the maximum over the unit sphere in $\mathbb{C}^{N}$, we have

[TABLE]

Observe that the left-hand-side can equivalently be represented as a Frobenius inner product via lifting as

[TABLE]

Having Lemmas 5.8 and 5.9 hold for any element $\bm{\rho}_{t}\in\mathbb{C}^{N}$ with $\|\bm{\rho}_{t}\|=1$ via unitary invariance, for any $\bm{\rho}\in\mathbb{C}^{N}$ we obtain

[TABLE]

which yields

[TABLE]

for any $\bm{\rho}\in\mathbb{C}^{N}$. Therefore, RIP+ with RIC-$\delta$ is established with probability at least $p$ for mapping $\frac{1}{\sqrt{M}}\mathcal{F}$ with $\mathbf{L}_{i}^{m},\mathbf{L}_{j}^{m}\sim\mathcal{N}(0,\frac{1}{2}\mathbf{I})+\mathrm{i}\mathcal{N}(0,\frac{1}{2}\mathbf{I})$ where $M=C(\delta)\cdot N\log N$.

Furthermore, we know that the true lifted unknown $\bm{\rho}_{t}\bm{\rho}_{t}^{H}$ lies in the positive semi-definite cone, which is a convex set. Since the spectral matrix $\hat{\mathbf{X}}$ is the projection of $\frac{1}{M}\mathcal{F}^{H}\mathcal{F}(\bm{\rho}_{t}\bm{\rho}_{t}^{H})$ onto the set of Hermetian symmetric matrices, from the non-expansiveness property of projections onto convex sets we have

[TABLE]

which completes the proof of Theorem 4.2.

6 Numerical Simulations

6.1 Signal Recovery From Random Gaussian Measurements

We begin by considering recovery of random signals from cross-correlations of complex random Gaussian measurements, $\mathbf{L}^{m}_{i},\mathbf{L}^{m}_{j}\sim\mathcal{N}(\mathbf{0},\frac{1}{2}\mathbf{I})+\mathrm{i}\mathcal{N}(\mathbf{0},\frac{1}{2}\mathbf{I})$. For our numerical evaluations of the Gaussian model, we conduct an experiment similar to that of [1]. We set $N=128$, and run $100$ instances of interferometric inversion by GWF with independently sampled Gaussian measurement vectors on two types of signals: random low-pass signals, $\bm{\rho}^{LP}$, and random Gaussian signals, $\bm{\rho}^{G}$. The entries of the signals are generated independently of the measurement vectors at each instance by

[TABLE]

where $P=N/8$, and $X_{l}$ and $Y_{l}$ are i.i.d. $\mathcal{N}(0,1)$. As described in [1], random low-pass signal corresponds to a bandlimited version of this random model and variances are adjusted so that the expected signal power is the same.

We implement the GWF algorithm with the learning rate heuristic of WF in [1] such that the descent algorithm takes smaller steps initially due to higher inaccuracy of the iterates. The step size is gradually increased such that $\mu_{k}=\text{min}(1-\mathrm{e}^{k/\tau_{0}},\mu_{max})$, where $\tau_{0}=33000$, and $\mu_{max}=0.2$. For $2500$ iterations, the learning parameter corresponds to a nearly linear regime and attains the maximum value of $0.073$.

In the experimentation, we compute the empirical probability of success after $2500$ iterations by counting the exact recovery instances of GWF recovery from different realizations of Gaussian measurements for $\{\mathbf{L}_{i}^{m},\mathbf{L}_{j}^{m}\}_{m=1}^{M}$. We evaluate the exact recovery by the relative normalized error of the final estimate, $\bm{\rho}_{GWF}$, such that $\text{dist}(\bm{\rho}_{GWF},\bm{\rho}_{t})/\|\bm{\rho}_{t}\|\leq err=10^{-5}$. In addition, we evaluate the probability of moderately precise recovery by setting $err=10^{-3}$. As shown in Figure 2, our experimentation indicates that beginning with $3N$ interferometric Gaussian measurements, GWF achieves exact recovery with high probabilities. Furthermore, the method provides robust recovery with as low as $2.3N$ interferometric Gaussian measurements, with a relative error of $10^{-3}$ and below at over $95$ percent.

6.2 Multistatic Passive Radar Imaging

An interferometric inversion problem of great interest is multistatic passive radar imaging. We consider an imaging set-up in which several static, terrestrial receivers are placed in a circle of radius around the scene of interest, which is illuminated by a transmitter of opportunity. An exemplary multistatic imaging geometry is illustrated in Figure 3. At receiver $i$, the back-scattered signal is collected at a fixed location by a linear map $\mathbf{L}_{i}$ parameterized by the temporal frequency variable $\omega$. The linear measurements collected at two receivers $i$ and $j$ are then pairwise correlated in time to yield the cross-correlation model (2) defined in the temporal frequency domain.

The key advantage of the interferometric model is the elimination of the dependence of measurements on the transmitter location and phase of the transmitted waveform, both of which are unknown in the passive scenario. In prior studies, the interferometric wave-based imaging was approached by low rank recovery methods [12, 11, 70]. We postulate that GWF framework provides a computationally and memory-wise efficient alternative to LRMR based passive radar imaging.

6.2.1 Received Signal Model

Let $\mathbf{a}_{i}^{r}\in\mathbb{R}^{3}$ denote the spatial locations of the receivers, and assume $S$ number of receivers such that $i=1,2,\cdots S$. We assume that scattered signals are due to a single source of opportunity located at $\mathbf{a}^{t}$. The location on the surface of the earth is denoted by $\mathbf{x}=(\bi x,\psi(\bi x))\in\mathbb{R}^{3}$, where $\bi x=(x_{1},x_{2})\in\mathbb{R}^{2}$ and $\psi:\mathbb{R}^{2}\rightarrow\mathbb{R}$ is a known ground topography; and ${\rho}:\mathbb{R}^{2}\rightarrow\mathbb{R}$ denotes the ground reflectivity.

Under the flat topography assumption and Born approximation, and assuming waves propagate in free space, the fast-time temporal Fourier transform of the received signal at the $i$-th receiver can be modeled as [5]

[TABLE]

where $\omega$ is the temporal frequency variable,$c_{0}$ is the speed-of-light in free space, $p(\omega,s)$ is the transmitted waveform, $\alpha_{i}(\bi x,\mathbf{a}^{t})$ is the azimuth beam pattern, and

[TABLE]

is the bi-static delay term.

Following the derivations in [11] and discretizing the domain $D$ of $\rho$ into $N$ samples, the cross-correlated data model for the multi-static scenario is obtained as

[TABLE]

where $\bm{\rho}_{k}=\rho(\mathbf{x}_{k})$ is the $k^{th}$ entry of the discretized scene reflectivity vector $\bm{\rho}\in\mathbb{C}^{N}$, and $\hat{\mathbf{a}}^{t}$ it the unit vector in the direction of the $\mathbf{a}^{t}$. We next discretize the temporal frequency domain $\Omega$ into $M^{\prime}$ samples and define the measurement vectors with entries

[TABLE]

Rewriting (79) using (80), we obtain the interferometric measurement model as:

[TABLE]

which corresponds to total of $M=M^{\prime}\binom{S}{2}$ cross-correlated measurements666Note that since we have multiple receivers to cross-correlate, the $m=1,\cdots M$ summation of the objective function in (3) now consists of two sums, one over all unique cross-correlations $i\neq j$, and one over the temporal frequency index $m^{\prime}$. In [52], we show that the model in (79) satisfies the sufficient condition of Theorem 4.1, and hence GWF can provide exact image reconstruction for multi-static passive radar.

6.2.2 Simulation Setup and Results

We assume isotropic transmit and receive antennas, and simulate a transmitted signal with $20$MHz bandwidth and $1$GHz center frequency. We assume isotropic transmit and receive antennas, and simulate a transmitted signal with $20$MHz bandwidth and $1$GHz center frequency. We place the center of the scene at the origin of the coordinate system, and generate a phantom, which consists of multiple point and extended targets as depicted in Figure 4. The transmitter is fixed and located at coordinates $\mathbf{a}^{t}=[11.5,11.5,0.5]$ km. We simulate a flat spectrum waveform, and sample the temporal frequency domain into $32$ samples.

i) GWF Reconstruction: To evaluate the performance of GWF for interferometric multistatic radar imaging, we simulate a $300\times 300\mathrm{m}^{2}$ scene, and discretize it by $10$m, which corresponds to $31\times 31$ pixels, hence, $N=961$. We use $16$ receiver antennas that are placed in a circle of radius $10$km around the scene at height of $0.5$km, which corresponds to $120$ unique cross-correlations at each temporal frequency sample. We generate the back-scattered signals at each receiver by the linear measurement map of (77), and correlate linear measurements of the each pairwise combination of receivers to generate interferometric data. In reconstruction, we deploy the approximate measurement vectors in which the transmitter distance is removed as defined in (79) and (80), hence only the transmitter look-direction is used at recovery by GWF.

In Figure 5, we demonstrate the reconstruction obtained by GWF after $10000$ iterations, using the gradually increasing step size heuristic described in Section 6.1. In addition, we plot the relative error of GWF iterates, which displays a geometric rate of convergence as stated in Theorem 4.1. The numerical results on our phantom image indicate GWF has the capacity to form accurate imagery of complex extended scenes.

ii) Comparison to LRMR by Uzawa’s Method: We next perform smaller sized experiments to compare with Uzawa’s method for LRMR by lifting. Namely, we compare GWF to three different formulations that yield Uzawa’s iterations of the form in (14). These are the convex-relaxed trace regularization problem solved by [11], the rank-1 constrained non-convex version with projections defined in (31), and the convex formulation obtained by dropping the rank constraint fully, i.e., a projected gradient descent analogue of [48]777Note that the original proposed method in [48] follows the Douglas-Rachford splitting framework to solve the non-trace regularized convex program in the lifted domain, which has different exact recovery arguments to Uzawa’s method. To avoid the computational burden required to solve a linear system in the lifted domain, we use the projected gradient approach under the identical problem formulation.. Notably, we assess the reconstruction performance with respect to amount of computations, and run the algorithms for identical number of flops. For our problem size, we run GWF for $8000$ iterations, which roughly corresponds to $56$ iterations of Uzawa’s method. As complexity is only known upto an order, to avoid from positive bias towards GWF we run the variations of the Uzawa’s method for $110$ iterations, i.e., double of what the mismatch in computation from lifting the unknown indicates.

The simulated scene corresponds to $144\times 144$m, discritized uniformly by $12$m, to yield a $12\times 12$ unknown, hence $N=144$. We use $12$ receiver antennas that are placed in a circle of radius $10$km around the scene at height of $0.5$km, yielding $66$ unique cross-correlations at each temporal frequency sample. The transmitted waveform bandwidth is set at $10$MHz, around $1.9$GHz central frequency, which corresponds to imaging below the range resolution limit of the system [52]. We provide mean squared error with respect to number of flops in computation over the lifted and the signal domains, and the final images reconstructed by each algorithm in Figures 6 and 7, respectively. Overall, our results indicate that GWF achieves exact reconstruction considerably faster than LRMR by Uzawa’s method when compared at identical number of flops. Furthermore, the convergent behavior of the convex LRMR methods under consideration supports our observation that their exact recovery guarantees are more stringent than those of GWF. Thereby, our numerical results demonstrate its superior performance at identical number of flops to lifting based approaches, and promote GWF as a favorable alternative to state of the art interferometric wave-based imaging methods.

7 Conclusion

In this paper, we present a novel framework for exact interferometric inversion. We approach the interferometric inversion problem from the perspective of phase retrieval techniques. We examine two of the most prominent phase retrieval methods, namely LRMR based PhaseLift, and non-convex optimization based WF, and bridge the theory between the two frameworks. We then generalize WF and formulate the GWF framework for interferometric inversion, and extend the exact recovery guarantees to arbitrary measurement maps with properties that are characterized in the equivalent lifted problem. Thereby, we establish exact recovery conditions for a larger class of problems than that of standard WF.

We identify the sufficient conditions for exact interferometric inversion on the lifted forward model as the restricted isometry property on rank-1, PSD matrices, with a restricted isometry constant of $\delta\leq 0.214$. In developing our theory, we use the special structure of the rank-1, PSD set of matrices to show that the RIP directly implies the regularity condition of WF. Furthermore, we show that the concentration bound of the spectral matrix directly implies the RIP over rank-1, PSD matrices for cross-correlations of the complex Gaussian model. Hence, in generalizing the theory of WF for interferometric inversion in the complex Gaussian case, we demonstrate that the regularity condition becomes redundant. We illustrate that the empirical probability of exact interferometric inversion by GWF require smaller oversampling factors than that of phase retrieval in the Gaussian model. Finally, we demonstrate the applicability of GWF in a deterministic, passive multi-static radar imaging problem using realistic imaging parameters. In conclusion, our paper shows that the computational and theoretical advantages promote GWF as a practical technique in real-world imaging applications.

Acknowledgement

This work was supported by the Air Force Office of Scientific Research (AFOSR) under the agreement FA9550-16-1-0234, Office of Naval Research (ONR) under the agreement N0001418-1-2068 and by the National Science Foundation (NSF) under Grant No ECCS-1809234.

Appendix A Derivations

A.1 Derivation of $\nabla\mathcal{J}$

Recall the objective function $\mathcal{J}$ in (3)

[TABLE]

Letting $e^{m}=\left(\mathbf{L}_{i}^{m}\right)^{H}\bm{\rho}\bm{\rho}^{H}\mathbf{L}_{j}^{m}-d_{ij}^{m}$, we write

[TABLE]

Having $\bm{\rho}$ and $\bm{\rho}^{H}$ independent by properties of Wirtinger derivatives, we compute the partial derivatives in (83) as

[TABLE]

Using the definition of complex gradient provided in (17), we finally get

[TABLE]

A.2 Proof Of Lemma 4.2

Assuming the RIP on $\mathcal{F}$ over rank-1, PSD matrices with RIC-$\delta$, we write

[TABLE]

Equivalently, from the definition of the Frobenius inner product and the adjoint operator $\mathcal{F}^{H}$, we re-express (86) as

[TABLE]

hence for any $\bm{\rho}\in\mathbb{C}^{N}$ we have

[TABLE]

Now consider the definition of the spectral norm, with

[TABLE]

From the Hermitian property of $\bm{\delta}$, we have

[TABLE]

where $\sqrt{\bm{\delta}}:\mathbb{C}^{N\times N}\rightarrow\mathbb{C}^{N\times N}$ with $\sqrt{\bm{\delta}}\sqrt{\bm{\delta}}=\bm{\delta}$, and the inequality follows from matrix Cauchy-Schwartz property, and the fact that Frobenius norm is unaltered by conjugation. Observe that using the Hermitian property of $\bm{\delta}$ via Definition 6.1 in [71], for any $\bm{\rho}\in\mathbb{C}^{N}$, we have $\|\sqrt{\bm{\delta}}(\bm{\rho}\bm{\rho}^{H})\|_{F}^{2}/\|\bm{\rho}\bm{\rho}\|_{F}^{2}\leq\delta$, hence

[TABLE]

Using the definition of the spectral norm, we simply obtain

[TABLE]

for any $\bm{\rho}\in\mathbb{C}^{N}$ if (87) holds. Furthermore since $\delta$ is the smallest constant such that (87) is satisfied, let

[TABLE]

Revisiting (91), since from (93) we know for any $\bm{\rho}$, $\|\bm{\delta}(\bm{\rho}\bm{\rho}^{H})\|\leq\delta\|\bm{\rho}\bm{\rho}^{H}\|_{F}$, the maximal value of $\delta$ is reached at $\mathbf{v}=\bm{\rho}/\|\bm{\rho}\|$, as

[TABLE]

and by definition

[TABLE]

Hence the proof is complete.

Appendix B Lemmas for Theorem 4.1

B.1 Proof of Lemma 5.1

Since $\hat{\bm{\rho}}_{t}=e^{\mathrm{i}\Phi(\mathbf{v}_{0})}\bm{\rho}_{t}$, we have

[TABLE]

Knowing that $\Phi(\mathbf{v}_{0})$ achieves $\text{Re}\langle\hat{\bm{\rho}}_{t},\mathbf{v}_{0}\rangle=|\langle{\bm{\rho}}_{t},\mathbf{v}_{0}\rangle|$, we have

[TABLE]

Since $\hat{\bm{\rho}}_{t}$ and $\mathbf{v}_{0}$ are unit norm, the geometric angle between them can be written as

[TABLE]

Invoking the representation theorem in Hilbert spaces, there exists a unit vector $\mathbf{v}_{0}^{\perp}$ that lies in the plane whose normal is $\mathbf{v}_{0}$ such that

[TABLE]

The inner product is zero only when: 1) $\hat{\bm{\rho}}_{t}=\cos(\theta)\mathbf{v}_{0}+\sin(\theta)\mathbf{v}_{0}^{\perp}$, 2) $\hat{\bm{\rho}}_{t}-(\cos(\theta)\mathbf{v}_{0}+\sin(\theta)\mathbf{v}_{0}^{\perp})$ is perpendicular to $\mathbf{v}_{0}$. The latter case occurs iff $\hat{\bm{\rho}}_{t}-\cos(\theta)\mathbf{v}_{0}=0$, which is true only for $\theta=0$, which indicates the identical solution as the former hence we have the unique representation

[TABLE]

Using the same representation, it is straightforward to see that the unit length element $\hat{\bm{\rho}}_{t}^{\perp}=-\sin(\theta)\mathbf{v}_{0}+\cos(\theta)\mathbf{v}_{0}^{\perp}$ satisfies

[TABLE]

and hence lies, in the plane whose normal is $\hat{\bm{\rho}}_{t}$.

B.2 Proof of Lemma 5.2

Recalling the representation of the spectral estimate in the lifted problem we have

[TABLE]

where $\mathcal{P}_{S}$ is the projection onto the set of symmetric matrices. $\hat{\mathbf{X}}_{PSD}$ is similarly written as

[TABLE]

since the $PSD$-cone is a subset of the set of symmetric matrices $S$, and $\mathcal{P}_{PSD}$ projection consists of projection onto $S$ by $\mathcal{P}_{S}$, followed by suppression of negative eigenvalues. From Lemma 4.2, denoting $\tilde{\bm{\rho}}_{t}=\bm{\rho}_{t}\bm{\rho}_{t}^{H}$ for the solution $\bm{\rho}_{t}$ obeying $\|\bm{\rho}_{t}\|=1$ we have

[TABLE]

And since $\tilde{\bm{\rho}}_{t}$ is Hermitian symmetric and that the spectral norm is unaffected by adjoint operation, we have

[TABLE]

Hence the spectral norm of the lifted error is upper bounded as follows:

[TABLE]

For the maximal eigenvalue $\lambda_{0}$ of $\hat{\mathbf{X}}$, we can write

[TABLE]

On the other hand, we have

[TABLE]

Hence we obtain

[TABLE]

Since the positive semi-definite estimate $\hat{\mathbf{X}}_{PSD}$ only differs from $\hat{\mathbf{X}}$ by the suppression of negative eigenvalues, and since spectral initialization only preserves the leading eigenvalue-eigenvector pair $\lambda_{0}$, $\mathbf{v}_{0}$ and $\lambda_{0}\geq 1-\delta$ where $\delta>0$, we have the identical $\bm{\rho}_{0}=\sqrt{\lambda_{0}}\mathbf{v}_{0}$.

B.3 Proof of Lemma 5.3

We consider the case where the spectral estimate (36) is projected onto the positive semi-definite cone as in LRMR. In this case, the estimate obtained in spectral initialization is projected onto the positive semi-definite cone to obtain $\hat{\mathbf{X}}_{PSD}$, which yields the identical initial point $\bm{\rho}_{0}$ per Lemma 5.2. Since its a symmetric, PSD matrix, we can decompose $\hat{\mathbf{X}}_{PSD}$ such that

[TABLE]

where $\mathbf{S}_{0}$ is a positive semi-definite matrix with its eigenvalues as $\sqrt{\lambda_{\hat{\mathbf{X}}_{PSD}}}$. Since

[TABLE]

where $\mathbf{v}^{H}\hat{\mathbf{X}}_{PSD}\mathbf{v}=\mathbf{v}^{H}\mathbf{S}_{0}\mathbf{S}_{0}\mathbf{v}=\|\mathbf{S}_{0}\mathbf{v}\|^{2}$. Then using the definitions of $\hat{\bm{\rho}}_{t}$ and $\hat{\bm{\rho}}_{t}^{\perp}$ from Lemma 5.1, we have

[TABLE]

and from Pythogorian theorem, since we have orthogonal components, we get

[TABLE]

Consider the expression $f=\sin^{2}(\theta)\|\mathbf{S}_{0}\hat{\bm{\rho}}_{t}\|^{2}-\|\mathbf{S}_{0}\hat{\bm{\rho}}_{t}^{\perp}\|^{2}$. Following the algebra in [43], we have

[TABLE]

which finally yields

[TABLE]

Since $\mathbf{v}_{0}$ is the unit vector that maximizes $\|\mathbf{S}_{0}\mathbf{v}\|^{2}$, with the fact that $\mathbf{v}_{0}^{\perp}$ is also a unit vector, this expression is always non-positive, hence

[TABLE]

[TABLE]

Equivalently, expressing the upper bound with the spectral estimate, we obtain

[TABLE]

Now we consider the spectral estimate $\hat{\mathbf{X}}_{PSD}$. We know that $\hat{\mathbf{X}}_{PSD}$ is obtained by projecting the intermediate estimate $\mathbf{Y}=\mathcal{F}^{H}\mathcal{F}(\tilde{\bm{\rho}}_{t})$ onto the feasible set of PSD matrices as defined in Lemma 5.2. From Lemma 4.2, we know that $\mathcal{F}^{H}\mathcal{F}$ is approximately an identity on the domain of rank-1, PSD matrices, hence we can write the perturbation model

[TABLE]

Since projection onto the PSD cone is non-expansive under the spectral norm, we have the following

[TABLE]

Setting $\hat{\mathbf{X}}_{PSD}=\tilde{\bm{\rho}}_{t}+\tilde{\mathbf{e}}$, the upper bound in (120) can be written as

[TABLE]

Since we have that $\tilde{\bm{\rho}}_{t}=\bm{\rho}_{t}\bm{\rho}_{t}^{H}$, the bound reduces to

[TABLE]

as $\hat{\bm{\rho}}_{t}^{\perp}$ is orthogonal to $\hat{\bm{\rho}}_{t}$ and the global phase component in $\hat{\bm{\rho}}_{t}=e^{j\Phi(v_{0})}\bm{\rho}_{t}$ vanishes in the quadratic form. Moreover, we know that $(\hat{\bm{\rho}}_{t}^{\perp})^{H}\tilde{\mathbf{e}}\hat{\bm{\rho}}_{t}^{\perp}$ is non-negative from positive semi-definitivity of $\hat{\mathbf{X}}_{PSD}$. Hence we further upper bound the numerator by the spectral norm of $\tilde{\mathbf{e}}$, and follow with the series of bounds

[TABLE]

where the last two inequalities follow from the non-expansiveness property of the projection operator onto the PSD cone, and the definition of $\bm{\delta}$ from Lemma 4.2, such that the spectral norm of $\bm{\delta}$ is upper bounded by the RIC-$\delta$ of $\mathcal{F}$. For the term in the denominator, since we established that $\|\tilde{\mathbf{e}}\|_{2}\leq\delta$, we have that $-\delta\leq\hat{\bm{\rho}}_{t}^{H}\tilde{\mathbf{e}}\hat{\bm{\rho}}_{t}\leq\delta$ to finally obtain

[TABLE]

B.4 Proof of Lemma 5.4

B.4.1 The Upper Bound

Noting that $\hat{\bm{\rho}}_{t}$ is the closest solution in $\mathit{P}$ to an estimate $\bm{\rho}$, we define $\bm{\rho}_{e}=\bm{\rho}-\hat{\bm{\rho}}_{t}$. Since $\bm{\rho}\in E(\epsilon)$, $\text{dist}(\bm{\rho},\bm{\rho}_{t})=\|\bm{\rho}_{e}\|\leq\epsilon$, and from reverse triangle inequality we have

[TABLE]

Since $\|\hat{\bm{\rho}}_{t}\|=\|{\bm{\rho}}_{t}\|=1$, we have that $1-\epsilon\leq\|{\bm{\rho}}\|\leq 1+\epsilon$. Having $\hat{\bm{\rho}}_{t}\hat{\bm{\rho}}_{t}^{H}={\bm{\rho}}_{t}{\bm{\rho}}_{t}^{H}$, we let $\tilde{\mathbf{e}}=\bm{\rho}\bm{\rho}^{H}-\hat{\bm{\rho}}_{t}\hat{\bm{\rho}}_{t}^{H}$ denote the error in the lifted problem. Then, writing the lifted error $\tilde{\bm{\rho}}_{e}=\bm{\rho}_{e}\bm{\rho}_{e}^{H}$ as

[TABLE]

finally yields the expression

[TABLE]

for the error in the lifted domain. We can then write the upper bound for $\|\tilde{\mathbf{e}}\|_{F}$ as

[TABLE]

Since all the arguments of the Frobenius norm in the right-hand-side have rank-1, we have $\|\cdot\|_{2}=\|\cdot\|_{F}$. Knowing that $\|\bm{\rho}_{t}\|=1$ we get

[TABLE]

Hence we obtain the upper bound

[TABLE]

B.4.2 The Lower Bound

Expanding the error in the lifted domain, we get

[TABLE]

Since we have the rank-1 lifted signals, the Frobenius norms and the inner product reduce to

[TABLE]

Having $\text{dist}^{2}(\bm{\rho},\bm{\rho}_{t})=\|\bm{\rho}\|^{2}+\|\bm{\rho}_{t}\|^{2}-2|\langle\bm{\rho},\bm{\rho}_{t}\rangle|=\|\bm{\rho}-\hat{\bm{\rho}}_{t}\|\geq 0$, we can lower bound (132) using (133) as

[TABLE]

Knowing that $\text{dist}^{2}(\bm{\rho},\bm{\rho}_{t})\leq\epsilon^{2}$ and the result from the reverse triangle inequality on $\|\bm{\rho}\|$, the terms within the minimization are further lower bounded using

[TABLE]

We then get the bound on the scalar multiplying $\text{dist}^{2}(\bm{\rho},\bm{\rho}_{t})$ as

[TABLE]

which yields the final lower bound as

[TABLE]

B.5 Proof of Lemma 5.5

From the adjoint property of the inner product we have

[TABLE]

Splitting the linear operator $\mathcal{F}^{H}\mathcal{F}$ over the rank-1 inputs, we can use the perturbation model from Lemma 4.2 such that

[TABLE]

Using the representation of $\tilde{\mathbf{e}}=\bm{\rho}\bm{\rho}^{H}-\bm{\rho}_{t}\bm{\rho}_{t}^{H}$ in (128), and the linearity of $\bm{\delta}$ we get

[TABLE]

Notably, the domain of $\bm{\delta}$ is by definition the set of rank-1, PSD matrices. Having $\bm{\rho}_{e}\hat{\bm{\rho}}_{t}^{H}+\hat{\bm{\rho}}_{t}\bm{\rho}_{e}^{H}$, a symmetric matrix of at most rank-2 in the argument of $\bm{\delta}$ on the right-hand-side, we can represent (139) by an eigenvalue decomposition and use the linearity of $\bm{\delta}$ to obtain

[TABLE]

where $\bm{\rho}_{e}\hat{\bm{\rho}}_{t}^{H}+\hat{\bm{\rho}}_{t}\bm{\rho}_{e}^{H}=\sum_{i=1}^{2}\lambda_{i}\mathbf{v}_{i}\mathbf{v}_{i}^{H}$, with $\|\mathbf{v}_{i}\|=1$. Plugging (140) into (138) and applying the triangle inequality, we have

[TABLE]

Furthermore, knowing that $\tilde{\mathbf{e}}$ is symmetric and at most rank-2, let $\tilde{\mathbf{e}}=\sum_{i=1}^{2}\sigma_{i}\mathbf{u}_{i}\mathbf{u}_{i}^{H}$, where $\|\mathbf{u}_{i}\|=1$. Then, using the triangle inequality on the right-hand-side terms in (141), and the outcome of Lemma 4.2 to RIP over rank-1, PSD matrices with RIC-$\delta_{1}$, we have

[TABLE]

Furthermore, since $\tilde{\mathbf{e}}$ is rank-2 by definition, and $\sum_{i=1}^{2}|\sigma_{i}|=\|\tilde{\mathbf{e}}\|_{*}\leq\sqrt{2}\|\tilde{\mathbf{e}}\|_{F}$, we obtain

[TABLE]

where again we use the fact that $\|\bm{\rho}_{e}\hat{\bm{\rho}}_{t}^{H}+\hat{\bm{\rho}}_{t}\bm{\rho}_{e}^{H}\|_{*}=|\lambda_{1}|+|\lambda_{2}|$ by definition, from on which the triangle inequality follows. Invoking the rank-1 property on the terms inside the parenthesis in (143), the Frobenius norm and the nuclear norms can be computed by the spectral norm, since $\|\cdot\|=\|\cdot\|_{F}=\|\cdot\|_{*}$ for a rank-1 argument. Then, having $\|\bm{\rho}_{t}\|=1$, for the right-hand-side we obtain

[TABLE]

Using the bound

[TABLE]

from Lemma 5.4, we finally obtain the upper bound on the perturbation on $\|\bm{\rho}\bm{\rho}^{H}-\bm{\rho}_{t}\bm{\rho}_{t}^{H}\|_{F}^{2}$ as

[TABLE]

Hence, setting $\delta_{2}=\delta_{1}\frac{(2+\epsilon)\sqrt{2}}{\sqrt{(1-\epsilon)(2-\epsilon)}}$, we have the local RIP-2 condition satisfied with RIC-$\delta_{2}$.

B.6 Proof of Lemma 5.6

Having $\tilde{\bm{\rho}}_{t}=\bm{\rho}_{t}\bm{\rho}_{t}^{H}=\hat{\bm{\rho}}_{t}\hat{\bm{\rho}}_{t}^{H}$, we re-write the gradient term in (61) as

[TABLE]

where $\tilde{\mathbf{e}}=\bm{\rho}\bm{\rho}^{H}-\bm{\rho}_{t}\bm{\rho}_{t}^{H}$. Then, simply from the triangle inequality and the definition of projection operator on the set of symmetric matrices, we get

[TABLE]

since the spectral norm is unchanged by Hermitian transpose operation. For the spectral norm of $\bm{\delta}(\tilde{\mathbf{e}})$, we use the representation in (140), and apply the triangle inequality such that

[TABLE]

where again, $\bm{\rho}_{e}\hat{\bm{\rho}}_{t}^{H}+\hat{\bm{\rho}}_{t}\bm{\rho}_{e}^{H}=\sum_{i=1}^{2}\lambda_{i}\mathbf{v}_{i}\mathbf{v}_{i}^{H}$, with $\|\mathbf{v}_{i}\|=1$. We obtain the identical form in (144) such that

[TABLE]

which overall yields

[TABLE]

Using the fact that $\bm{\rho}\in E(\epsilon)$, and $\hat{\bm{\rho}}_{t}=e^{\mathrm{i}\phi(\bm{\rho})}\bm{\rho}_{t}$, we know that $1-\epsilon\leq\|{\bm{\rho}}\|\leq 1+\epsilon$ from triangle and reverse triangle inequalities. Hence,

[TABLE]

Thereby, setting $c=(1+\epsilon)(1+\delta_{1})(2+\epsilon)$ the right-hand-side of the regularity condition is upper bounded as

[TABLE]

and the regularity condition (63) is established if the following condition holds:

[TABLE]

B.7 Proof of Lemma 5.7

Indeed, the form that the regularity condition is nothing but the restricted strong convexity condition. Since $\nabla\mathcal{J}(\bm{\rho}_{t}e^{\mathrm{i}\Phi(\bm{\rho})})=0$ for any $\Phi(\bm{\rho})=\Phi_{0}\in[0,2\pi)$, by definition, one can equivalently write (67) as

[TABLE]

where $\eta=(\frac{1}{\alpha}+\frac{c^{2}}{\beta})$. Reorganizing the terms, we have

[TABLE]

Letting $g(\bm{\rho})=\mathcal{J}(\bm{\rho})-\frac{\eta}{2}\|\bm{\rho}\|^{2}$, we can write (157) as:

[TABLE]

For any $\bm{\rho}$ in the $\epsilon$-ball of $\bm{\rho}_{t}e^{\mathrm{i}\Phi(\bm{\rho})}$, (158) is merely the local convexity condition for $g$ at point $\hat{\bm{\rho}}_{t}=\bm{\rho}_{t}e^{\mathrm{i}\Phi(\bm{\rho})}$. Since the $\epsilon$-ball around $\hat{\bm{\rho}}_{t}$ is a convex set, we can use the equivalent condition

[TABLE]

Plugging in the definition for $g$, we obtain

[TABLE]

Since $\bm{\rho}_{t}e^{\mathrm{i}\Phi(\bm{\rho})}$ is a global minimizer of $\mathcal{J}$, it satisfies the first order optimality condition with $\nabla\mathcal{J}(\bm{\rho}_{t}e^{\mathrm{i}\Phi(\bm{\rho})})=0$, and the minimum it attains is [math]. Hence the condition reduces to

[TABLE]

Appendix C Lemmas for Theorem 4.2

C.1 Proof of Lemma 5.8

For the intermediate stage$,\mathbf{Y}$, of the spectral estimate, we write

[TABLE]

Reorganizing the terms in (162) and taking the expectation, we have

[TABLE]

For fixed $m$, and having $\bm{\rho}_{t}$ is independent of sampling vectors, the $N\times N$ matrix inside the summation has the entries of the form,

[TABLE]

where $k,l$ denote the row and column indexes, respectfully. Since $\mathbf{L}_{i}^{m}$ and $\mathbf{L}_{j}^{m}$ are independent of each other and have i.i.d. entries, the $4^{th}$ moments of Gaussian entries are removed as $\mathbb{E}[(\mathbf{L}_{i}^{m})_{k}(\overline{\mathbf{L}_{i}^{m}})_{n}(\overline{\mathbf{L}_{j}^{m}})_{l}(\mathbf{L}_{j}^{m})_{n^{\prime}}]=\mathbb{E}[(\mathbf{L}_{i}^{m})_{k}(\overline{\mathbf{L}_{i}^{m}})_{n}]\mathbb{E}[(\overline{\mathbf{L}_{j}^{m}})_{l}(\mathbf{L}_{j}^{m})_{n^{\prime}}]$, in which the the expectations are only non-zero for $n=k$, $n^{\prime}=l$, yielding

[TABLE]

where $\mathbf{L}_{i}^{m},\mathbf{L}_{j}^{m}\sim\mathcal{N}(0,\frac{1}{2}\mathbf{I})+\mathrm{i}\mathcal{N}(0,\frac{1}{2}\mathbf{I})$ have unit variance. Hence,

[TABLE]

which is precisely equal to $\mathbb{E}[\mathbf{Y}]=\bm{\rho}_{t}\bm{\rho}_{t}^{H}$.

C.2 Proof of Lemma 5.9

We use the machinery in the proof of the concentration bound of the Hessian in [1].888Corresponds to Lemma 7.4 in the source material. By unitary invariance, we take $\bm{\rho}_{t}=\mathbf{e}_{1}$, where $\mathbf{e}_{1}$ is the first standard basis vector. We want to establish that

[TABLE]

The inequality in (167) is equivalent to

[TABLE]

for any $\mathbf{y}\in\mathbb{C}^{N}$ obeying $\|\mathbf{y}\|=1$. Letting $\mathbf{y}=[\mathbf{y}_{1},\tilde{\mathbf{y}}]$ where $\tilde{\mathbf{y}}\in\mathbb{C}^{N-1}$, and similarly partitioning the sampling vectors as $\mathbf{L}_{i}^{m}=[(\mathbf{L}_{i}^{m})_{1},\tilde{\mathbf{L}}_{i}^{m}]$ we have

[TABLE]

This yields

[TABLE]

Due to the independence $\mathbf{L}_{i}^{m}$ and $\mathbf{L}_{j}^{m}$, and from the fact that they are zero mean, unit variance i.i.d. Gaussian entry vectors, all four terms on the right-hand-side are measures of distance to expected values.

For the second and third terms, we have

[TABLE]

where the inequalities follow from Cauchy-Schwartz and $\ell_{2}\leq\ell_{1}$. We then invoke Hoeffding’s inequality from Proposition 10 in [72]. For any $\delta_{0}$ and $\gamma$, there exists a constant $C(\delta_{0},\gamma)$ such that $M\geq C(\delta_{0},\gamma)\sqrt{N\sum_{m=1}^{M}|(\mathbf{L}_{i}^{m})_{1}|^{4}|(\mathbf{L}_{j}^{m})_{1}|^{2}}$ where

[TABLE]

with probability $1-3e^{-2\gamma N}$. For the final term, we invoke the Bernstein-type inequality of Proposition 16 in [72] per [1], such that for any positive constants $\delta_{0}$, $\gamma$, there exists the constant $C(\delta_{0},\gamma)$ with

[TABLE]

where

[TABLE]

with probability $1-2e^{-2\gamma N}$.

To control the remaining terms, we use Chebyshev’s inequality, per [1]. For any $\epsilon_{0}>0$ there exists a constant $C$ with $M\geq C\cdot N$ such that the following hold

[TABLE]

[TABLE]

with probability at least $1-3N^{-2}$. Moreover, from union bound we have

[TABLE]

with probability at least $1-2N^{-2}$. As in [1], we denote the event that the results from Chebyshev’s inequality hold by $E_{0}$. Then, on event $E_{0}$, combining all the terms, the inequality

[TABLE]

holds with probability at least $1-8e^{-2\gamma N}$. We then follow by the $\epsilon$-net argument of [1] via Lemma 5.4 in [72] to bound the operator norm such that

[TABLE]

where $\mathcal{S}_{\mathbb{C}^{N}}$ is the unit sphere in $\mathbb{C}^{N}$ and $\mathcal{N}$ is an $1/4$-net of $\mathcal{S}_{\mathbb{C}^{N}}$. Then, choosing appropriate $\epsilon_{0}$, $\delta_{0}$ and $\gamma$, and applying the union bound we have

[TABLE]

with probability $1-8e^{-\gamma N}$ for $\delta=2\epsilon_{0}+4\delta_{0}$, and

[TABLE]

From $E_{0}$, we have $M\geq C\cdot N$ which gives $M=\mathcal{O}(N\log N)$, where overall event holds with probability at least $1-8e^{-\gamma N}-5N^{-2}$, hence, the proof is complete.