Solving Complex Quadratic Systems with Full-Rank Random Matrices

TL;DR

This paper demonstrates that complex quadratic systems with full-rank random matrices can be efficiently solved using spectral initialization and Wirtinger flow, extending phase retrieval techniques to more general and practical measurement models.

Contribution

It introduces a novel approach to solve full-rank complex quadratic systems with rotation-invariant sub-Gaussian matrices, broadening applicability beyond rank-1 cases.

Findings

High-probability success of spectral initialization and Wirtinger flow

Effective solution for complex full-rank measurement models

Numerical experiments confirm theoretical results

Abstract

We tackle the problem of recovering a complex signal $\boldsymbol x\in\mathbb{C}^n$ from quadratic measurements of the form $y_i=\boldsymbol x^*\boldsymbol A_i\boldsymbol x$, where $\boldsymbol A_i$ is a full-rank, complex random measurement matrix whose entries are generated from a rotation-invariant sub-Gaussian distribution. We formulate it as the minimization of a nonconvex loss. This problem is related to the well understood phase retrieval problem where the measurement matrix is a rank-1 positive semidefinite matrix. Here we study the general full-rank case which models a number of key applications such as molecular geometry recovery from distance distributions and compound measurements in phaseless diffractive imaging. Most prior works either address the rank-1 case or focus on real measurements. The several papers that address the full-rank complex case adopt the computationally-demanding semidefinite relaxation approach. In this paper we prove that the general class of problems with rotation-invariant sub-Gaussian measurement models can be efficiently solved with high probability via the standard framework comprising a spectral initialization followed by iterative Wirtinger flow updates on a nonconvex loss. Numerical experiments on simulated data corroborate our theoretical analysis.

Figures

Tables

Table 1. TABLE I: Comparison of SVD and the power method.

$n$	Average time [seconds]			Average rel. distance
	SVD	Power		SVD	Power
2000	4.35	0.34		0.522	0.529
5000	54.83	4.35		0.518	0.524
7000	115.48	7.25		0.517	0.521
10000	381.94	37.62		0.517	0.558

Equations

$$y_i={\mathbf{x}}^{\ast }{\mathbf{A}}_i\mathbf{x},i=1,\dots ,m.$$

$$\begin{array}{r}f(\mathbf{z})=\frac{1}{m}\sum _{i=1}^m{\left|{\mathbf{z}}^{\ast }{\mathbf{A}}_i\mathbf{z}-y_i\right|}^2.\end{array}$$

$${\mathbf{z}}^{(t+1)}={\mathbf{z}}^{(t)}-\eta \nabla f(\mathbf{z}),$$

$$\begin{array}{rl}\nabla f(\mathbf{z})={\left(\frac{\partial f}{\partial \mathbf{z}}\right)}^{\ast }=\frac{1}{m}\sum _{i=1}^m& \left({\mathbf{z}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{z}-{\mathbf{x}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{x}\right){\mathbf{A}}_i\mathbf{z}\\ & +\left({\mathbf{z}}^{\ast }{\mathbf{A}}_i\mathbf{z}-{\mathbf{x}}^{\ast }{\mathbf{A}}_i\mathbf{x}\right){\mathbf{A}}_i^{\ast }\mathbf{z}.\end{array}$$

$$\begin{array}{rl}{\mathrm{dist}}^2\left(\mathbf{z},\mathbf{x}\right)& =\underset{\varphi \in (0,2\pi ]}{\min }{\Vert \mathbf{z}-\mathbf{x}{e}^{\mathbf{j}\varphi }\Vert }_2^2\\ & =\Vert \mathbf{z}{\Vert }_2^2+\Vert \mathbf{x}{\Vert }_2^2-2|{\mathbf{z}}^{\ast }\mathbf{x}|,\end{array}$$

$$\mathbf{S}=\frac{1}{m}\sum _{i=1}^m{\overline{y}}_i{\mathbf{A}}_i,$$

$$\begin{array}{r}\mathbb{E}\left[\mathbf{S}\right]=\sum _{rc}{x}_r{\overline{x}}_c\cdot \mathbb{E}\left[\frac{1}{m}\sum _{i=1}^m{\overline{\mathbf{A}}}_{i,rc}\cdot {\mathbf{A}}_i\right].\end{array}$$

$$\mathrm{Pr}\left(A_{i,rc}^{(R)}\in \mathbb{R}\right)$$

$$\mathrm{Pr}\left(A_{i,rc}^{(R)}{A}_{i,kl}^{(R)}\in \mathbb{R}\right)$$

$$\mathbb{E}\left[A_{i,rc}^{(R)}{A}_{i,kl}^{(R)}\right]=\mathbb{E}\left[A_{i,rc}^{(R)}\right]\cdot \mathbb{E}\left[A_{i,kl}^{(R)}\right]=0,$$

$$\mathbb{E}\left[{\overline{\mathbf{A}}}_{i,rc}\cdot {\mathbf{A}}_{i,rc}\right]$$

$$\mathbb{E}\left[{\overline{\mathbf{A}}}_{i,rc}\cdot {\mathbf{A}}_{i,kl}\right]$$

$$\mathbb{E}\left[\mathbf{S}\right]=2\mathbf{x}{\mathbf{x}}^{\ast }.$$

$$\Vert \frac{1}{m}\sum _{i=1}^m{\mathbf{p}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{q}\cdot {\mathbf{A}}_i-2\mathbf{q}{\mathbf{p}}^{\ast }\Vert <\nu ,$$

$${\mathbf{z}}^{(0)}=\Vert \mathbf{x}{\Vert }_2\cdot {\mathbf{v}}_0.$$

$$\mathbb{E}\left[\frac{1}{2m}\sum _{i=1}^m{\overline{y}}_i{y}_i\right]=\mathbb{E}\left[\frac{1}{2}{\mathbf{x}}^{\ast }\mathbf{S}\mathbf{x}\right]=\Vert \mathbf{x}{\Vert }_2^4.$$

$${\mathbf{z}}^{(0)}={\left(\frac{1}{2m}\sum _{i=1}^m{\overline{y}}_i{y}_i\right)}^{\frac{1}{4}}\cdot {\mathbf{v}}_0.$$

$${\mathbf{v}}_0^{(t+1)}=\frac{{\mathbf{S}}^{\ast }\mathbf{S}{\mathbf{v}}_0^{(t)}}{\Vert {\mathbf{S}}^{\ast }\mathbf{S}{\mathbf{v}}_0^{(t)}{\Vert }_2},$$

$${\mathrm{dist}}^2\left({\mathbf{z}}^{(0)},\mathbf{x}\right)\le \frac{51}{24}\delta \Vert \mathbf{x}{\Vert }_2^2,$$

$$\mathcal{P}=\left\{{\mathbf{x}}_0{e}^{\mathbf{j}\varphi }:\varphi \in (0,2\pi ]\right\}.$$

$$E(\rho )=\left\{\mathbf{z}|\mathrm{dist}(\mathbf{z},\mathbf{x})\le \rho \Vert \mathbf{x}{\Vert }_2,\mathbf{x}\in \mathcal{P}\right\},$$

$$\text{Re}\left(⟨\nabla f(\mathbf{z}),\mathbf{z}-\mathbf{x}{e}^{\mathbf{j}{\varphi }_{\min }}⟩\right)\ge \frac{1}{\alpha }{\text{dist}}^2(\mathbf{z},\mathbf{x})+\frac{1}{\beta }\Vert \nabla f(\mathbf{z}){\Vert }_2^2,$$

$${\mathrm{dist}}^2\left({\mathbf{z}}^{(t)},\mathbf{x}\right)\le {\left(1-\frac{2\eta }{\alpha }\right)}^t{\mathrm{dist}}^2\left({\mathbf{z}}^{(0)},\mathbf{x}\right).$$

$$⟨\nabla f(\mathbf{z}),\mathbf{z}-\mathbf{x}⟩\ge \frac{1}{\alpha }{\text{dist}}^2(\mathbf{z},\mathbf{x})+\frac{1}{\beta }\Vert \nabla f(\mathbf{z}){\Vert }_2^2.$$

$$\Vert \frac{1}{m}\sum _{i=1}^m{\mathbf{p}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{q}\cdot {\mathbf{A}}_i-2\mathbf{q}{\mathbf{p}}^{\ast }\Vert <\nu ,$$

$$\mathbb{E}\left[\overline{a_{i,r}}{a}_{i,c}\cdot \overline{a_{i,c}}{a}_{i,r}\right]=\mathbb{E}\left[|a_{i,r}{|}^2|a_{i,c}{|}^2\right]>0.$$

$${\mathrm{dist}}^2\left({\mathbf{z}}^{(t)},\mathbf{x}\right)\le {\left(1-\frac{2\eta }{\alpha }\right)}^t\cdot {\rho }^2\Vert \mathbf{x}{\Vert }_2^2,$$

$${\mathrm{dist}}^2({\mathbf{z}}^{(t+1)},\mathbf{x})\le \left(1-\frac{2\eta }{\alpha }\right){\mathrm{dist}}^2({\mathbf{z}}^{(t)},\mathbf{x}),$$

$$\mathrm{Re}\left(⟨\mathbf{h},\nabla f(\mathbf{z})⟩\right)=\mathrm{Re}(⟨\nabla f(\mathbf{z}),\mathbf{z}-\mathbf{x}{e}^{\mathbf{j}{\varphi }_{\min }}⟩),$$

$$\begin{array}{rl}& \text{Re}\left(⟨\mathbf{h},\nabla f(\mathbf{z})⟩\right)\\ & =\frac{1}{m}\sum _{i=1}^m[2\cdot \text{Re}\left({\mathbf{h}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{h}\cdot {\mathbf{h}}^{\ast }{\mathbf{A}}_i\mathbf{h}\right)+3\cdot \text{Re}\left({\mathbf{h}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{h}\cdot {\mathbf{h}}^{\ast }{\mathbf{A}}_i\mathbf{x}\right)\\ & +3\cdot \text{Re}\left({\mathbf{h}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{x}\cdot {\mathbf{h}}^{\ast }{\mathbf{A}}_i\mathbf{h}\right)+2\cdot \text{Re}\left({\mathbf{h}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{x}\cdot {\mathbf{h}}^{\ast }{\mathbf{A}}_i\mathbf{x}\right)\\ & +\text{Re}\left({\mathbf{h}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{x}\cdot {\mathbf{x}}^{\ast }{\mathbf{A}}_i\mathbf{h}\right)+\text{Re}\left({\mathbf{x}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{h}\cdot {\mathbf{h}}^{\ast }{\mathbf{A}}_i\mathbf{x}\right)].\end{array}$$

$$\begin{array}{rl}& \left|\frac{1}{m}\sum _{i=1}^m\text{Re}({\mathbf{h}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{x}\cdot {\mathbf{h}}^{\ast }{\mathbf{A}}_i\mathbf{x})-2\cdot \text{Re}({\mathbf{h}}^{\ast }\mathbf{x}{\mathbf{h}}^{\ast }\mathbf{x})\right|\\ & \le \left|\frac{1}{m}\sum _{i=1}^m{\mathbf{h}}^{\ast }{\mathbf{A}}_i^{\ast }\mathbf{x}\cdot {\mathbf{h}}^{\ast }{\mathbf{A}}_i\mathbf{x}-2{\mathbf{h}}^{\ast }\mathbf{x}{\mathbf{h}}^{\ast }\mathbf{x}\right|\\ & \le \Vert \mathbf{h}{\Vert }_2^2\Vert \mathbf{x}{\Vert }_2^2\cdot \Vert \frac{1}{m}\sum _{i=1}^m{\hat{\mathbf{h}}}^{\ast }{\mathbf{A}}_i^{\ast }\hat{\mathbf{x}}\cdot {\mathbf{A}}_i-2\hat{\mathbf{x}}{\hat{\mathbf{h}}}^{\ast }\Vert \\ & \le \nu \Vert \mathbf{h}{\Vert }_2^2\Vert \mathbf{x}{\Vert }_2^2,\end{array}$$

Full text

I Introduction

Systems of quadratic equations model many problems in applied science, including phase retrieval [1, 2, 3, 4], the unlabeled distance geometry problem (uDGP) [5, 6], the turnpike and beltway problems [7, 8, 6], unknown view tomography [9, 10, 11], blind channel estimation [12, 13], power flow analysis and power system state estimation [14]. Phase retrieval, in particular, has motivated considerable recent research on quadratic equations. The phaseless measurements are given as $y_{i}=|{\boldsymbol{a}}_{i}^{*}{\boldsymbol{x}}|^{2}={\boldsymbol{x}}^{*}{\boldsymbol{a}}_{i}{\boldsymbol{a}}_{i}^{*}{\boldsymbol{x}}$, with the measurement matrices ${\boldsymbol{a}}_{i}{\boldsymbol{a}}_{i}^{*}$ being rank-1 positive semidefinite matrices. In this paper we study a different measurement model with full-rank measurement matrices. Such measurements arise in a number of the aforementioned applications.

In combinatorial optimization problems such as the uDGP[5, 6] and the nanostructure problem [15, 16], the goal is to reconstruct the relative locations of a set of points from their unlabeled pairwise distances. The distribution $y_{i}$ of the distance $i$ can be formulated as a quadratic form with respect to the point density [6], with the measurement matrix being high-rank. Unknown view tomography aims to reconstruct a 3D density map from a collection of 2D projection images with unknown view angles. When the view angles are assumed to be uniformly distributed on the unit sphere, rotation invariant features can be estimated from 2D projection images and later used for reconstruction [11]. Many of these features can be written as quadratic forms with respect to high-rank measurement matrices.

A variety of combinatorial graph problems can be formulated as quadratic problems with high-rank measurement matrices [17], including the problem of finding the minimum energy spin configuration of atoms located on a grid in quantum physics [18], and the problem of minimizing the number of connections between layers of a circuit board in very-large-scale-integrated (VLSI) circuit design [19]. Furthermore, in an electric transmission network consisting of nodes (buses) and edges (transmission lines), power flow analysis tries to compute the complex voltage at all nodes given measured or specified system variables at selected nodes and edges. This task is then cast as solving a system of quadratic equations with respect to the complex voltages [14] with the measurement matrix being sparse and with rank greater than one. For noisy measurements the task is known as power system state estimation.

These problems can be modeled as systems of quadratic equations where the measurement matrices are not necessarily rank-1 or real. Recovery of a signal from its complex quadratic measurements is naturally formulated as a nonconvex optimization problem, where solving for the globally optimal solution is in general intractable. Recent works on nonconvex quadratic problems such as phase retrieval [20, 21], phase synchronization [22, 23], and low-rank matrix recovery [24] have shown that a globally optimal solution can be recovered from sufficient measurements with high probability when iid Gaussian measurement vectors or matrices are used. Light transport in random media can be modeled by an iid complex Gaussian matrix, and dedicated hardwares like optical processing units (OPU) have been used to produce rapid (20kHz) random projections of high-dimensional data in the million range [25]. In the calibration experiment where the calibration signals are controlled and known, the phases of the calibration measurements can be recovered and used to estimate the complex Gaussian measurement vectors [26]. This application bridges the gap between well-understood random measurement theories and real applications. It further motivates the question of what other types of measurement vectors or matrices also enjoy such favorable properties. In particular, extending the Gaussian measurement model to the sub-Gaussian case has recently attracted considerable interests [27, 28, 29, 30]. Deriving the theory for other ensembles such as Bernoulli [28, 30] has the potential to explain a wider range of applications.

Building upon our earlier work on random Gaussian measurements in [31], we show that our results hold for a slightly larger class of measurement matrices which we term rotation-invariant sub-Gaussian matrices. In this case the signal of interest ${\boldsymbol{x}}\in\mathbb{C}^{n}$, the full-rank sub-Gaussian measurement matrices ${\boldsymbol{A}}_{i}\in\mathbb{C}^{n\times n}$, the measurement $y_{i}\in\mathbb{C}$ are all in the complex domain and we have

[TABLE]

We propose to recover a globally optimal solution via the standard framework comprising a spectral initialization and iterative Wirtinger flow (WF) updates. We prove that when the number of measurements $m$ exceeds the signal length $n$ by some sufficiently large constant $C$, i.e. $m>Cn$, the signal ${\boldsymbol{x}}$ can be recovered up to a global phase shift with high probability.

I-A Prior art

Similar quadratic equation problems have been studied in other contexts. Candès et al. [21] cast the phase retrieval problem as a system of structured quadratic equations and solved it via WF with a linear convergence rate. As this is a non-convex problem, they used a suitably constructed spectral initializer, ${\boldsymbol{z}}^{(0)}$, for the Gaussian measurement model. Spectral initialization for phase retrieval was originally proposed in [20]. It produces ${\boldsymbol{z}}^{(0)}$ which is close to a globally optimal solution with high probability when sufficient measurements are available. The works of [32, 33] subsequently showed that adapting the loss and truncating the measurements adaptively in the initialization and gradient stages could lead to improved performance. Although the spectral initialization was originally developed for random Gaussian measurements, it can also be adapted to work with other types of measurements [34, 6].

Additionally, as shown in the proofs of the WF approach [21], in phase retrieval some of the entries of the measurement matrix are correlated. This makes it impossible to use some of the well-established random matrix theory [35]. In our measurement model, the matrix entries are pairwise uncorrelated and have zero mean, which leads to a different and much simplified proof to establish the convergence guarantees.

Lu and Li [36] studied generalizations of spectral initialization in the real case and focused on the asymptotic behavior of the initializer with respect to the sampling ratio $m/n$ in the high-dimensional limit. Moving beyond the Gaussian measurement model, Ghods et al. [37] proposed a linear spectral estimator for general nonlinear measurement systems. The works of Wang and Xu [38, 39] addressed a generalized phase retrieval problem where ${\boldsymbol{A}}_{i}$ is a Hermitian matrix. They used algebraic methods [40] to find the number of measurements needed for a successful recovery. Here we build upon these results and show that both the initialization and convergence proofs can be derived using the Bernstein-type inequalities for the full-rank rotation-invariant sub-Gaussian measurement model in (1).

Solving systems of quadratic equations is closely related to low-rank matrix recovery—it is equivalent to recovering a rank-1 positive semidefinite (PSD) matrix ${\boldsymbol{X}}={\boldsymbol{x}}{\boldsymbol{x}}^{*}$ with $y_{i}=\langle{\boldsymbol{A}}_{i},{\boldsymbol{X}}\rangle$ [41, 42]. The nonconvex low-rank constraint on ${\boldsymbol{X}}$ can be relaxed to the convex minimum nuclear norm constraint. Alternatively, Carlsson and Gerosa [43] proposed to directly search for a PSD matrix with known rank $k$ and used it to perform phase retrieval from Fourier measurements. The works of [44, 45, 46, 47] focused on establishing sufficient conditions that warrant such a semidefinite convex relaxation. However, recovering the relaxed ${\boldsymbol{X}}$ is computationally demanding even for moderate-size problems. To address this issue, Yurtsever et al. [48] modified the conditional gradient method by performing a small random sketch of ${\boldsymbol{X}}$ that can be used to recover ${\boldsymbol{X}}$ later. On the other hand, the works in [24, 49, 50, 51] approached the problem directly by factoring the unknown matrix as ${\boldsymbol{M}}={\boldsymbol{U}}{\boldsymbol{V}}^{T}$, where ${\boldsymbol{U}}\in\mathbb{R}^{n_{1}\times k},{\boldsymbol{V}}\in\mathbb{R}^{n_{2}\times k}$ and ${\boldsymbol{M}}\in\mathbb{R}^{n_{1}\times n_{2}}$ are all real matrices, and searching for ${\boldsymbol{U}}$ and ${\boldsymbol{V}}$ instead. When ${\boldsymbol{M}}$ is positive semidefinite and the measurement operator satisfies the restricted isometry property [52], Zheng and Lafferty focused on the Gaussian measurement model and proved that ${\boldsymbol{M}}$ could be recovered with guarantees [50]. Tu et al. improved upon [50] and proposed the Procrustes Flow approach for the more general case where $n_{1}\neq n_{2}$ [51]. Our contribution to this line of work lies not just in going from the real case to the complex case in the rotation-invariant sub-Gaussian model (which requires technical interventions at every step), but also in establishing new proofs based on (and including) Lemma 3 that offer a simpler way to obtain the recovery guarantees.

Furthermore, recent works on low-rank matrix recovery showed that the nonconvex problem enjoys a globally optimal landscape [53, 54, 55, 56]. For low-rank square, PSD matrix recovery, Bhojanapalli et al. [53] proved that there are no spurious local minima when using incoherent linear measurements defined by the RIP condition. In this case, a global convergence guarantee can be obtained for stochastic gradient descent from random initialization. Park et al. [54] later generalized it to the non-square case. Via a geometric analysis on the Hessian of the objective function, Ge et al. [55] proposed a way to find directions in which local minima can be “improved”. Recent developments on using various nonconvex optimization approaches to solve low-rank matrix factorization can be found in the review paper by Chi et al. [57].

In this work we complement these results by studying the complex signal recovery problem: we aim to recover a complex signal ${\boldsymbol{x}}\in\mathbb{C}^{n}$ from its complex quadratic measurements up to a global phase shift.

I-B Paper outline

The paper is organized as follows. In Section II-B, we extend the derivations from [20] to the rotation-invariant sub-Gaussian measurement model. We show that the spectral initialization concentrates around a global optimum with high probability and compute the associated concentration bounds. In Section III, we analyze the regularity condition and derive new results for the rotation-invariant sub-Gaussian measurement model. The two results are then combined to give the main theorem of this paper. Computational experiments are presented in Section IV. The proofs of the lemmas are given in the Appendix.

II Problem formulation

II-A Rotation-invariant sub-Gaussian measurement model

For convenience, let ${\boldsymbol{r}}_{i}\in\mathbb{R}^{2n^{2}}$ denote the real and imaginary coefficients of the entries of ${\boldsymbol{A}}_{i}$. The $m$ coefficient vectors ${\boldsymbol{r}}_{i}$ for $1\leq i\leq m$ are independent and identically distributed following a multivariate rotation-invariant sub-Gaussian distribution [58, 35].

• •

The distribution of ${\boldsymbol{r}}_{i}$ does not change under unitary transforms [58]. It follows that the probability density function $p\left({\boldsymbol{r}}_{i}\right)$ depends only on the norm $\|{\boldsymbol{r}}_{i}\|_{2}$. Without loss of generality, we also assume $\mathbb{E}[r_{ik}^{2}]=1,\ \forall\ k=1,\cdots,2n^{2}$.

• •

The coefficient vector ${\boldsymbol{r}}_{i}$ is a sub-Gaussian random vector [35, Definition 5.22] such that the one-dimensional marginals ${\boldsymbol{g}}^{T}{\boldsymbol{r}}_{i}$ are sub-Gaussian random variables for all ${\boldsymbol{g}}\in\mathbb{R}^{2n^{2}}$. In particular, every single entry $r_{ik}$ is also sub-Gaussian by definition.

We shall refer to the model given by (1) as the rotation-invariant sub-Gaussian measurement model, which is a generalization of the Gaussian measurement model from our earlier work [31]. Some examples of matrices in this model are given in Section IV-A.

Instead of solving the system of quadratic equations in (1) directly, we formulate it as a minimization problem. Namely, we minimize the following loss function $f({\boldsymbol{z}})$ to obtain the recovered signal ${\boldsymbol{z}}$:

[TABLE]

Clearly, for any solution ${\boldsymbol{z}}_{0}$ to (1) we have $f({\boldsymbol{z}}_{0})=0$. Although we do not prove it here, one can expect that with sufficiently many “generic” measurements the map from ${\boldsymbol{z}}$ to $[y_{1},\ldots,y_{m}]^{T}$ is injective up to a global phase. Such results have been rigorously proven for the phase retrieval problem [59, 60].

Suppose we are given a “good” initialization point ${\boldsymbol{z}}^{(0)}$ (finding such a point is discussed in Section II-B). The solution is then updated iteratively via Wirtinger flow (WF):

[TABLE]

where $\eta>0$ is some suitable step size, and $\nabla f({\boldsymbol{z}})$ can be computed as

[TABLE]

If ${\boldsymbol{x}}$ is a global minimum of $f({\boldsymbol{z}})$, then ${\boldsymbol{x}}e^{{\boldsymbol{j}}\phi}$ is also a global minimum for any $\phi\in(0,2\pi]$. Consequently, it is standard to define the squared distance between the recovered solution ${\boldsymbol{z}}$ and the true solution ${\boldsymbol{x}}$ as

[TABLE]

where ${\boldsymbol{z}}^{*}{\boldsymbol{x}}=|{\boldsymbol{z}}^{*}{\boldsymbol{x}}|e^{{\boldsymbol{j}}\phi_{{\boldsymbol{z}}^{*}{\boldsymbol{x}}}}$ and the minimum is achieved when $\phi=\phi_{\text{min}}:=-\phi_{{\boldsymbol{z}}^{*}{\boldsymbol{x}}}$.

II-B Spectral initialization

Spectral initialization is widely used in problems with quadratic measurements to obtain an initialization that is close to a global optimum. Similar to [20, 21], we show that the spectral initializer, ${\boldsymbol{z}}^{(0)}$, is close to a global optimizer ${\boldsymbol{x}}$ with high probability and can be used to initialize the WF update in (3). The rationale behind the spectral initialization strategy is that we can get a good estimate ${\boldsymbol{S}}$ of $2{\boldsymbol{x}}{\boldsymbol{x}}^{*}$ using sufficient measurements. The spectral initializer ${\boldsymbol{z}}^{(0)}$ can then be constructed from the leading left or right singular vectors of ${\boldsymbol{S}}$.

Unlike in the phase retrieval problem, which uses the Hermitian matrix $\frac{1}{m}\sum_{i=1}^{m}y_{i}{\boldsymbol{a}}_{i}{\boldsymbol{a}}_{i}^{*}$ as the estimate of ${\boldsymbol{I}}+2{\boldsymbol{x}}{\boldsymbol{x}}^{*}$ (${\boldsymbol{I}}$ is the identity matrix), we propose the following estimation of $2{\boldsymbol{x}}{\boldsymbol{x}}^{*}$ under the rotation-invariant sub-Gaussian measurement model

[TABLE]

where $\overline{y}_{i}={\boldsymbol{x}}^{*}{\boldsymbol{A}}_{i}^{*}{\boldsymbol{x}}$ is the complex conjugate of $y_{i}$. To understand this intuitively, note that the expectation of ${\boldsymbol{S}}$ is

[TABLE]

Let $A_{i,rc}^{(R)}$ and $A_{i,rc}^{(I)}$ denote the real and imaginary coefficients of the $(r,c)$-th entry ${\boldsymbol{A}}_{i,rc}$. Since all the coefficients of the matrix entries in ${\boldsymbol{A}}_{i}$ are generated from a rotation-invariant distribution, any two coefficients are also rotation-invariant, pairwise uncorrelated and “by assumption” have unit variance. Using [58, Proposition 4.1.1], we have:

[TABLE]

where $(r,c)\neq(k,l)$. It is easy to verify that $\mathbb{E}\left[A_{i,rc}^{(R)}\right]=\mathbb{E}\left[A_{i,kl}^{(R)}\right]=0$ and $\mathbb{E}\left[A_{i,rc}^{(R)}A_{i,kl}^{(R)}\right]=0$ so that

[TABLE]

meaning that the coefficients $A_{i,rc}^{(R)}$ and $A_{i,kl}^{(R)}$ are pairwise uncorrelated, and so are “any” two different coefficients of the entries in ${\boldsymbol{A}}_{i}$. Using (10), it is easy to check that

[TABLE]

where $r\neq k$ or $c\neq l$. Hence $\mathbb{E}\left[\frac{1}{m}\sum_{i=1}^{m}\overline{{\boldsymbol{A}}}_{i,rc}{\boldsymbol{A}}_{i}\right]$ is a matrix with the $(r,c)$-th entry equaling to $2$ and the rest of the entries being zeros. Thus

[TABLE]

The following lemma implies that the matrix ${\boldsymbol{S}}$ concentrates around $\mathbb{E}\left[{\boldsymbol{S}}\right]$ in spectral norm with high probability when $m$ is sufficiently large.

Lemma 1**.**

Under the rotation-invariant sub-Gaussian measurement model given by (1), for every $\nu>0$, when the number of measurements satisfies $m>Cn$ for some sufficiently large constant $C:=C(\nu)$, we have for fixed unit vectors ${\boldsymbol{p}},{\boldsymbol{q}}\in\mathbb{C}^{n}$ that

[TABLE]

with probability at least 1-20\exp\big{(}-m\cdot C_{1}(C,\nu)\big{)}, where $C_{1}(C,\nu)>0$ is some constant depending on $C$ and $\nu$.

Note that the statement of Lemma 1 is slightly more general than what we need right now, since it allows ${\boldsymbol{p}}\neq{\boldsymbol{q}}$ (here we set both to $\frac{{\boldsymbol{x}}}{\|{\boldsymbol{x}}\|_{2}}$). This will be useful in the later sections. As derived in Appendix A-A, the concentration proof hinges on the rotation invariance of the measurement matrices ${\boldsymbol{A}}_{i}$ [58]: If we define ${\boldsymbol{B}}={\boldsymbol{R}}{\boldsymbol{A}}$ with ${\boldsymbol{R}}\in\mathbb{C}^{n\times n}$ being a complex unitary matrix, the real and imaginary coefficients of the entries in ${\boldsymbol{B}}$ have the same joint distribution as those of ${\boldsymbol{A}}$.

Lemma 1 implies that ${\boldsymbol{S}}$ is close to $2{\boldsymbol{x}}{\boldsymbol{x}}^{*}$ already for “reasonable” finite values of $m$ (see Fig. 1), and the likelihood that this is not the case decays exponentially with $m$. Let $\{{\boldsymbol{u}}_{0},{\boldsymbol{v}}_{0}\}$ be the “leading” left and right singular vectors of ${\boldsymbol{S}}$, both ${\boldsymbol{u}}_{0}$ and ${\boldsymbol{v}}_{0}$ are highly correlated with $\frac{{\boldsymbol{x}}}{\|{\boldsymbol{x}}\|_{2}}$ with high probability, which is made precise in the proof of the below Lemma 2. We can then use either ${\boldsymbol{u}}_{0}$ or ${\boldsymbol{v}}_{0}$ to construct the spectral initializer ${\boldsymbol{z}}^{(0)}$. Here we shall pick ${\boldsymbol{v}}_{0}$ in the following discussion.

When the norm of ${\boldsymbol{x}}$ is known and fixed, the spectral initializer is

[TABLE] 2. 2.

When the norm of the signal is unknown, we can estimate it from the quadratic measurements. Using (13), we compute the following:

[TABLE]

When $m$ is sufficiently large, we prove that $\frac{1}{2m}\sum_{i=1}^{m}\overline{y}_{i}y_{i}$ is close to its expectation $\|{\boldsymbol{x}}\|_{2}^{4}$ with high probability (see the proof of Lemma 2 in Appendix A-B). Based on this result, we can scale one of the leading singular vectors ${\boldsymbol{v}}_{0}$ of ${\boldsymbol{S}}$ to get our spectral initializer,

[TABLE]

Since ${\boldsymbol{v}}_{0}$ is also the leading eigenvector of ${\boldsymbol{S}}^{*}{\boldsymbol{S}}$ we can use the power iteration to compute it and avoid a full singular value decomposition (SVD) of ${\boldsymbol{S}}$. We initialize it with some random unit-norm vector, ${\boldsymbol{v}}_{0}^{(0)}$, and compute the following power iteration until convergence

[TABLE]

whose computational complexity is $\mathcal{O}(n^{2})$, as opposed to the $\mathcal{O}(n^{3})$ complexity of a full SVD.

One interpretation of spectral initialization is that we are computing an approximate least-squares estimate of the matrix ${\boldsymbol{x}}{\boldsymbol{x}}^{*}$ in the subspace spanned by the measurement matrices. With this guiding principle we can adapt the spectral initialization to other types of measurements. For example, in our work on the uDGP [6] we constructed an orthonormal basis for the matrix subspace spanning the measurement matrices. The measurements can be interpreted as projections of ${\boldsymbol{x}}{\boldsymbol{x}}^{*}$ on to this basis. The least-squares estimate of ${\boldsymbol{x}}{\boldsymbol{x}}^{*}$ and the corresponding spectral initialization can then be easily obtained.

The following lemma states that the distance between the spectral initializer ${\boldsymbol{z}}^{(0)}$ and a global optimizer ${\boldsymbol{x}}$ is small with high probability when $m$ is sufficiently large.

Lemma 2**.**

Under the rotation-invariant sub-Gaussian measurement model given by (1), when the number of complex quadratic measurements satisfies $m>Cn$ for some sufficiently large constant $C$, for every $\delta\in(0,24)$, there exists a global optimizer ${\boldsymbol{x}}$ of the loss in (2) such that the distance between the spectral initializer ${\boldsymbol{z}}^{(0)}$ and ${\boldsymbol{x}}$ obeys

[TABLE]

with probability at least 1-20\exp\big{(}-m\cdot C_{1}(C,\delta)\big{)}, where $C_{1}(C,\delta)>0$ is some constant depending on $C$ and $\delta$.

As a consequence of Lemma 2, no matter how small a $\delta$ we choose, we can guarantee that ${\boldsymbol{z}}^{(0)}$ is $\mathcal{O}(\delta)$-close to a global optimizer ${\boldsymbol{x}}$ with high probability by increasing the number of measurements. As illustrated in Fig. 2, suppose $\rho\|{\boldsymbol{x}}\|_{2}$ is the size of the neighborhood around ${\boldsymbol{x}}$ where a simple local optimization method such as WF can be used to recover ${\boldsymbol{x}}$ with high probability. Such a neighborhood around ${\boldsymbol{x}}$ is referred to as “basin of attraction” in [32]. Our goal is then to balance the trade-off between making $\delta$ small enough so that ${\boldsymbol{z}}^{(0)}$ falls within the basin of attraction and reducing the required number of measurements.

III Convergence analysis

Let ${\boldsymbol{x}}_{0}$ denote a global optimizer, and $\mathcal{P}$ the set of all vectors that differ from ${\boldsymbol{x}}_{0}$ by some phase shift $\phi$:

[TABLE]

In order to determine the neighborhood size $\rho\|{\boldsymbol{x}}\|_{2}$, we study the convergence behavior of the WF iterates in the neighborhood $E(\rho)$ of $\mathcal{P}$, defined as

[TABLE]

where $\mathrm{dist}({\boldsymbol{z}},{\boldsymbol{x}})=\|{\boldsymbol{z}}-{\boldsymbol{x}}e^{{\boldsymbol{j}}\phi_{\min}}\|_{2}$ is computed as in (5). The objective function $f({\boldsymbol{z}})$ is said to satisfy the regularity condition $RC(\alpha,\beta,\rho)$ if the following holds for all ${\boldsymbol{z}}\in E(\rho)$ [21],

[TABLE]

for the choice of constants $\alpha>0$, $\beta>0$, $\rho>0$. The regularity condition $RC(\alpha,\beta,\rho)$ can be derived straightforwardly by demanding that the WF step takes us closer to a global optimizer. More precisely, it ensures that the WF iterate (3) with a step size $\eta\in\left(0,\frac{2}{\beta}\right]$ converges linearly to a global optimizer ${\boldsymbol{x}}$ when the descent is initialized within the neighborhood $E(\rho)$ [21, Lemma 7.10]:

[TABLE]

One of the main challenges in going from the real to the complex case lies in the more complicated definition of the regularity condition and the related convergence analysis. Note that in the real case the regularity condition reads simply

[TABLE]

Following the strategy from the Wirtinger flow paper [21] to lower-bound the left-hand side of (20) would result in complicated derivations involving the computation of the Hessian matrix. We show below how, thanks to our measurement model, we can greatly simplify these derivations using the central Lemma 3.

III-A Establishing the convergence criterion

We now show that there exist choices of parameters $\alpha,\beta,\rho$ such that the objective function $f({\boldsymbol{z}})$ introduced in (2) satisfies the regularity condition $RC(\alpha,\beta,\rho)$ in (20) with high probability, and choose a set of parameter values such that spectral initialization followed by WF succeeds with high probability. The existence of good parameters is shown in three steps according to (20) and (21):

Finding a positive lower bound on $\mathrm{Re}\left(\langle\nabla f({\boldsymbol{z}}),{\boldsymbol{z}}-{\boldsymbol{x}}e^{{\boldsymbol{j}}\phi_{\min}}\rangle\right)$; 2. 2.

Finding an upper bound on $\|\nabla f({\boldsymbol{z}})\|_{2}^{2}$; 3. 3.

Choosing a suitable set of $(\alpha,\beta,\rho)$-values to obtain the regularity condition $RC(\alpha,\beta,\rho)$ in (20).

The main tool in proving these steps is a matrix concentration bound (a high-probability spectral norm bound) similar to the one in Lemma 1. However, Lemma 1 is stated for a particular, fixed choice of the unit vectors ${\boldsymbol{p}},{\boldsymbol{q}}$. Since we want the above bounds which imply the regularity condition to hold for all vectors in $E(\rho)$, it will be useful to strengthen Lemma 1 so that it holds simultaneously for all choices of ${\boldsymbol{p}}$ and ${\boldsymbol{q}}$.

Lemma 3**.**

Under the rotation-invariant sub-Gaussian measurement model given by (1), for every $\nu>0$, when the number of measurements $m$ satisfies $m>Cn$ for some sufficiently large constant $C:=C(\nu)$, we have for all ${\boldsymbol{p}},{\boldsymbol{q}}\in\mathbb{C}^{n}$ satisfying $\|{\boldsymbol{p}}\|_{2}=1$, $\|{\boldsymbol{q}}\|_{2}=1$ that

[TABLE]

with probability at least 1-20\exp\big{(}-m\cdot C_{2}(C,\nu)\big{)}, where $C_{2}(C,\nu)>0$ is some constant depending on $C$ and $\nu$.

In the standard phase retrieval measurement model [21], rows and columns of the measurement matrices are correlated. This is in particular the case for the $(r,c)$-th entry $\overline{a_{i,r}}a_{i,c}$ and the $(c,r)$-th entry $\overline{a_{i,c}}a_{i,r}$. We have

[TABLE]

As a consequence, the distribution of $\overline{a_{i,r}}a_{i,c}\cdot\overline{a_{i,c}}a_{i,r}$ is not centered. This precludes a result parallel to Lemma 3 which in our case lets us establish the regularity condition in a straightforward way.

With Lemma 2 and Lemma 3 in hand, we can now state our main result.

Theorem 1**.**

Under the rotation-invariant sub-Gaussian measurement model given by (1), when the number of complex quadratic measurements $m>Cn$ for some sufficiently large constant $C$,

There exists a choice of $1>\nu>0,\ 1>\rho>0,\ \alpha>0$, and $\beta>0$, such that $RC(\alpha,\beta,\rho)$ holds with probability at least 1-\kappa^{\prime}\cdot\exp\big{(}-m\cdot C_{2}(C,\nu)\big{)}, where $C_{2}(C,\nu)$ is some constant depending $C$ and $\nu$, and $\kappa^{\prime}>0$ is an absolute constant. 2. 2.

Under this choice of parameters, if the step size $\eta$ is chosen so that $0<\eta\leq\frac{2}{\beta}$, the WF iterates (3) initialized at the spectral initializer ${\boldsymbol{z}}^{(0)}$ converge linearly to a global optimizer ${\boldsymbol{x}}$,

[TABLE]

with probability at least 1-\kappa\cdot\exp\big{(}-m\cdot C_{3}(C,\nu,\rho)\big{)}, where $\kappa>0$ is some absolute constant and $C_{3}(C,\nu,\rho)$ depends on $C,\nu$ and $\rho$, but not on $m$..

Proof.

The main task in proving the theorem is to prove Part 1 (the regularity condition). Once we establish that there exists a choice of parameters such that the regularity condition holds in the neighborhood of a global minimizer, Part 2 (linear convergence to a global minimizer) follows simply by noting that $RC(\alpha,\beta,\rho)$ implies

[TABLE]

whenever ${\boldsymbol{z}}^{(t)}\in E(\rho)$.

Part 1: Establishing the regularity condition. We work as follows:

Finding a positive lower bound on $\mathrm{Re}\left(\langle\nabla f({\boldsymbol{z}}),{\boldsymbol{z}}-{\boldsymbol{x}}e^{{\boldsymbol{j}}\phi_{\min}}\rangle\right)$ in the neighborhood of ${\boldsymbol{x}}$

We let ${\boldsymbol{h}}={\boldsymbol{z}}e^{-{\boldsymbol{j}}\phi_{\min}}-{\boldsymbol{x}}$ so that

[TABLE]

and equivalently look for a lower bound on $\mathrm{Re}\left(\langle{\boldsymbol{h}},\nabla f({\boldsymbol{z}})\rangle\right)$ for all ${\boldsymbol{h}},{\boldsymbol{x}}\in\mathbb{C}^{n}$ satisfying $\frac{\|{\boldsymbol{h}}\|_{2}}{\|{\boldsymbol{x}}\|_{2}}\leq\rho$. We proceed by showing that $\mathbb{E}\left[\mathrm{Re}\left(\langle{\boldsymbol{h}},\nabla f({\boldsymbol{z}})\rangle\right)\right]>0$, and that $\mathrm{Re}\left(\langle{\boldsymbol{h}},\nabla f({\boldsymbol{z}})\rangle\right)$ is close to $\mathbb{E}\left[\mathrm{Re}\left(\langle{\boldsymbol{h}},\nabla f({\boldsymbol{z}})\rangle\right)\right]$ with high probability, so that a strictly positive lower bound can be established with high probability when $m$ is sufficiently large. The expression (25) can be expanded as

[TABLE]

We rely on the spectral norm bound in Lemma 3 to lower-bound (26), by bounding each of the six terms in turn. Since the exact same logic applies to all terms, we give details only for the second one, $\frac{1}{m}\sum_{i=1}^{m}\textnormal{Re}({\boldsymbol{h}}^{*}{\boldsymbol{A}}_{i}^{*}{\boldsymbol{x}}\cdot{\boldsymbol{h}}^{*}{\boldsymbol{A}}_{i}{\boldsymbol{x}})$. Let ${\boldsymbol{h}}=\|{\boldsymbol{h}}\|_{2}\cdot\widehat{{\boldsymbol{h}}}$, and ${\boldsymbol{x}}=\|{\boldsymbol{x}}\|_{2}\cdot\widehat{{\boldsymbol{x}}}$, where $\|\widehat{{\boldsymbol{h}}}\|_{2}=1$ and $\|\widehat{{\boldsymbol{x}}}\|_{2}=1$. By Lemma 3 we have

[TABLE]

for all ${\boldsymbol{h}}$ and ${\boldsymbol{x}}$ with probability at least 1-20\exp\big{(}-m\cdot C_{2}(C,\nu)\big{)}. From (5) we see that ${\boldsymbol{h}}^{*}{\boldsymbol{x}}={\boldsymbol{z}}^{*}{\boldsymbol{x}}e^{{\boldsymbol{j}}\phi_{\min}}-\|{\boldsymbol{x}}\|_{2}^{2}=|{\boldsymbol{z}}^{*}{\boldsymbol{x}}|-\|{\boldsymbol{x}}\|_{2}^{2}$ is a real number. We then have the following

[TABLE]

also holds for all ${\boldsymbol{h}}$ and ${\boldsymbol{x}}$ with probability111The probability that (28) holds is no less than the probability that (27) holds. at least 1-20\exp\big{(}-m\cdot C_{2}(C,\nu)\big{)}. Repeating for every term in (26), for all ${\boldsymbol{h}},{\boldsymbol{x}}\in\mathbb{C}^{n}$ satisfying $\frac{\|{\boldsymbol{h}}\|_{2}}{\|{\boldsymbol{x}}\|_{2}}\leq\rho$, we find that

[TABLE]

holds with probability at least 1-\kappa_{1}\exp\big{(}-m\cdot C_{2}(C,\nu)\big{)} where $\kappa_{1}>0$ is some absolute constant.

Finding an upper bound on $\|\nabla f({\boldsymbol{z}})\|_{2}^{2}$ in the neighborhood of ${\boldsymbol{x}}$

We can rewrite $\|\nabla f({\boldsymbol{z}})\|_{2}$ as follows,

[TABLE]

To upper-bound (30) we again rely on the spectral norm bound in Lemma 3. Let ${\boldsymbol{h}}=\|{\boldsymbol{h}}\|_{2}\cdot\widehat{{\boldsymbol{h}}}$, and ${\boldsymbol{x}}=\|{\boldsymbol{x}}\|_{2}\cdot\widehat{{\boldsymbol{x}}}$, where $\|\widehat{{\boldsymbol{h}}}\|_{2}=1$ and $\|\widehat{{\boldsymbol{x}}}\|_{2}=1$. We bound the second term (say) in (30) as follows

[TABLE]

holds with probability at least 1-20\exp\big{(}-m\cdot C_{2}(C,\nu)\big{)}, implying that

[TABLE]

holds with at least the same probability. Repeating for all terms in (30), we get

[TABLE]

with probability at least 1-\kappa_{2}\exp\big{(}-m\cdot C_{2}(C,\nu)\big{)} for all ${\boldsymbol{h}},{\boldsymbol{x}}\in\mathbb{C}^{n}$ satisfying $\frac{\|{\boldsymbol{h}}\|_{2}}{\|{\boldsymbol{x}}\|_{2}}\leq\rho$, where $\kappa_{2}>0$ is some absolute constant.

Choosing suitable $\nu,\rho$ and $\alpha,\beta$

We have that (29) and (33) hold simultaneously with probability at least 1-\kappa^{\prime}\cdot\exp\big{(}-m\cdot C_{2}(C,\nu)\big{)} where $\kappa^{\prime}>0$ is some absolute constant.

In (29), $c_{1}(\nu,\rho)>0$ is a sufficient condition to make $\mathrm{Re}(\langle{\boldsymbol{h}},\nabla f({\boldsymbol{z}})\rangle)>0$ so that we could establish the regularity condition.

If $\nu\geq 1$, we have $c_{1}(\nu,\rho)\leq-\frac{5}{2}\rho-\frac{1}{2}\rho^{2}<0$. Hence $\nu$ needs to be less than $1$. 2. 2.

Given some $\nu\in(0,1)$, there exists a matching $\rho\in(0,1)$ to ensure $c_{1}(\nu,\rho)>0$. The chosen $\nu,\rho$ always lead to $c_{2}(\nu,\rho)>0$ in (33).

In other words, the radius of the convergence neighbourhood $E(\rho)$ could grow as large as $\|{\boldsymbol{x}}\|_{2}$ when there are sufficient measurements.

It remains to show that there exist $\alpha>0$, $\beta>0$ so that

[TABLE]

holds with probability at least 1-\kappa^{\prime}\cdot\exp\big{(}-m\cdot C_{2}(C,\nu)\big{)} for all ${\boldsymbol{h}},{\boldsymbol{x}}\in\mathbb{C}^{n}$ satisfying $\frac{\|{\boldsymbol{h}}\|_{2}}{\|{\boldsymbol{x}}\|_{2}}\leq\rho$. Note that (34) is equivalent to the regularity condition (20).

The parameters $\alpha$ and $\beta$ can be chosen as follows.

From (21) we need $1-\frac{2\eta}{\alpha}\geq 0$. Since $\eta\in(0,\frac{2}{\beta}]$, according to [21, Lemma 7.10], $\alpha$ and $\beta$ should satisfy $\frac{4}{\alpha\beta}\leq 1$. Making the change of variable $\alpha=a/{\|{\boldsymbol{x}}\|_{2}^{2}}$, $\beta=b\cdot c_{2}(\nu,\rho)\|{\boldsymbol{x}}\|_{2}^{2}$, we then require

[TABLE] 2. 2.

From (29) and (33), the regularity condition (20) will hold if

[TABLE]

or equivalently

[TABLE]

There exist many choices of $a$ and $b$ (and thus $\alpha$ and $\beta$) that simultaneously satisfy (35) and (37), and consequently the regularity condition (20). For example, to get the best convergence rate, we can choose $a$, $b$ that maximize $\frac{4}{\alpha\beta}=\frac{4}{c_{2}(\nu,\rho)}\cdot\frac{1}{ab}$ subject to (35) and (37):

• •

If $c_{2}(\nu,\rho)\geq c_{1}(\nu,\rho)^{2}$, we can choose

[TABLE]

• •

If $c_{2}(\nu,\rho)<c_{1}(\nu,\rho)^{2}$, we can choose $a$, $b$ that solve

[TABLE]

Part 2: Linear convergence to a global minimizer with spectral initialization.

After $\nu,\rho,\alpha$, and $\beta$ are chosen in Part 1, it remains to ensure that the spectral initializer ${\boldsymbol{z}}^{(0)}$ falls inside the neighbourhood $E(\rho)$ with high probability. Using Lemma 2, we have that

[TABLE]

holds with probability at least $1-20\exp\left(-m\cdot C_{1}\left(C,\frac{24}{51}\rho^{2}\right)\right)$. Note that choosing $\rho\in(0,1)$ naturally satisfies the constraint imposed on $\frac{24}{51}\rho^{2}$ in Lemma 2, i.e. $\frac{24}{51}\rho^{2}\in\left(0,24\right)$. Combining (34) and (40), we have that the WF update (3) linearly converges to a global minimizer with probability at least 1-\kappa\cdot\exp\big{(}-m\cdot C_{3}(C,\nu,\rho)\big{)}, where $\kappa>0$ is some absolute constant and $C_{3}(C,\nu,\rho)=\min\left\{C_{2}(C,\nu),\ C_{1}(C,\frac{24}{51}\rho^{2})\right\}$ depends on $\{C,\nu,\rho\}$, but not on $m$.

∎

III-B Choosing the $\{\alpha,\beta,\rho\}$-values for the regularity condition

Generally, for some $\nu\in(0,1)$, one begins by choosing a suitable $\rho\in(0,1)$ to obtain a positive lower bound on $\mathrm{Re}\left(\langle\nabla f({\boldsymbol{z}}),{\boldsymbol{z}}-{\boldsymbol{x}}e^{{\boldsymbol{j}}\phi_{\min}}\rangle\right)$ as in (29). The values of $\alpha>0,\ \beta>0$ should be large enough to ensure the regularity condition holds for all ${\boldsymbol{z}}\in E(\rho)$, and satisfy $\frac{4}{\alpha\beta}\leq 1$ so that the iterates ${\boldsymbol{z}}^{(t)}$ come closer to ${\boldsymbol{x}}$ with each iteration.

As derived in the proof of Theorem 1, it is clear that there are many choices for the values of $\{\alpha,\beta,\rho\}$ such that the initializer ${\boldsymbol{z}}^{(0)}\in E(\rho)$ and the objective function $f({\boldsymbol{z}})$ satisfies the regularity condition (20) with high probability. For concreteness, we shall showcase a particular choice of “good” parameter values.

If we choose $\nu=0.01$ and $\rho=0.2$, it is easy to verify that

[TABLE]

where ${\boldsymbol{h}}={\boldsymbol{z}}e^{-{\boldsymbol{j}}\phi_{\min}}-{\boldsymbol{x}}$. We can then choose $\alpha=\frac{2}{3}\cdot\frac{1}{\|{\boldsymbol{x}}\|_{2}^{2}}$ and $\beta=\frac{2}{3}\cdot 120\|{\boldsymbol{x}}\|_{2}^{2}$ to obtain

[TABLE]

When $m$ is sufficiently large, the regularity condition holds for all ${\boldsymbol{h}}$ satisfying $\rho=\frac{\|{\boldsymbol{h}}\|_{2}}{\|{\boldsymbol{x}}\|_{2}}\leq 0.2$ with high probability. We can show that the following update converges linearly to a global optimizer:

[TABLE]

where $0<\eta\leq\frac{2}{\beta}$.

IV Experimental results

We perform numerical experiments to corroborate the theoretical results. Theorem 1 states that the step size is upper-bounded by $\frac{2}{\beta}$ where $\beta$ is one of the regularity condition parameters in (20). In Section III-B $\beta$ is proportional to the squared norm of the signal to recover. Hence, in all experiments the step size is chosen as $\frac{0.1}{\|{\boldsymbol{x}}\|_{2}^{2}}$ where the signal norm, $\|{\boldsymbol{x}}\|_{2}^{2}$, is estimated using (16). The value of $0.1$ was experimentally found to give a suitable balance between convergence speed and reliability.222Code available at https://github.com/swing-research/random_quadratic_equations under the MIT License.

IV-A Rotation-invariant sub-Gaussian distributions

We next illustrate examples of rotation-invariant sub-Gaussian distributions which we use in our experiments. Let

[TABLE]

where $q\in[-1,0]$, and ${\boldsymbol{s}}\sim\mathcal{N}({\boldsymbol{0}},{\boldsymbol{I}}_{d\times d})$, ${\boldsymbol{s}}\neq{\boldsymbol{0}}$. If $q=-1$, the variable ${\boldsymbol{u}}$ is uniformly distributed on the sphere $\|{\boldsymbol{u}}\|_{2}=1$. If $q\in(-1,0)$, the pdf of ${\boldsymbol{u}}$ is

[TABLE]

(see Appendix A-D). If $q=0$, then ${\boldsymbol{u}}={\boldsymbol{s}}$ follows the standard multivariate Gaussian distribution. In the experiments, we further scale ${\boldsymbol{u}}\in\mathbb{R}^{d}$ with the scaling parameter $\gamma>0$: ${\boldsymbol{r}}=\gamma\cdot{\boldsymbol{u}}$ so that the coefficient vector ${\boldsymbol{r}}$ satisfies $\mu^{2}=\text{Var}(r_{i})=\mathbb{E}\left[r_{i}^{2}\right]=1$.

IV-B Closeness of spectral initializer

In this experiment we monitor how the distance between the initialization and the true solution varies with the number of measurements. We fix $n=100$ and try different values of $m$ with $\frac{m}{n}$ uniformly sampled between $1$ and $10$. We run $100$ random trials for each $\frac{m}{n}$ value and calculate the average relative distance between the initialization and a global optimizer. In each trial we generate a random signal ${\boldsymbol{x}}\in\mathbb{C}^{n}$ and $m$ complex random rotation-invariant sub-Gaussian matrices from the same distribution to produce $m$ complex quadratic measurements. We repeat this experiment with multiple distributions by varying the $q$ parameter in Section IV-A.

Distance between complex signals is defined in (5). We define relative distance as $\frac{\mathrm{dist}\left({\boldsymbol{x}},{\boldsymbol{z}}^{(0)}\right)}{\|{\boldsymbol{x}}\|_{2}}$ where ${\boldsymbol{x}}$ is the original signal and ${\boldsymbol{z}}^{(0)}$ is the initialization. In Fig. 3 we can see that the spectral initializer comes closer to a global optimizer as $\frac{m}{n}$ increases. The behavior is not affected by varying $q$.

IV-C Phase transition behavior

In this experiment we evaluate how the proposed approach transits from a failure phase to a success phase as we increase the number of measurements. We fix $n=100$ and try different values of $m$ with $\frac{m}{n}$ sampled uniformly between $1.5$ and $5.5$. We again run 100 random trials for each $\frac{m}{n}$ value and calculate the success rate. Success is declared if the relative distance between the recovered and true signal is less than $10^{-5}$. Again, in each trial a random signal ${\boldsymbol{x}}\in\mathbb{C}^{n}$ is reconstructed and multiple distributions from Section IV-A are used to generate the measurement matrices.

The iterative WF reconstruction is terminated if the relative distance between successive iterations is less than $10^{-6}$ or if 2500 iterations are completed. Here we define relative distance between successive iterates as $\frac{\mathrm{dist}\left({\boldsymbol{z}}^{(t-1)},\ {\boldsymbol{z}}^{(t)}\right)}{\|{\boldsymbol{x}}\|_{2}}$ where ${\boldsymbol{z}}^{(t-1)}$ and ${\boldsymbol{z}}^{(t)}$ are two solutions recovered from successive WF iterates and ${\boldsymbol{x}}$ is the original signal. As the true signal and its norm are unknown during the WF updates, (16) is used to estimate $\|{\boldsymbol{x}}\|_{2}$. Fig. 4 shows that approximately $4n$ measurements are needed to successfully recover the signal for all tested values of $q$.333We note that changing the numerical tolerance in the algorithm stopping criterion can shift the curve.

IV-D Reconstruction of an image

In this experiment we reconstruct an image via its complex quadratic measurements given by (1). For image reconstruction the iterative WF reconstruction is terminated if the distance (5) between successive iterations is less than $10^{-6}$ or if 2500 iterations are completed.444Note that we use distance rather than relative distance. This is a stricter termination criteria when $\|{\boldsymbol{x}}\|>1$. The measurement matrices are from Section IV-A with $q=0$, which corresponds to a complex random Gaussian measurement model. We reconstruct the three color channels of an image of size $n=22\times 15=330$ pixels separately when $\frac{m}{n}=4$. Fig. 5 shows the absolute value of the spectral initialization and the corresponding successful reconstruction.

We define the relative error as $\frac{\||{\boldsymbol{z}}|-{\boldsymbol{x}}\|}{\|{\boldsymbol{x}}\|}$, where $|{\boldsymbol{z}}|$ is the absolute value of the recovered image and ${\boldsymbol{x}}$ is the original image. We further define relative distance as $\frac{\mathrm{dist}\left({\boldsymbol{x}},{\boldsymbol{z}}\right)}{\|{\boldsymbol{x}}\|_{2}}$. When $\frac{m}{n}=4$, the relative error of the spectral initialization is $0.34$. The relative distances between the three channels of the original image and their respective spectral initializations are $0.53$, $0.49$ and $0.51$. The reconstruction relative error is $4.78\times 10^{-7}$. The relative distances between the three channels of the original and their respective reconstructions are $5.45\times 10^{-7}$, $7.73\times 10^{-7}$ and $8.66\times 10^{-7}$.

We also run our algorithm from the beginning to reconstruct the same image for varying number of measurements. For each value of $\frac{m}{n}$ we draw a new set of measurement matrices, calculate a new spectral initialization and use the drawn measurement matrices and initialization for Wirtinger flow updates. Fig. 5 shows a failure case when $\frac{m}{n}=1$ and Fig. 6 shows the relative distances of the recovered images for each channel.

IV-E Comparison between spectral and random initializations

Random initialization has been shown to be a viable alternative to spectral initialization when solving the phase retrieval problem [61]. We next compare the two initialization strategies using the relative distance between the recovered signal and the true signal. Once again we run 100 random trials for each $\frac{m}{n}$ value with $n=100$. We separately analyze the behaviors of the two initializations in the low-oversampling regime, $\frac{m}{n}\in[1,2]$, and the high-oversampling regime, $\frac{m}{n}\in[2,5]$. In each trial we generate a random signal ${\boldsymbol{x}}\in\mathbb{C}^{n}$ and $m$ standard complex multivariate Gaussian measurement matrices to produce $m$ complex quadratic measurements. We compare our proposed spectral initialization against a random initialization that is generated from a standard complex Gaussian distribution. The spectral initialization is computed using 10 power iterations. Both initializations are scaled so that their norms match the estimated signal norm given by (16). The same termination criteria is used as in section IV-C.

The relative distances between the recovered signal and the true signal using different initialization strategies are shown in Fig. 7. We can see that when $1.1\leq\frac{m}{n}\leq 1.7$, the spectral initialization with only 10 power iterations performs better than the random initialization. When $\frac{m}{n}\geq 1.8$, the two initialization strategies perform almost equally well.

IV-F Computational efficiency of spectral initialization

When computing the spectral initialization, the power method is more computationally efficient than a full SVD. In our experiments, we observed that using only 10 power iterations is generally enough to get a good spectral initializer that matches the performance of the exact SVD. In Table I we compare the full SVD and the power method when $\frac{m}{n}=4$ across 5 random trials in terms of the runtimes and the relative distances between the spectral initializers and the true signal. We used up to 24 cores of a system with two 20-core IBM 2.4GHz POWER9 CPUs and up to 115.2GB of RAM. We can see that the power method is preferable since we are only interested in obtaining the leading singular vector.

V Conclusion

We addressed the problem of recovering a signal ${\boldsymbol{x}}\in\mathbb{C}^{n}$ from a system of complex random quadratic equations $y_{i}={\boldsymbol{x}}^{*}{\boldsymbol{A}}_{i}{\boldsymbol{x}}$, for rotation-invariant sub-Gaussian measurement matrices $\{{\boldsymbol{A}}_{i}\}_{i=1}^{m}$. Our analysis complements the existing results on quadratic equations with real measurements and rank-1 positive semidefinite measurement matrices, and extends them to full-rank complex matrices. Since our measurements matrices have uncorrelated entries, the new proofs based on (and including) Lemma 3 can be made much simpler than those for phase retrieval, where the entries of measurement matrices are correlated. Our main result has a standard form: we show that when the number of complex measurements exceeds the length of ${\boldsymbol{x}}$ multiplied by some sufficiently large $C$, then with high probability: 1) the spectral initializer ${\boldsymbol{z}}^{(0)}$ is close to a global optimizer; 2) the WF iterates initialized with ${\boldsymbol{z}}^{(0)}$ converge linearly to a global optimizer. Numerical experiments corroborate the theoretical analysis and show that a global optimum can be successfully recovered when $m$ is sufficiently large.

Recent phase retrieval works showed that a regularized spectral initialization and WF update could improve the robustness and performance of the recovery algorithm [32, 33, 62]. Chen et al. further proved that vanilla gradient descent with random initialization enjoys favorable convergence guarantees in solving the phase retrieval problem [61]. Recent works [53, 54, 55, 56] on the optimization landscape of the low-rank matrix recovery problem give us reason to believe similar optimization landscape could also exist in this case. Our ongoing work involves extending the latest developments in phase retrieval to the aforementioned rotation-invariant sub-Gaussian measurement model. Perhaps more importantly, we hope to adapt our approach to work with quadratic measurements obtained using general high-rank complex matrices that arise in key applications discussed in the introduction.

Appendix A Proofs for rotation-invariant sub-Gaussian measurement model

A-A Proof of Lemma 1

Proof.

Recall that we assume the real and imaginary coefficients of the entries have unit variance, \mathbb{E}\big{[}{A_{i,rc}^{(R)}}^{2}\big{]}=\mathbb{E}\big{[}{A_{i,rc}^{(I)}}^{2}\big{]}=1\,.

For fixed ${\boldsymbol{p}},{\boldsymbol{q}}$, by rotation invariance, we can choose

[TABLE]

where $r_{1},\ r_{2}$ are non-negative real numbers satisfying $r_{1}^{2}+r_{2}^{2}=1$. Let $b_{i}=r_{1}e^{{\boldsymbol{j}}\phi_{1}}\overline{A}_{i,11}+r_{2}e^{{\boldsymbol{j}}\phi_{2}}\overline{A}_{i,21}$, and $\widetilde{{\boldsymbol{A}}}_{i}$ denote the matrix ${\boldsymbol{A}}_{i}$ with the $(1,1)$-th and $(2,1)$-th entries replaced by [math]s. Then

[TABLE]

where ${\boldsymbol{H}}=\frac{1}{m}\sum_{i=1}^{m}b_{i}\widetilde{{\boldsymbol{A}}}_{i}$, $B_{1}=\frac{1}{m}\sum_{i=1}^{m}b_{i}A_{i,11}-2r_{1}e^{{\boldsymbol{j}}\phi_{1}}$, $B_{2}=\frac{1}{m}\sum_{i=1}^{m}b_{i}A_{i,21}-2r_{2}e^{{\boldsymbol{j}}\phi_{2}}$, $\mathrm{Re}(B_{1})$ denotes the real coefficient of $B_{1}$ and $\mathrm{Im}(B_{1})$ denotes the imaginary coefficient of $B_{1}$.

For the first term of (48), we have:

[TABLE]

where

[TABLE]

One can verify that $r_{1}\cos\phi_{1}\left(|A_{i,11}|^{2}-2\right)$ is a centered subexponential random variable (cf. [35, Lemma 5.14]). Similarly, $F_{i}$ is also a centered subexponential random variable. Using the Bernstein-type inequality [35, Proposition 5.16], we have:

[TABLE]

where $c>0$ is some absolute constant and $K_{1}$ and $K_{2}$ are the respective subexponential norms. Combining (49), (51) and (52), we then have

[TABLE]

where $\widehat{C}_{1}(\nu)$ is a constant depending on $\nu$,

[TABLE]

For the second, third and fourth term of (48), we obtain similarly:

[TABLE]

where $\widehat{C}_{2}(\nu),\widehat{C}_{3}(\nu),\widehat{C}_{4}(\nu)$ are some constants depending on $\nu$.

To compute an upper bound on the spectral norm $\|{\boldsymbol{H}}\|$ in (48), we adapt an approach from [35, Theorem 5.39]. The idea is to bound $\left|{\boldsymbol{u}}^{*}{\boldsymbol{H}}{\boldsymbol{v}}\right|$ uniformly for all ${\boldsymbol{u}},{\boldsymbol{v}}\in\mathbb{C}^{n}$ on the unit sphere $\mathscr{S}^{n-1}$. In order to take the union bound over ${\boldsymbol{u}}$s and ${\boldsymbol{v}}$s, the unit sphere $\mathscr{S}^{n-1}$ is first discretized using an $\epsilon$-net $\mathscr{N}_{\epsilon}$ [35, Definition 5.1], for $\epsilon\in[0,1)$. For every fixed pair $({\boldsymbol{u}},{\boldsymbol{v}})$, we establish a high-probability upper bound on $\left|{\boldsymbol{u}}^{*}{\boldsymbol{H}}{\boldsymbol{v}}\right|$, and then take the union bound over $({\boldsymbol{u}},{\boldsymbol{v}})\in\mathscr{N}_{\epsilon}\times\mathscr{N}_{\epsilon}$, taking care of the adjustments so that the result holds over $\mathscr{S}^{n-1}\times\mathscr{S}^{n-1}$.

Approximation

We first bound the error of approximating $\|{\boldsymbol{H}}\|$ using a $({\boldsymbol{u}},{\boldsymbol{v}})$ from $\mathscr{N}_{\epsilon}\times\mathscr{N}_{\epsilon}$. Suppose ${\boldsymbol{u}}_{1},{\boldsymbol{v}}_{1}\in\mathscr{S}^{n-1}$ is chosen such that $\|{\boldsymbol{H}}\|=|\langle{\boldsymbol{H}}{\boldsymbol{v}}_{1},{\boldsymbol{u}}_{1}\rangle|$, and choose ${\boldsymbol{u}}_{2},{\boldsymbol{v}}_{2}\in\mathscr{N}_{\epsilon}$ that approximate ${\boldsymbol{u}}_{1},{\boldsymbol{v}}_{1}$ as $\|{\boldsymbol{u}}_{1}-{\boldsymbol{u}}_{2}\|_{2}\leq\epsilon$, $\|{\boldsymbol{v}}_{1}-{\boldsymbol{v}}_{2}\|_{2}\leq\epsilon$. We get

[TABLE]

It follows that

[TABLE]

Taking the maximum over all ${\boldsymbol{u}}_{2},{\boldsymbol{v}}_{2}$ in the above inequality, we obtain the bound

[TABLE]

We now choose $\epsilon=\frac{1}{4}$. According to [35, Lemma 5.2], there exists a $\frac{1}{4}$-net with cardinality $\left|\mathscr{N}_{\nicefrac{{1}}{{4}}}\right|\leq 9^{n}$. Since we are maximizing over $({\boldsymbol{u}},{\boldsymbol{v}})\in\mathscr{N}_{\nicefrac{{1}}{{4}}}\times\mathscr{N}_{\nicefrac{{1}}{{4}}}$, the total cardinality is bounded as $\left|\mathscr{N}_{\nicefrac{{1}}{{4}}}\right|^{2}\leq 81^{n}$ so that

[TABLE]

Concentration

For a fixed $({\boldsymbol{u}},{\boldsymbol{v}})\in\mathscr{N}_{\nicefrac{{1}}{{4}}}\times\mathscr{N}_{\nicefrac{{1}}{{4}}}$, we have:

[TABLE]

where $G_{i}=\sum_{kl}\overline{u}_{k}v_{l}\cdot\widetilde{{\boldsymbol{A}}}_{i,kl}$. For the first term of (62), $b_{i}^{(R)}$ and $b_{i}^{(I)}$ are linear combinations of the real and imaginary coefficients of $A_{i,11}$ and $A_{i,21}$. On the other hand, $G_{i}^{(R)}$ and $G_{i}^{(I)}$ are linear combinations of the coefficients of the entries in $\widetilde{{\boldsymbol{A}}}_{i}$ that do not contain $A_{i,11}$ and $A_{i,21}$. We can check that $b_{i}^{(R)}G_{i}^{(R)}-b_{i}^{(I)}G_{i}^{(I)}$ is a centered subexponential random variable as before. Using the Bernstein-type inequality[35, Proposition 5.16], we have:

[TABLE]

where $c>0$ is some absolute constant and $K_{3}$ is the subexponential norm of $b_{i}^{(R)}G_{i}^{(R)}-b_{i}^{(I)}G_{i}^{(I)}$. We get a similar result for the second term of (62),

[TABLE]

where $c>0$ is some absolute constant and $K_{4}$ is the subexponential norm of $b_{i}^{(R)}G_{i}^{(I)}+b_{i}^{(I)}G_{i}^{(R)}$.

Combining (62),(63),(64), we have

[TABLE]

where $\widehat{C}_{5}(\nu)$ is a constant depending on $\nu$ and the measurement model,

[TABLE]

Union bound

Taking the union bound over all unit vectors $({\boldsymbol{u}},{\boldsymbol{v}})\in\mathscr{N}_{\nicefrac{{1}}{{4}}}\times\mathscr{N}_{\nicefrac{{1}}{{4}}}$ with cardinality $\left|\mathscr{N}_{\nicefrac{{1}}{{4}}}\right|^{2}\leq 81^{n}$, if $m\geq Cn$,

[TABLE]

Using (60), we have

[TABLE]

If $m\geq Cn$, using (48) and combining all the bounds so far (53), (55), (56), (57), (68), we get

[TABLE]

where $C_{1}(C,\nu)$ depends on $C$ and $\nu$, but not on $m$.

[TABLE]

When $C$ is sufficiently large, we have that $C_{1}(C,\nu)>0$ so that both sides of (69) go to $1$ as $m\to\infty$.

∎

A-B Proof of Lemma 2

Proof.

Let $\{{\boldsymbol{u}}_{0},{\boldsymbol{v}}_{0}\}$ be the leading left and right singular vectors of ${\boldsymbol{S}}$, and $\tau_{0}$ be its largest singular value. Using Lemma 1, the following holds with probability at least 1-20\exp\big{(}-m\cdot C_{1}(C,\delta)\big{)}:

[TABLE]

Hence, on this event, ${\boldsymbol{u}}_{0}^{*}(2{\boldsymbol{x}}{\boldsymbol{x}}^{*}){\boldsymbol{v}}_{0}\geq\tau_{0}-\delta\|{\boldsymbol{x}}\|_{2}^{2}$ and we have with at least the same probability that

[TABLE]

In the following proof the spectral initializer ${\boldsymbol{z}}^{(0)}$ is constructed from ${\boldsymbol{v}}_{0}$. When ${\boldsymbol{z}}^{(0)}$ is constructed from ${\boldsymbol{u}}_{0}$, the proof is similar.

When the signal norm $\|{\boldsymbol{x}}\|_{2}$ is known, we can choose ${\boldsymbol{z}}^{(0)}=\|{\boldsymbol{x}}\|_{2}\cdot{\boldsymbol{v}}_{0}$ as the spectral initializer. Since $|{\boldsymbol{u}}_{0}^{*}{\boldsymbol{x}}|\leq\|{\boldsymbol{x}}\|_{2}$, we have

[TABLE]

Using (5), the squared distance between the spectral initializer ${\boldsymbol{z}}^{(0)}=\|{\boldsymbol{x}}\|_{2}{\boldsymbol{v}}_{0}$ and ${\boldsymbol{x}}$ is then bounded as

[TABLE]

with probability at least 1-20\exp\big{(}-m\cdot C_{1}(C,\delta)\big{)}. 2. 2.

When the signal norm $\|{\boldsymbol{x}}\|_{2}$ is unknown, we estimate it from $R=\frac{1}{2m}\sum_{i=1}^{m}\overline{y}_{i}y_{i}$. By rotation invariance of the sub-Gaussian matrix ${\boldsymbol{A}}_{i}$, we can simply assume ${\boldsymbol{x}}=\|{\boldsymbol{x}}\|_{2}{\boldsymbol{e}}_{1}$ so that

[TABLE]

Using (51) in the proof of Lemma 1 ($\phi_{1}=0$ in this case), we know that on the same event on which (71) holds (of probability \geq 1-20\exp\big{(}-m\cdot C_{1}(C,\delta)\big{)}), it also holds that

[TABLE]

We thus have

[TABLE]

and choose the spectral initializer as ${\boldsymbol{z}}^{(0)}=\sqrt[4]{R}{\boldsymbol{v}}_{0}$. Assuming $\delta\in(0,24)$ and using (73), we have that

[TABLE]

with probability at least 1-20\exp\big{(}-m\cdot C_{1}(C,\delta)\big{)}.

∎

A-C Proof of Lemma 3

Proof.

Let ${\boldsymbol{G}}({\boldsymbol{p}},{\boldsymbol{q}}):=\frac{1}{m}\sum_{i=1}^{m}{\boldsymbol{p}}^{*}{\boldsymbol{A}}_{i}^{*}{\boldsymbol{q}}\cdot{\boldsymbol{A}}_{i}-2{\boldsymbol{q}}{\boldsymbol{p}}^{*}$. We prove that the bound on the spectral norm $\|{\boldsymbol{G}}({\boldsymbol{p}},{\boldsymbol{q}})\|$ in (22) holds with high probability for all unit vectors ${\boldsymbol{p}},{\boldsymbol{q}}$ by combining Lemma 1 with yet another union bound. Let ${\boldsymbol{p}}_{1},{\boldsymbol{q}}_{1}\in\mathscr{S}^{n-1}$ such that

[TABLE]

Let ${\boldsymbol{p}}_{2},{\boldsymbol{q}}_{2}\in\mathscr{N}_{\epsilon}$ further obey ${\boldsymbol{p}}_{1},{\boldsymbol{q}}_{1}$ as $\|{\boldsymbol{p}}_{1}-{\boldsymbol{p}}_{2}\|_{2}\leq\epsilon$, $\|{\boldsymbol{q}}_{1}-{\boldsymbol{q}}_{2}\|_{2}\leq\epsilon$ (they exist by the definition of an $\epsilon$-net). We can write

[TABLE]

so that

[TABLE]

Taking the maximum over all ${\boldsymbol{p}}_{2},{\boldsymbol{q}}_{2}$ in the above inequality, we get

[TABLE]

We again choose $\epsilon=\frac{1}{4}$ so that as in Lemma 1, $|\mathscr{N}_{\nicefrac{{1}}{{4}}}\times\mathscr{N}_{\nicefrac{{1}}{{4}}}|=|\mathscr{N}_{\nicefrac{{1}}{{4}}}|^{2}\leq 81^{n}$ and

[TABLE]

By Lemma 1, if $m>Cn$ for some sufficiently large $C$,

[TABLE]

for fixed unit vectors ${\boldsymbol{p}}_{2},{\boldsymbol{q}}_{2}$. To get a result which holds for all ${\boldsymbol{p}}$ and ${\boldsymbol{q}}$, we take the union bound over $\mathscr{N}_{\nicefrac{{1}}{{4}}}\times\mathscr{N}_{\nicefrac{{1}}{{4}}}$. For $m\geq Cn$,

[TABLE]

Using (83), we have

[TABLE]

where $C_{2}(C,\nu):=C_{1}\left(C,\frac{\nu}{2}\right)-C^{-1}\ln 81>0$ for $C$ sufficiently large.

∎

A-D Example of rotation-invariant distributions

For completeness, we exhibit here one family of rotation-invariant sub-Gaussian distributions. Consider the random variable ${\boldsymbol{u}}=\|{\boldsymbol{s}}\|_{2}^{q}\cdot{\boldsymbol{s}}$, where $q\in(0,1)$, and ${\boldsymbol{s}}\in\mathbb{R}^{d}\sim\mathcal{N}({\boldsymbol{0}},{\boldsymbol{I}})$, ${\boldsymbol{s}}\neq{\boldsymbol{0}}$. We have:

[TABLE]

The entries of the Jacobian matrix ${\boldsymbol{J}}=d{\boldsymbol{s}}/d{\boldsymbol{u}}$ are given as

[TABLE]

The Jacobian matrix ${\boldsymbol{J}}$ is thus

[TABLE]

with the determinant given by

[TABLE]

We obtain the expression for the pdf as

[TABLE]

We can see that $p({\boldsymbol{u}})$ only depends on the norm $\|{\boldsymbol{u}}\|_{2}$, and is thus invariant under unitary transform. For $q=-1$ we obtain a uniform distribution on the sphere, while $q=0$ gives the Gaussian distribution. In all cases, when $q\in[-1,0]$, it is easy to check that the moments of the random variables generated as above are suitably bounded (the tails decay faster than the Gaussian) so that these random variables are sub-Gaussian.