Phase Transitions of Spectral Initialization for High-Dimensional Nonconvex Estimation

I Introduction

where $f(\cdot\,|\,\cdot)$ is a conditional density function modeling the acquisition process. This model arises in many problems in signal processing and statistical learning. Examples include photon-limited imaging , phase retrieval , signal recovery from quantized measurements , and various single-index and generalized linear regression problems .

The standard method for recovering $\boldsymbol{\xi}$ is to use the estimator

In light of these issues, there is strong recent interest in developing and analyzing efficient iterative methods that directly solve nonconvex forms of (2). Examples include the alternating minimization scheme for phase retrieval , the Wirtinger Flow algorithm and its variants , iterative projection methods , and recent schemes for phase retrieval using linear programming . A common ingredient that contributes to the success of these algorithms for nonconvex estimation is that they all use some carefully-designed spectral method as an initialization step, which is then followed by further (iterative) refinement. Beyond the signal estimation problem considered in this paper, related spectral methods have also been successfully applied to initialize algorithms for solving other nonconvex problems such as matrix completion , low-rank matrix recovery , blind deconvolution , sparse coding , and joint alignment from pairwise differences .

In this paper, we present an exact high-dimensional analysis of a widely-used spectral method for estimating $\boldsymbol{\xi}$. The method consists of only two steps: First, construct a data matrix from the sensing vectors and measurements as

The idea of this spectral method can be traced back to the early work of Li , under the name of Principal Hessian Directions for general multi-index models. Similar spectral techniques were also proposed in , for initializing algorithms for matrix completion. In , Netrapalli, Jain, and Sanghavi used this method to address the problem of phase retrieval. Under the assumption that the sensing vectors consist of i.i.d. Gaussian random variables, these authors show that the leading eigenvector $\boldsymbol{x}_{1}$ is aligned with the target vector $\boldsymbol{\xi}$ in direction when there are sufficiently many measurements. More specifically, they show that the squared cosine similarity

which measures the degree of the alignment between the two vectors, approaches $1$ with high probability, when the number of samples $m\geq c_{1}n\log^{3}n$. This sufficient condition on sample complexity was later improved to $m\geq c_{2}n\log n$ in , and further improved to $m\geq c_{3}n$ in with an additional trimming step on the measurements. In these expressions, $c_{1},c_{2},c_{3}$ stand for some unspecified numerical constants.

In this paper, we provide a precise asymptotic characterization of the performance of the spectral method under Gaussian measurements. Our analysis considers general acquisition models under arbitrary conditional distributions $f(y\,|\,\boldsymbol{a}_{i}^{\top}\boldsymbol{\xi})$, of which the phase retrieval problem is a special case. Unlike previous work, which only provides bounds for $\rho(\boldsymbol{\xi},\boldsymbol{x}_{1})$, we derive the exact high-dimensional limit of this value. In particular, we show that, as $n$ and $m$ both tend to infinity with the sampling ratio $\alpha\overset{\text{def}}{=}m/n$ kept fixed, the squared cosine similarity $\rho$ converges in probability to a limit value $\rho(\alpha)$. Explicit formulas are provided for computing $\rho(\alpha)$.

Geometrically, the squared cosine similarity $\rho(\boldsymbol{\xi},\boldsymbol{x}_{1})$ as defined in (4) specifies the angle $\theta$ between $\boldsymbol{\xi}$ and $\boldsymbol{x}_{1}$. The values of $\rho$ vary from to $1$: $\rho=1$ means perfect alignment, i.e., $\theta=0$ or $\pi$; and $\rho=0$ is the opposite case, meaning $\boldsymbol{x}_{1}$ is orthogonal to (i.e. uncorrelated with) $\boldsymbol{\xi}$. That the spectral method can yield an estimate $\boldsymbol{x}_{1}$ with a positive $\rho$ in high dimensional settings is a nontrivial property. To see this, assume that $\boldsymbol{\xi}$ is pointing towards the “north pole” in the unit $(n-1)$-sphere $\mathcal{S}^{n-1}$, as illustrated in Figure 1. If we choose $\boldsymbol{x}_{1}$ uniformly at random from $\mathcal{S}^{n-1}$, then with high probability, the resulting correlation $\sqrt{\rho(\boldsymbol{\xi},\boldsymbol{x}_{1})}$ will be of order $\mathcal{O}(1/\sqrt{n})$. In other words, for large $n$, most of the uniform measure on $\mathcal{S}^{n-1}$ is concentrated within a very thin band of width $\mathcal{O}(1/\sqrt{n})$ near the “equator” of the sphere (see Figure 1).

Our analysis reveals a phase transition phenomenon that occurs at certain critical values of the sampling ratio. In particular, there exist a lower and an upper threshold, denoted by $\alpha_{c,\min}$ and $\alpha_{c,\max}$, respectively, that mark the transitions between two very different phases.

(a) An uncorrelated phase takes place when the sampling ratio $\alpha<\alpha_{c,\min}$. Within this phase, the limiting value $\rho(\alpha)=0$, meaning that the estimate from the spectral method is asymptotically uncorrelated with the target vector $\boldsymbol{\xi}$. In this case, the spectral method is not effective, as its estimate $\boldsymbol{x}_{1}$ is no better than a random guess drawn uniformly from the hypersphere $\mathcal{S}^{n-1}$.

(b) A correlated phase takes place when $\alpha>\alpha_{c,\max}$, with $\alpha_{c,\max}$ being the upper threshold. Within this phase, the limiting value $\rho(\alpha)>0$. Geometrically, the estimate $\boldsymbol{x}_{1}$ (or its negative version $-\boldsymbol{x}_{1}$) will be concentrated on the surface of a right-circular cone (see Figure 1) whose generating lines make an angle \theta=\arccos\big{(}\sqrt{\rho(\alpha)}\big{)} to the target vector $\boldsymbol{\xi}$. Moreover, $\rho(\alpha)$ tends to 1 as $\alpha\rightarrow\infty$.

In many signal estimation models that we have studied so far, the two thresholds coincide, i.e. $\alpha_{c,\min}=\alpha_{c,\max}$, meaning that the phase transition happens at a single critical value of the sampling ratio. However, it is indeed possible that $\alpha_{c,\min}<\alpha_{c,\max}$, in which case a finite number of correlated and uncorrelated phases alternative when $\alpha$ varies within the interval $(\alpha_{c,\min},\alpha_{c,\max})$. A concrete example demonstrating this more complicated situation can be found in Section IV-C.

The above phase transition phenomenon also has implications in terms of the computational complexity of the spectral method. In a correlated phase, there is a nonzero gap between the largest and the second largest eigenvalues of $\boldsymbol{D}_{m}$. As a result, the leading eigenvector $\boldsymbol{x}_{1}$ can be efficiently computed by using power iterations on $\boldsymbol{D}_{m}$. In contrast, within an uncorrelated phase, the gap of the eigenvalues converges to zero, making power iterations inefficient.

The rest of the paper is organized as follows. After precisely laying out the various technical assumptions, we present in Section II the main results of this work, stated as Theorem 1 and Proposition 1. Examples and numerical simulations are also provided there to demonstrate and verify these analytical results. In particular, as a worked example, we derive a universal closed-form expression for the limiting values $\rho(\alpha)$ for all acquisition models that generate one-bit $\left\{0,1\right\}$ measurements. We prove Theorem 1 in Section III. Key to our proof is a deterministic, fixed-point characterization of the squared cosine similarity $\rho(\boldsymbol{\xi},\boldsymbol{x}_{1})$, which is valid for any finite dimension $n$ and for any deterministic sensing vectors $\left\{\boldsymbol{a}_{i}\right\}$. When specialized to Gaussian measurements, this fixed-point characterization allows us to connect our problem to a generalized version of the spiked population model (see, e.g., ) studied in random matrix theory. In Section IV, we look more closely at the phase transition phenomenon predicted by our asymptotic results and prove Proposition 1. Section V concludes the paper with discussions on possible generalizations and improvements of our results as well as their connections to related work in the literature.

II Main Results

In what follows, we first state the basic assumptions under which our results are proved.

The sensing vectors are independent Gaussian random vectors. Specifically, let $(a_{ij})$, for $i,j\geq 1$, be a doubly infinite array of i.i.d. standard normal random variables. Then the $i$th sensing vector $\boldsymbol{a}_{i}=[a_{i1},a_{i2},\ldots,a_{in}]^{\top}$.

$m=m(n)$ with $\alpha_{n}=m(n)/n\rightarrow\alpha>0$ as $n\rightarrow\infty$.

$\mathinner{\!\left\lVert\boldsymbol{\xi}_{n}\right\rVert}=\kappa>0$.

Let $s,y$ and $z$ be three random variables such that

where $f(\cdot\,|\,\cdot)$ is the conditional density function (1) associated with the observation model, and $\mathcal{T}(\cdot)$ is the preprocessing step used in the construction of $\boldsymbol{D}_{m}$ in (3). We shall assume that the probability measure of the random variable $z$ is supported within a finite interval $[0,\tau]$. Throughout the paper, we always take $\tau$ to be the tightest such upper bound.

As $\lambda$ approaches $\tau$ from the right,

The last three assumptions require some explanations. First, we note that assumption (A.4) requires that $z$ should take values within a finite interval on the positive axis. This can be enforced by choosing a suitable function $\mathcal{T}(\cdot)$. For example, in the problem of phase retrieval, the measurement model ($y=s^{2}$) leads to unbounded $\left\{y_{i}\right\}$. We can set

where $t>0$ is some parameter and $\mathds{1}_{\mathinner{\!\left\lvert y\right\rvert}\leq t^{2}}$ denotes the indicator function for the condition $\mathinner{\!\left\lvert y\right\rvert}\leq t^{2}$. This is indeed the trimming strategy proposed in . As shown there, this boundedness condition on the support of $z$ is an essential ingredient in achieving linear sample complexities. The assumption that $z$ be nonnegative is largely made to simplify our analysis, but this restriction can be removed. In a recent work , Mondelli and Montanari extended our results by showing that the same asymptotic predictions presented in this paper still hold under cases where $z$ can take negative values. See also Remark 2 in Section II-B.

The inequality in (A.6) is also a natural requirement. To see this, we note that the data matrix $\boldsymbol{D}_{m}$ in (3) is the sample average of $m$ i.i.d. random rank-one matrices $\left\{\mathcal{T}(y_{i})\boldsymbol{a}_{i}\boldsymbol{a}_{i}^{\top}\right\}_{i\leq m}$. When the number of samples $m$ is large, this sample average should be “close” to the statistical expectation, i.e.,

so that $\boldsymbol{a}_{i}^{\top}\boldsymbol{\xi}=\kappa s_{i}$ and the conditional density of $y_{i}$ given $s_{i}$ is $f(y\,|\,\kappa s_{i})$. Since $s_{i},y_{i}$ and $z_{i}$ are all independent of $\boldsymbol{u}_{i}$,

The above argument provides an intuitive but nonrigorous explanation for why the spectral initialization method would work. The approximation in (9) can be made exact if the signal dimension $n$ is kept fixed and the number of measurement $m$ goes to infinity. However, we consider the case when $m$ and $n$ both tend to infinity, at a constant ratio $\alpha=m/n$ bounded away from and $\infty$. In this regime, the approximation in (9) will not become an equality even if $m\rightarrow\infty$. As we will show, the correlation $\rho(\boldsymbol{\xi}_{n},\boldsymbol{x}_{1}^{n})$ between the target vector $\boldsymbol{\xi}_{n}$ and the sample eigenvector $\boldsymbol{x}_{1}^{n}$ will converge to a deterministic value $\rho(\alpha)$ that depends on the sampling ratio $\alpha$.

A final remark before we present our main results: Since the eigenvector $\boldsymbol{x}_{1}^{n}$ is always normalized, the spectral method cannot provide any information about the norm of $\boldsymbol{\xi}_{n}$. However, in many cases where the sensing vectors are drawn from certain random ensembles, there are simple methods to accurately estimate $\kappa=\mathinner{\!\left\lVert\boldsymbol{\xi}_{n}\right\rVert}$. We provide some discussions on how to do this in Appendix -B.

II-B Main Results: Asymptotic Characterizations

In this section, we summarize the main results of our work on an asymptotic characterization of the spectral method with Gaussian measurements. To state our results, we first need to introduce several helper functions. Let $s,z$ be the random variables defined in (5). We consider two functions

both defined on the open interval $(\tau,\infty)$, where $\tau$ is the bound in assumption (A.4). Within their domains, it is easy to check that both functions are convex. In particular, $\psi_{\alpha}(\lambda)$ achieves its minimum at a unique point denoted by

be a modification of $\psi_{\alpha}(\lambda)$. This new function is again defined for $\lambda\in(\tau,\infty)$.

There is a unique solution, denoted by $\lambda^{\ast}_{\alpha}$, to the equation

where $\psi^{\prime}_{\alpha}(\cdot)$ and $\phi^{\prime}(\cdot)$ denote the derivatives of the two functions.

Let $\lambda_{1}^{\boldsymbol{D}_{m}}\geq\lambda_{2}^{\boldsymbol{D}_{m}}$ be the top two eigenvalues of $\boldsymbol{D}_{m}$.

as $n\rightarrow\infty$. Moreover, $\zeta_{\alpha}(\lambda^{\ast}_{\alpha})\geq\zeta_{\alpha}(\overline{\lambda}_{\alpha})$, with the inequality becoming strict if and only if $\psi^{\prime}_{\alpha}(\lambda^{\ast}_{\alpha})>0$.

The above theorem, whose proof is given in Section III, provides a complete asymptotic characterization of the performance of the spectral method. In particular, the theorem shows that the squared cosine similarity $\rho(\boldsymbol{\xi}_{n},\boldsymbol{x}_{1}^{n})$ converges in probability to a deterministic value in the high-dimensional limit. Moreover, there exists a generic phase transition phenomenon: depending on the sign of the derivative $\psi^{\prime}_{\alpha}(\cdot)$ at $\lambda^{\ast}_{\alpha}$, the limiting value can be either zero (i.e., the uncorrelated phase) or strictly positive (i.e., the correlated phase). The computational complexity of the spectral method is also very different in the two phases. Within the uncorrelated phase, the gap between the top two leading eigenvalues, $\lambda_{1}^{\boldsymbol{D}_{m}}$ and $\lambda_{2}^{\boldsymbol{D}_{m}}$, diminishes to zero, making iterative algorithms such as power iterations increasingly difficult to converge. In contrast, within the correlated phase, the spectral gap converges to a positive value.

The results of this work were first reported in . When this paper was under review, the results given in Theorem 1 were further extended by Mondelli and Montanari in . In particular, these authors extended our asymptotic predictions from the real-valued case to the complex-valued case, and more importantly, they showed that the same predictions still hold under cases where the variable $z$ defined in assumption (A.4) can take negative values. See [33, Lemma 2] for details.

denote the smallest and the largest elements in $\Lambda$, respectively.

Under (A.1) – (A.6), and as $n\rightarrow\infty$,

and $\rho(\alpha)$ is a function with the following parametric representation in terms of a parameter $\lambda$:

for all $\lambda>\lambda_{c,\max}$. Moreover, $\rho(\alpha)\rightarrow 1$ as $\alpha\rightarrow\infty$.

In many of the signal acquisition models we have studied, the set $\Lambda$ contains exactly one element. In this case, $\lambda_{c,\min}=\lambda_{c,\max}$ and hence $\alpha_{c,\min}=\alpha_{c,\max}$. Consequently, the phase transition of the spectral method takes place at a single threshold value $\alpha_{c}$, which separates the uncorrelated phase from the correlated one. However, it is indeed possible to find cases for which $\alpha_{c,\min}<\alpha_{c,\max}$. This leads to a more complicated scenario, where a finite number of correlated and uncorrelated phases can alternatively take place within the interval $(\alpha_{c,\min},\alpha_{c,\max})$. One such example is given in Section IV-C.

II-C Worked-Example: Binary Models

To illustrate the results presented above, we consider here a special case where $z_{i}$ takes only binary values $\left\{0,1\right\}$. This situation naturally appears in problems such as logistic regression and one-bit quantized sensing, where the measurements $y_{i}\in\left\{0,1\right\}$ and we can set $z_{i}=y_{i}$. For cases where the measurements $\left\{y_{i}\right\}$ are not necessarily binary, this type of one-bit model is still relevant whenever the preprocessing function $z=\mathcal{T}(x)$ generates binary outputs. The simplicity of this setting allows us to obtain closed-form expressions for the various quantities in Theorem 1 and Proposition 1.

To proceed, we first explicitly compute the functions $\phi(\lambda)$ and $\psi_{\alpha}(\lambda)$ defined in Section II-B as

and both functions are defined on the interval $\lambda>1$. The minimum of $\psi_{\alpha}(\lambda)$ is achieved as $\overline{\lambda}_{\alpha}=1+\sqrt{\alpha d}$, and thus

Solving equation (17) and using (18), we get

where $\alpha_{c}=\frac{d}{(c-d)^{2}}$. (Note that this result can also be obtained by invoking the parametric characterization of $\rho(\alpha)$ given in Proposition 1.) Finally, the asymptotic predictions (19) for the top two eigenvalues can be computed as

It is interesting to note that the asymptotic characterizations given in (24), (25) and (26) are universal, in the sense that they only depend on the two constants $c$ and $d$ defined in (23) but not on the exact details of the joint probability distributions of $s,y$ and $z$. Thus, for one-bit models, it suffices to compute the constants in (23), which then completely determine the asymptotic performance of the spectral method.

II-D Numerical Simulations

Consider the case where $\left\{y_{i}\right\}$ are binary random variables generated according to the following conditional distribution:

where $\beta$ is some constant. Let $z_{i}=\mathcal{T}(y_{i})=y_{i}$. Since $z_{i}\in\left\{0,1\right\}$, we just need to compute the constants $c$ and $d$ in (23), after which we can use the closed-form expressions (24), (25) and (26) to obtain the asymptotic predictions. In Figure 2(a) we compare the analytical prediction (24) of the squared cosine similarity with results of numerical simulations. In our experiment, we set the signal dimension to $n=4096$. The norm of $\boldsymbol{\xi}_{n}$ is $\kappa=3$, and $\beta=6$. The sample averages and error bars (corresponding to one standard deviation) shown in the figure are calculated over 16 independent trials. We can see that the analytical predictions match numerical results very well. Figure 2(b) shows the top two eigenvalues. When $\alpha<\alpha_{c}$, the two eigenvalues are asymptotically equal, but they start to diverge as $\alpha$ becomes larger than $\alpha_{c}$. To clearly illustrate this phenomenon, we plot in the insert the eigengap $\lambda_{1}-\lambda_{2}$ as a function of $\alpha$.

In the second example, we consider the problem of phase retrieval, where

Here, $\omega_{i}\sim_{\text{i.i.d.}}\mathcal{N}(0,1)$ and $\sigma\geq 0$ is the standard deviation of the noise. In , the authors show that it is important to omit large values of $\left\{y_{i}\right\}$, and they propose to use the scheme in (8) when constructing the data matrix $\boldsymbol{D}_{m}$. A different strategy can be found in , where the authors propose to use

In what follows, we shall refer to (8) and (28) as the trimming algorithm and the subset algorithm, respectively. Figure 3(a) shows the asymptotic performance of these two algorithms and compare them with numerical results ($n=4096$ and 16 independent trials). The performance of the subset algorithm (for which we choose the parameter $t=1.5$) can be characterized by the closed-form formula (24). The trimming algorithm (for which we use $t=3$) is more complicated as $z_{i}$ is no longer binary. We use the parametric characterization in Proposition 1 to obtain its asymptotic performance. Again, our analytical predictions match numerical results. The performance of both algorithms clearly depends on the choice of the thresholding parameter $t$. To show this, we plot in Figure 3(b) the critical phase transition points $\alpha_{c}$ of both algorithms as functions of $t$, at two different noise levels: $\sigma=0$ and $\sigma=2$. This points to the possibility of using our analytical prediction to optimally tune the algorithmic parameters and, more generally, to optimize the functional form of the preprocessing function $\mathcal{T}(\cdot)$. Indeed, the optimal design of $\mathcal{T}(\cdot)$ was obtained in a recent work , which leverages the asymptotic characterizations given here. Interestingly, under a mild technical condition, it is shown that there exists a simple fixed design that is uniformly optimal over all sampling ratios; see [36, Theorem 1].

III Proof of the Main Results

In this section, we prove Theorem 1, which provides an exact characterization of the asymptotic performance of the spectral method for signal estimation.

We first rewrite the data matrix $\boldsymbol{D}_{m}$ in (3) as

where $\boldsymbol{A}=[\boldsymbol{a}_{1},\boldsymbol{a}_{2},\ldots,\boldsymbol{a}_{m}]$ is an $n\times m$ matrix of i.i.d. normal random variables and

is a diagonal matrix with entries $z_{i}=\mathcal{T}(y_{i})$. Our goal boils down to studying the largest eigenvalue of $\boldsymbol{D}_{m}$ and the associated eigenvector $\boldsymbol{x}_{1}^{n}$. To simplify notation, we shall first assume that $\boldsymbol{\xi}_{n}=\kappa\boldsymbol{e}_{1}$, with $\boldsymbol{e}_{1}$ being the first vector in the canonical basis.

The non-null eigenvalues of $\boldsymbol{D}_{m}$ are equal to those of a companion matrix

which bears strong resemblance to a sample covariance matrix. Limiting spectral distributions (LSDs) of sample covariance matrices have been extensively studied in random matrix theory; see for instance and the references given there. As a special case, when $\boldsymbol{Z}$ is the identity matrix, the LSD of $\widetilde{\boldsymbol{D}}_{m}$ is given by the classical Marčenko-Pastur law . Results for more general diagonal matrices $\boldsymbol{Z}$ are also available . However, in these studies, $\boldsymbol{Z}$ and $\boldsymbol{A}$ need to be independent. A challenge in our problem is that $\boldsymbol{Z}$ and $\boldsymbol{A}$ are correlated. To see this, we partition each sensing vector $\boldsymbol{a}_{i}$ into two parts as in (10). We can then write

where $\boldsymbol{s}\overset{\text{def}}{=}[s_{1},s_{2},\ldots,s_{m}]^{\top}$ is an $m$-dimensional Gaussian random vector, and $\boldsymbol{U}$ is an $(n-1)\times m$ matrix consisting of i.i.d. standard normal random variables. Since $\boldsymbol{\xi}_{n}=\kappa\boldsymbol{e}_{1}$, the diagonal elements of $\boldsymbol{Z}$ are independent of $\boldsymbol{U}$ but they do depend on $\boldsymbol{s}$ through $y_{i}\sim f(y\,|\,\kappa s_{i})$. Consequently, we cannot apply existing results on the LSD of sample covariance matrices to our case.

Our proof of Theorem 1 consists of two main ingredients. First, we will show in Proposition 2 that $\lambda_{1}^{\boldsymbol{D}_{m}}$ and $\rho(\boldsymbol{\xi}_{n},\boldsymbol{x}_{1}^{n})$ can be obtained from a fixed-point equation involving a function $L_{m}(\mu)$, to be defined in (36), where $\mu>0$ is an auxiliary variable. The main benefit of introducing the variable $\mu$ and the function $L_{m}(\mu)$ is that, for each $\mu>0$, the above-mentioned correlation between $\boldsymbol{A}$ and $\boldsymbol{Z}$ can be effectively decoupled. This then allows us to obtain the second ingredient of our proof: using results from random matrix theory , we show in Section III-C that $L_{m}(\mu)$, under the assumption of Gaussian sensing vectors, will converge almost surely to a deterministic limit function as the dimension $n\rightarrow\infty$ (see Proposition 4).

III-B A Fixed-Point Characterization

By substituting (31) into (29), we can write $\boldsymbol{D}_{m}$ in a more compact block-partitioned form as

Next, we consider a parametric family of matrices $\left\{\boldsymbol{P}_{m}+\mu\hskip 0.5pt\boldsymbol{q}_{m}\boldsymbol{q}_{m}^{\top}\mathrel{\mathop{\mathchar 58\relax}}{\mu>0}\right\}$, and let $L_{m}(\mu)$ denote their largest eigenvalues, i.e.,

In what follows, we show how to compute $\lambda_{1}^{\boldsymbol{D}_{m}}$ and $\rho(\boldsymbol{\xi}_{n},\boldsymbol{x}_{1}^{n})$ via a fixed-point equation involving $L_{m}(\mu)$. Since we assume that $\boldsymbol{\xi}=\kappa\,\boldsymbol{e}_{1}$ and that the leading eigenvector $\boldsymbol{x}_{1}^{n}$ is normalized, the quantity $\rho(\boldsymbol{\xi}_{n},\boldsymbol{x}_{1}^{n})$ is equal to $(\boldsymbol{e}_{1}^{\top}\boldsymbol{x}_{1}^{n})^{2}$, the squared magnitude of the first element of the eigenvector.

Our discussions below are general and they apply to any block-partitioned matrix in the form

Let $\lambda_{1}^{\boldsymbol{P}}\geq\lambda_{2}^{\boldsymbol{P}}\geq\ldots\lambda_{n-1}^{\boldsymbol{P}}$ be the set of eigenvalues of $\boldsymbol{P}$, and let $\boldsymbol{w}_{1},\boldsymbol{w}_{2},\ldots,\boldsymbol{w}_{n-1}$ be a corresponding set of orthonormal eigenvectors. Consider a function

which has poles on those eigenvalues for which $\boldsymbol{w}_{i}^{\top}\boldsymbol{q}\neq 0$. In what follows, we restrict the domain of $R(\lambda)$ to

Within this open interval, $R(\lambda)$ is a well-defined smooth function. It increases monotonically from $-\infty$ to , and thus it admits a functional inverse, denoted by $R^{-1}(x)$, for all $x<0$. Similar to (36), we define

Let $\boldsymbol{P}$ be a symmetric matrix and $\boldsymbol{q}$ a nonzero vector. Then, for each $\mu>0$,

Moreover, $L(\mu)$ is a nondecreasing convex function with $\lim_{\mu\rightarrow\infty}L(\mu)=\infty$. It is differentiable everywhere on $(0,\infty)$ except at (up to) one point.

Since $\boldsymbol{P}$ is diagonalizable by an orthonormal matrix, we can assume without loss of generality that $\boldsymbol{P}$ is a diagonal matrix. In this case, we can simply write $R(\lambda)=\sum_{i}\frac{q_{i}^{2}}{\lambda_{i}^{\boldsymbol{P}}-\lambda}$, and this function is defined on the open interval $(\max\left\{\lambda_{i}^{\boldsymbol{P}}\mathrel{\mathop{\mathchar 58\relax}}q_{i}\neq 0\right\},\infty)$.

Using the matrix determinant lemma , we can compute the characteristic polynomial of $\boldsymbol{P}+\mu\hskip 0.5pt\boldsymbol{q}\boldsymbol{q}^{\top}$ as

In (40), $\operatorname{adj}(\cdot)$ stands for the adjugate of a matrix. To reach (41), we have used the fact that, for any diagonal matrix $\boldsymbol{A}=\operatorname{diag}\left\{d_{1},d_{2},\ldots,d_{n-1}\right\}$, $\operatorname{adj}(\boldsymbol{A})=\operatorname{diag}\left\{\prod_{j\neq 1}d_{j},\prod_{j\neq 2}d_{j},\ldots,\prod_{j\neq n-1}d_{j}\right\}$.

Partition the set $\left\{1,2,\ldots,n-1\right\}$ into two subsets:

We observe that the characteristic polynomial can be factored into $c(\lambda)=c_{1}(\lambda)c_{2}(\lambda)$, where $c_{2}(\lambda)=\prod_{i\in\mathcal{I}_{2}}(\lambda-\lambda_{i}^{\boldsymbol{P}})$ and

It is possible that the second subset $\mathcal{I}_{2}$ is empty, in which case $c_{2}(\lambda)$ is understood to be equal to $1$, but $\mathcal{I}_{1}$ is never empty, since $\boldsymbol{q}\neq\boldsymbol{0}$. Next, we study the largest root of the polynomial $c_{1}(\lambda)$. For any $\lambda>\max\left\{\lambda_{i}^{\boldsymbol{P}}\mathrel{\mathop{\mathchar 58\relax}}i\in\mathcal{I}_{1}\right\}$, we can write

Recall that $R(\lambda)$ is the function defined in (38) and $R^{-1}(\cdot)$ is its functional inverse. It follows from (43) that $R^{-1}(-1/\mu)$ is the only root of $c_{1}(\lambda)$ in the interval $(\max\left\{\lambda_{i}^{\boldsymbol{P}}\mathrel{\mathop{\mathchar 58\relax}}i\in\mathcal{I}_{1}\right\},\infty)$, and therefore it is also the largest root. Due to the factorization $c(\lambda)=c_{1}(\lambda)c_{2}(\lambda)$, we have

Finally, since $R^{-1}(-1/\mu)>\max\left\{\lambda_{i}\mathrel{\mathop{\mathchar 58\relax}}i\in\mathcal{I}_{1}\right\}$, we reach the formula in (39).

By construction, $R^{-1}(-1/\mu)$ is strictly increasing and $\lim_{\mu\rightarrow\infty}R^{-1}(-1/\mu)=\infty$. It is also differentiable everywhere on $(0,\infty)$. It follows that $L(\mu)$ is nondecreasing with $L(\infty)=\infty$, and that the function is differentiable everywhere except for at most one point $\mu_{0}$, which, if it exists, must satisfy the identity $R^{-1}(-1/\mu_{0})=\lambda_{1}^{\boldsymbol{P}}$. Finally, the convexity of $L(\mu)$ follows from the fact that it is the maximum of a set of linear functions, as $L(\mu)=\lambda_{1}(\boldsymbol{P}+\mu\hskip 0.5pt\boldsymbol{q}\boldsymbol{q}^{\top})=\max_{\boldsymbol{x}\mathrel{\mathop{\mathchar 58\relax}}\,\mathinner{\!\left\lVert\boldsymbol{x}\right\rVert}=1}\boldsymbol{x}^{\top}(\boldsymbol{P}+\mu\hskip 0.5pt\boldsymbol{q}\boldsymbol{q}^{\top})\boldsymbol{x}$. ∎

Given a block-partitioned matrix $\boldsymbol{D}$, the following proposition shows that its leading eigenvalue $\lambda_{1}^{\boldsymbol{D}}$ and the squared cosine similarity $(\boldsymbol{e}_{1}^{\top}\boldsymbol{x}_{1})^{2}$ can be obtained from the function $L(\mu)$.

Let $\mu^{\ast}>0$ be the unique solution to the fixed-point equation

Then, $\lambda_{1}^{\boldsymbol{D}}=L(\mu^{\ast})$ and

where $\partial_{-}L(\mu)$ and $\partial_{+}L(\mu)$ denote the left and right derivatives of $L(\mu)$, respectively. In particular, if $L(\mu)$ is differentiable at $\mu^{\ast}$, then

We prove this result in Appendix -C. Note that (44) is equivalent to

Since $L(\mu)$ is nondecreasing with $L(\infty)=\infty$ whereas $a+1/\mu$ decreases monotonically from $\infty$ to , the equation (47), and thus (44), always admits one and only one solution. Moreover, by Lemma 1, $L(\mu)$ is a convex function, and therefore its left and right derivatives always exist.

The characterization given in Proposition 2 is valid for any block-partitioned matrix in the form of (37). When applied to the specific case of our data matrix in (32), with its components $a_{m}$, $\boldsymbol{P}_{m}$ and $\boldsymbol{q}_{m}$ defined as in (33), (34) and (35), this result provides a very general deterministic characterization of the performance of the spectral method that is valid for any finite dimension $n$ and for any sensing vectors.

Next, we specialize to the case of i.i.d. Gaussian sensing vectors and show that $L_{m}(\mu)$ converges almost surely to a deterministic function as $m,n\rightarrow\infty$. To that end, we note that $L_{m}(\mu)$ is the leading eigenvalue of

is a rank-one perturbation of the diagonal matrix $\boldsymbol{Z}$ given in (30). Since $\boldsymbol{U}$ and $\boldsymbol{M}_{m}$ are independent, we first study the spectrum of $\boldsymbol{M}_{m}$.

Let $\lambda_{1}^{\boldsymbol{M}_{m}}\geq\lambda_{2}^{\boldsymbol{M}_{m}}\geq\ldots\geq\lambda_{m}^{\boldsymbol{M}_{m}}$ be the set of eigenvalues of $\boldsymbol{M}_{m}$ in descending order. Let

be the empirical spectral measure of the last $m-1$ eigenvalues.

Fix $\mu>0$. As $m,n\rightarrow\infty$, the empirical spectral measure $f^{\boldsymbol{M}_{m}}(\lambda)$ converges almost surely to the probability law of the random variable $z$. Meanwhile,

where $Q^{-1}(\cdot)$ is the functional inverse of the function

The domain of $Q(\lambda)$ is the open interval $(\tau,\infty)$, with $\tau$ being the upper bound of the support of the probability law of $z$.

By construction, $Q(\lambda)$ is a continuous and strictly decreasing function with $Q(\infty)=0$. Assumption (A.5) further guarantees that $\lim_{\lambda\rightarrow\tau^{+}}Q(\lambda)=\infty$. Thus, $Q(\lambda)$ admits a functional inverse and that $Q^{-1}(1/\mu)$ is well-defined for all $\mu>0$.

According to assumption (A.4) stated in Section II-A, the law of $z$ is supported within the interval $[0,\tau]$. The above proposition, whose proof can be found in Appendix -D, shows that the spectrum of $\boldsymbol{M}_{m}$ consists of two parts: a “bulk spectrum” of $m-1$ eigenvalues supported within $[0,\tau]$ and a single spiked eigenvalue $\lambda_{1}^{\boldsymbol{M}_{m}}$ well separated from the bulk. This setting is a generalization of the classical spiked population model . Adapting the results given in (see also for related results under more general settings), we thus reach the second important ingredient of our proof of Theorem 1, characterizing the asymptotic limit of $L_{m}(\mu)$.

where $\zeta_{\alpha}(\cdot)$ is the function defined in (16) and $Q^{-1}(1/\mu)$ is the limit value in (50).

Recall from (48) that $L_{m}(\mu)$ is the leading eigenvalue of $\tfrac{1}{m}\boldsymbol{U}\boldsymbol{M}_{m}\boldsymbol{U}^{\top}$. Since $\boldsymbol{U}$ and $\boldsymbol{M}_{m}$ are independent, and since $\boldsymbol{U}$ is a Gaussian random matrix with a rotationally invariant distribution, we can equivalently study the leading eigenvalue of the following matrix

Proposition 3 shows that $\left\{\lambda_{i}^{\boldsymbol{M}}\mathrel{\mathop{\mathchar 58\relax}}i\geq 2\right\}$ form a bulk spectrum, which converges to the law of $z$ as $m\rightarrow\infty$, whereas $\lambda_{1}^{\boldsymbol{M}}$ converges to a “spike” $\lambda_{\mu}=Q^{-1}(1/\mu)>\tau$, which is separated from the bulk.

The asymptotic limits of extreme sample eigenvalues of matrices in the form of (53) have been studied in . In our proof, we use the asymptotic characterization given in . Key to this asymptotic analysis is the function $\psi_{\alpha}(\lambda)$ definedWe have adapted the original definition of $\psi_{\alpha}(\lambda)$ in [32, eq. (3.2)] because our matrix in (53) has a slightly different scaling from the one considered in . in (14). The asymptotic behaviors of the leading sample eigenvalue turn out to depend on the sign of $\psi^{\prime}_{\alpha}(\lambda)$ at the point $\lambda_{\mu}$:

In particular, applying [32, Theorem 4.1], we have

The case when $\psi^{\prime}_{\alpha}(\lambda_{\mu})\leq 0$ is covered in [32, Theorem 4.2]. Adapting that result to our specific setting, we have

III-D Proof of Theorem 1

We are now ready to prove our asymptotic characterizations given in Theorem 1. Since the sensing vectors $\boldsymbol{a}_{i}$ are drawn from the rotationally invariant multivariate normal distribution, the quantity $\rho(\boldsymbol{\xi}_{n},\boldsymbol{x}_{1}^{n})$ for a general vector $\boldsymbol{\xi}_{n}$ (with $\mathinner{\!\left\lVert\boldsymbol{\xi}_{n}\right\rVert}=\kappa$) and $\rho(\kappa\boldsymbol{e}_{1},\boldsymbol{x}_{1}^{n})$ for the special case $\boldsymbol{\xi}_{n}=\kappa\boldsymbol{e}_{1}$ have exactly the same probability distribution. In what follows, we will carry out the proof by assuming that the target vector $\boldsymbol{\xi}_{n}=\kappa\boldsymbol{e}_{1}$. By showing that $\rho(\kappa\boldsymbol{e}_{1},\boldsymbol{x}_{1}^{n})$ converges to the right-hand side of (18) almost surely, the convergence to the same limit in probability for a general $\boldsymbol{\xi}_{n}$ then follows as an immediate consequence.

To start, we use the deterministic characterization given in Proposition 2. For each $m\geq 1$, let $\mu_{m}$ be the unique fixed-point of (44). Equivalently, $\mu_{m}$ satisfies the identity

To determine the asymptotic behavior of the leading eigenvector $\boldsymbol{x}_{1}^{n}$, we use the characterization given in (45). Since $\left\{L_{m}(\mu)\right\}$ are convex functions, we apply Lemma 4 in Appendix -E. In particular, if $\zeta_{\alpha}(Q^{-1}(1/\mu))$ is differentiable at $\mu=\mu^{\ast}$, that lemma gives us

Substituting these limits into (45), we get

To simplify the above expression, we introduce a change of variable, writing $\lambda=Q^{-1}(1/\mu)$. In particular, $\lambda^{\ast}=Q^{-1}(1/\mu^{\ast})$. Using the characterization (57) and recalling the definition of $Q(\lambda)$ in (51), we get

where $\phi(\cdot)$ is defined in (13). By their constructions, it is easily checked that $\zeta_{\alpha}(\lambda)$ is a nondecreasing continuous function on $(\tau,\infty)$ whereas $\phi(\lambda)$ is a strictly decreasing continuous function. Moreover, by assumption (A.5), $\lim_{\lambda\rightarrow\tau^{+}}\phi(\lambda)=\infty$. Thus, the existence of $\lambda^{\ast}$ satisfying (59) and its uniqueness are guaranteed. Substituting $\lambda^{\ast}=Q^{-1}(1/\mu^{\ast})$ into (58) gives us

where we have also used the fact that $Q^{\prime}(\lambda)=\phi^{\prime}(\lambda)$. To reach the characterization (18) given in the theorem, we just need to note that, by its definition in (16), $\zeta^{\prime}_{\alpha}(\lambda)=\psi^{\prime}_{\alpha}(\lambda)$ if $\psi^{\prime}_{\alpha}(\lambda)>0$ and $\zeta^{\prime}_{\alpha}(\lambda)=0$ if $\psi^{\prime}_{\alpha}(\lambda)<0$.

Next, we characterize the first two eigenvalues $\lambda_{1}^{\boldsymbol{D}_{m}}$ and $\lambda_{2}^{\boldsymbol{D}_{m}}$. By Proposition 2, the leading eigenvalue $\lambda_{1}^{\boldsymbol{D}_{m}}=L_{m}(\mu_{m})$. Since $\mu_{m}\overset{\text{a.s.}}{\longrightarrow}\mu^{\ast}$, applying Lemma 3 stated in Appendix -E leads to

Recall from (32) that $\boldsymbol{P}_{m}$ is a principal submatrix of $\boldsymbol{D}_{m}$ obtained by deleting the first row and column of $\boldsymbol{D}_{m}$. It follows from the standard Cauchy interlacing theorem (see, e.g., [42, Theorem 4.3.8]) that

Applying [32, Lemma 3.1] (which is due to ), the upper edge of the support of the limiting spectral density of $\boldsymbol{P}_{m}$ is given by

where $\overline{\lambda}_{\alpha}$ is the minimizing point defined in (15). It follows that $\lambda_{2}^{\boldsymbol{P}_{m}}\overset{\text{a.s.}}{\longrightarrow}\zeta_{\alpha}(\overline{\lambda}_{\alpha})$ and $\lambda_{1}^{\boldsymbol{P}_{m}}\overset{\text{a.s.}}{\longrightarrow}\zeta_{\alpha}(\overline{\lambda}_{\alpha})$, and thus

by the interlacing inequalities in (60). Finally, by the constructions of $\psi_{\alpha}(\lambda)$ and $\zeta_{\alpha}(\lambda)$, we have $\zeta_{\alpha}(\lambda)>\zeta_{\alpha}(\overline{\lambda}_{\alpha})$ if and only if $\psi^{\prime}_{\alpha}(\lambda)>0$, and the proof is complete.

IV Sampling Ratios and Phase Transitions

In this section, we study the phase transition phenomena characterized in Theorem 1 in more detail. In particular, we prove Proposition 1 (as stated in Section II-B), which specifies the phase transitions and the asymptotic limits of the cosine similarities in terms of the sampling ratio $\alpha$.

By Theorem 1, whether the leading eigenvector $\boldsymbol{x}_{1}^{n}$ is asymptotically correlated or uncorrelated with the target vector $\boldsymbol{\xi}_{n}$ depends on the sign of the derivative $\psi^{\prime}_{\alpha}(\lambda)$ evaluated at a point $\lambda^{\ast}_{\alpha}$. And this point is uniquely defined through the equation $\zeta_{\alpha}(\lambda^{\ast}_{\alpha})=\phi(\lambda^{\ast}_{\alpha})$. Let $\overline{\lambda}_{\alpha}$, defined in (15), be the point at which the strictly convex function $\psi_{\alpha}(\lambda)$ achieves its minimum. Calculating the derivative of $\psi_{\alpha}(\lambda)$ and setting it to zero, we get

By the construction of the function $\zeta_{\alpha}(\lambda)$ in (16) and by the monotonicity of $\phi(\lambda)$, we can conclude that $\psi^{\prime}(\lambda^{\ast}_{\alpha})>0$ if and only if

where $\Delta(\lambda)$ is obtained by removing a common factor $\overline{\lambda}_{\alpha}$ from the difference $\psi_{\alpha}(\overline{\lambda}_{\alpha})-\phi(\overline{\lambda}_{\alpha})$ and by writing $\overline{\lambda}_{\alpha}$ simply as $\lambda$. Let $\Lambda$ be the set consisting of all the zero-crossings of $\Delta(\lambda)$ within the open interval $(\tau,\infty)$. Using (61), we can then establish a one-to-one mapping between points in $\Lambda$ and a set of critical values of the sampling ratios.

The set $\Lambda$ is nonempty. It contains a finite number of points, denoted by $\lambda_{c,1}\leq\lambda_{c,2}\leq\ldots\leq\lambda_{c,r}$ for some $r\geq 1$. Moreover,

We first show that $\Lambda$ is nonempty. For $\lambda>\tau$, applying the Cauchy-Schwartz inequality gives us

To study the function $\Delta(\lambda)$ as $\lambda\rightarrow\infty$, we note that

Since $\Delta(\lambda)$ is a continuous function, (64) and (65) imply that there must exist at least one zero-crossing.

Next, we show the upper bound given in (63). For any $\lambda_{c}\in\Lambda$, we have from (62) that

By assumption (A.4), $z$ is bounded within $[0,\tau]$. It follows that

IV-B Proof of Proposition 1

Write $\lambda_{c,\min}=\lambda_{c,1}$ and $\lambda_{c,\max}=\lambda_{c,r}$. The corresponding critical sampling ratios $\alpha_{c,\min}$ and $\alpha_{c,\max}$, as defined in (20), are obtained through the one-to-one mapping given in (61).

Fix $\alpha<\alpha_{c,\min}$. By the monotonicity of the mapping (61), the corresponding $\overline{\lambda}_{\alpha}$ is strictly less than the smallest zero-crossing point $\lambda_{c,1}$. From the proof of Lemma 2, we conclude that $\Delta(\overline{\lambda}_{\alpha})>0$, and thus $\zeta_{\alpha}(\cdot)$ and $\phi_{\alpha}(\cdot)$ intersects at a point $\lambda^{\ast}_{\alpha}<\overline{\lambda}_{\alpha}$. This implies that $\psi^{\prime}_{\alpha}(\lambda^{\ast}_{\alpha})<0$ and thus Theorem 1 gives us

Now fix $\alpha>\alpha_{c,\max}$, in which case $\overline{\lambda}_{\alpha}>\lambda_{c,\max}$. Since $\Delta(\overline{\lambda}_{\alpha})<0$, we must have $\lambda^{\ast}_{\alpha}>\overline{\lambda}_{\alpha}$ and thus $\psi^{\prime}_{\alpha}(\lambda^{\ast}_{\alpha})>0$. To derive the parametric form of $\rho(\alpha)$ given in the statement of the proposition, we note that

Thus, the equation $\zeta_{\alpha}(\lambda^{\ast}_{\alpha})=\phi(\lambda^{\ast}_{\alpha})$ becomes $\psi_{\alpha}(\lambda^{\ast}_{\alpha})=\phi(\lambda^{\ast}_{\alpha})$. Using the explicit definitions of these functions given in (13) and (14), we get

Substituting (70) and (71) into the asymptotic characterization (18) gives us (22), which, together with (68), provides a parametric representation of the function $\rho(\alpha)$.

Finally, we show that $\rho(\alpha)\rightarrow 1$ as $\alpha\rightarrow\infty$. From (68) and after some simple manipulations, we have

Since $\lambda^{\ast}_{\alpha}\rightarrow\infty$ as $\alpha\rightarrow\infty$, the above formula gives us

When the set $\Lambda$ consists of a single element, which is the case for many signal acquisition models we have studied, $\lambda_{c,\min}=\lambda_{c,\max}$ and thus $\alpha_{c,\min}=\alpha_{c,\max}$. There then exists a single critical sampling ratio $\alpha_{c}$ separating the uncorrelated phase from the correlated one. For $\alpha<\alpha_{c}$, the estimates from the spectral method is asymptotically orthogonal to $\boldsymbol{\xi}_{n}$; for $\alpha>\alpha_{c}$, the estimates will be concentrated on the surface of a right-circular cone whose generating lines make an angle $\theta=\arccos(\sqrt{\rho(\alpha)})$ to the target vector $\boldsymbol{\xi}_{n}$. The situation is more complicated when $\Lambda$ contains multiple zero-crossings, in which case a finite number of correlated and uncorrelated phases can alternatively take place between $\alpha_{c,\min}$ and $\alpha_{c,\max}$. A concrete example demonstrating this situation is shown in the next subsection.

IV-C Multiple Phase Transitions: an Example

where $0<\theta<1$, and $\mathcal{I}_{1},\mathcal{I}_{2}$ are two nonoverlapping intervals on the positive real axis. We also set $\kappa=\mathinner{\!\left\lVert\boldsymbol{\xi}_{n}\right\rVert}=1$, and thus $\boldsymbol{a}_{i}^{\top}\boldsymbol{\xi}_{n}$ has the same distribution as a standard normal random variable, denoted by $s$. Define

Choose $\theta=0.48$, $\mathcal{I}_{1}=[4.7947,4.9847]$, and $\mathcal{I}_{2}=[0.8995,0.8998]$. We then have $\beta_{1}=1.0086\times 10^{-6}$, $\beta_{2}=1.5970\times 10^{-4}$, $\omega_{1}=2.3976\times 10^{-5}$ and $\omega_{2}=1.2926\times 10^{-4}$. In this case, $\Delta(\lambda)$ turns out to have three zero-crossings:

By the mapping in (61), they correspond to three critical sampling ratios:

Using the characterization given in Theorem 1, we obtain the limiting values of the squared cosine similarity as a function of the sampling ratio $\alpha$. Figure 4 illustrates this function $\rho(\alpha)$. We can see that, when $\alpha<\alpha_{c,1}$, the estimates given by the spectral method are asymptotically uncorrelated with $\boldsymbol{\xi}_{n}$. When $\alpha$ is in the interval $(\alpha_{c,1},\alpha_{c,2})$, however, the function $\rho(\alpha)$ has a small “bump” (see the insert for a zoomed-in view), meaning that the estimates become asymptotically correlated with $\boldsymbol{\xi}_{n}$. However, the correlation returns to zero as $\alpha$ moves past the second phase transition point $\alpha_{c,2}$. Finally, when $\alpha>\alpha_{c,3}$, the estimates become correlated with $\boldsymbol{\xi}_{n}$ again, and $\rho(\alpha)$ tends to one as $\alpha\rightarrow\infty$.

It would be desirable to obtain a deeper understanding of the above phenomenon involving multiple phase transitions. The example provided here is purely theoretical, as its phase transitions take place at very large values of $\alpha$. It will be interesting to explore other possible examples of multiple phase transitions with more practical values of $\alpha$. Moreover, as most signal acquisition models we have studied seem to involve only a single phase transition point, it will be interesting to seek easy-to-verify conditions for the function $\Delta(\lambda)$ defined in (62) to have only one zero-crossing. We leave these as interesting open questions.

V Discussion

In this paper, we have presented a precise asymptotic characterization of the performance of a spectral method for estimating signals from generalized linear measurements with Gaussian sensing vectors. Our analysis also reveals a phase transition phenomenon that takes place at certain critical sampling ratios. Below a minimum threshold, estimates given by the methods are nearly orthogonal to the true signal $\boldsymbol{\xi}$, thus carrying no information; above a maximum threshold, the estimates become increasingly aligned with $\boldsymbol{\xi}$. The computational complexity of the spectral method is also markedly different in the two phases. Within the uncorrelated phase, the gap between the top two leading eigenvalues diminishes to zero. In contrast, a nonzero spectral gap emerges within the correlated phase. In this section, we close the paper by discussing some possible directions for extending and improving our results as well as their connections to related work in the literature.

The rate of convergence and more refined analysis. The performance of the spectral method was first studied in for the problem of phase retrieval. In that paper, it is shown that, for each $\delta\in(0,1)$, there is a constant $c_{1}(\delta)$ such that $\rho(\xi_{n},x_{1}^{n})>1-\delta$ with high probability when

Another possible direction to further refine our analysis is to consider second-order asymptotics at the level of central limit theorems (CLTs). See for instance for a related CLT analysis for the extreme eigenvalues of spiked covariance models.

Alternative initialization schemes. The spectral method considered in this paper is certainly not the only choice for initialization purposes. For example, an interesting alternative is the simple linear estimator studied in :

By using the moment calculations in [46, Proposition 1.1] and bounding high-order moments, one can easily obtain that

where $s$ and $z$ are the random variables defined in (5).

For cases where the function $g(s)$ is neither odd nor even, the choice between the spectral method and the linear estimator is not as clear-cut. The spectral method exhibits phase transition behaviors with its estimates in the uncorrelated phase at small values of $\alpha$. In contrast, as shown in (73), the performance of the linear estimator increases as a monotonic function of $\alpha$. As a result, in the regime of very small $\alpha$, the linear estimator will be preferable. For (moderately) larger values of $\alpha$, the comparison between the spectral method and the linear estimator cannot be easily made, as their performance also depends on the preprocessing function $\mathcal{T}(\cdot)$ used in (3) and (72).

The incorporation of priors. In this work, we assume that the target signal $\boldsymbol{\xi}_{n}$ is an arbitrary unknown (deterministic) signal. In many applications, the underlying signals satisfy additional constraints (such as sparsity). In , the authors considered a two-step scheme, where the initial linear estimate given in (72) is further projected onto a set which encapsulates one’s prior knowledge about $\boldsymbol{\xi}$. It will be interesting to consider and analyze similar projection schemes for the estimates obtained by the spectral method.

Universality and more realistic sensing vectors. Our asymptotic analysis assumes that the sensing vectors are real-valued i.i.d. Gaussian random vectors. Numerical simulations seem to suggest that the theoretical predictions given in Theorem 1 remain valid for more general random measurement ensembles and for complex-valued sensing vectors. To demonstrate this, we show in Figure 5 the results of applying the spectral method to estimate a $64\times 64$ cameraman image from phaseless measurements under Poisson noise:

where the bound $\tau$ is set to 5 and $\mathinner{\!\left\lVert\boldsymbol{\xi}\right\rVert}$ is normalized to 1 in our simulations. Two measurement ensembles are considered: real-valued sensing vectors whose elements are independent Rademacher ($\pm 1$) random variables, and complex-valued sensing vectors with elements drawn from the complex Gaussian distribution $\mathcal{N}(0,\tfrac{1}{2})+j\mathcal{N}(0,\tfrac{1}{2})$. We see from the figure that the theoretical predictions (the solid lines) have excellent agreement with simulation results for this moderately-sized problem, even though the sensing vectors can be non-Gaussian. Rigorously establishing the validity of our asymptotic predictions without the Gaussian assumption will be an important future work. Thanks to the deterministic characterization given in Proposition 2, this task boils down to showing that the result of Proposition 4 still holds when the sensing matrix consists of i.i.d. entries drawn from more general distributions. A related but more ambitious line of work will be to characterize the performance of the spectral method for structured and more practical sensing ensembles such as the coded diffraction scheme for phase retrieval with random modulation patterns.

and try to solve it via projected gradient descent

where $\mathcal{P}_{r}$ denotes projection onto the set of rank-$r$ matrices, and $\mu>0$ is the step size. As pointed out in , the spectral method studied in this paper can be viewed as the very first iteration of (75), if we start the algorithm from $\boldsymbol{X}_{0}=\boldsymbol{0}_{n\times n}$ and consider the special case of recovering a symmetric rank-one matrix (i.e., $r=1$, $p=n$) with $\boldsymbol{A}_{i}=\boldsymbol{a}_{i}\boldsymbol{a}_{i}^{\top}$. Thus, an interesting line of future research is to extend the results of this work, notably the key characterization given in Proposition 2, to more general settings with $r>1$ and $p\neq n$ and to other sensing matrices. Such extensions will be useful in applications such as low-rank matrix recovery, covariance estimation, and blind deconvolution.

In this appendix, we provide two sufficient conditions for Assumption (A.5) to hold.

Case 1: Suppose that the probability law of the random variable $z$ contains a point mass $c\,\delta(z-\tau)$ at its upper boundary $\tau$, where $c$ is some positive constant. This applies to the logistic regression model in Example 1, the subset algorithm (28) in Example 2, the noisy phase retrieval model in (74), and the quantization model described in Section IV-C.

which tends to $\infty$ as $\lambda$ approaches $\tau$ from the right.

Case 2: Suppose that there exist some positive constants $c$ and $\varepsilon$ such that the probability density function $p_{Z}(z)$ of $z$ and the conditional moment $h(z)$ are both bounded below by $c$ for all $z\in[\tau-\varepsilon,\tau]$. The model in (8) represents one such case. Under this setting,

-B Norm Estimation

The spectral initialization method estimates the orientation of the vector $\boldsymbol{\xi}_{n}$ but it provides no information about its norm, as the eigenvector $\boldsymbol{x}_{1}^{n}$ is always normalized. In many cases where the sensing vectors come from certain random ensembles, the norm $\mathinner{\!\left\lVert\boldsymbol{\xi}_{n}\right\rVert}$ can be accurately estimated from the measurements.

As a simple illustrative example, we can consider the (noiseless) phase retrieval problem: $y_{i}=(\boldsymbol{a}_{i}^{\top}\boldsymbol{\xi}_{n})^{2}$, where $\boldsymbol{\xi}_{n}$ is a deterministic unknown vector with $\kappa=\mathinner{\!\left\lVert\boldsymbol{\xi}_{n}\right\rVert}$, and the sensing vectors $\left\{\boldsymbol{a}_{i}\right\}$ are i.i.d. standard normal random vectors. Since $\boldsymbol{a}_{i}^{\top}\boldsymbol{\xi}_{n}\sim\mathcal{N}(0,\kappa^{2})$, the measurement $y_{i}$ can be represented as

where $s_{i}$ (for $1\leq i\leq m$) are i.i.d. standard normal random variables. A simple estimator of the norm is then

which is asymptotically consistent as $m\rightarrow\infty$.

More generally, consider an observation model $y_{i}\sim f(y\,|\,\boldsymbol{a}_{i}^{\top}\boldsymbol{\xi}_{n})$, where $f(\cdot\,|\,\cdot)$ is a conditional probability density function and $\boldsymbol{a}_{i}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\boldsymbol{I}_{n})$. Again, writing $\boldsymbol{a}_{i}^{\top}\boldsymbol{\xi}_{n}=\kappa s_{i}$ for i.i.d. normal random variables $\left\{s_{i}\right\}$, we can represent the probability distributions of the measurements $\left\{y_{i}\right\}$ as

We note that the estimator in (76) is a special case of (77). More generally, one could also estimate $\kappa$ by using maximum likelihood

whose asymptotic consistency can be established under standard conditions on the parametric density function $p_{\kappa}(y)$.

-C Proof of Proposition 2

By a suitable choice of a transformation matrix

where $\mathcal{I}_{1},\mathcal{I}_{2}$ are the two sets of indices defined in (42) and $\widetilde{\boldsymbol{q}}$ is a vector consisting of all the nonzero elements of $\boldsymbol{W}^{\top}\boldsymbol{q}$. Let $\lambda_{1}^{\widetilde{\boldsymbol{D}}}$ and $\widetilde{\boldsymbol{x}}_{1}$ be the largest eigenvalue of $\widetilde{\boldsymbol{D}}$ and an associated unit-norm eigenvector, respectively. Clearly, $\lambda_{1}^{\boldsymbol{D}}=\lambda_{1}^{\widetilde{\boldsymbol{D}}}$ and $(\boldsymbol{e}_{1}^{\top}\boldsymbol{x}_{1})^{2}=(\boldsymbol{e}_{1}^{\top}\widetilde{\boldsymbol{x}}_{1})^{2}$. Thus, we just need to consider $\widetilde{\boldsymbol{D}}$ in our proof.

Due to its block-diagonal form, the eigenvalues of $\widetilde{\boldsymbol{D}}$ is the union of those of its top-left submatrix

and those of its bottom-right submatrix $\operatorname{diag}\left\{\lambda_{i}^{\boldsymbol{P}}\mathrel{\mathop{\mathchar 58\relax}}i\in\mathcal{I}_{2}\right\}$. In particular,

The eigenvectors associated with $\operatorname{diag}\left\{\lambda_{i}^{\boldsymbol{P}}\mathrel{\mathop{\mathchar 58\relax}}i\in\mathcal{I}_{2}\right\}$ are easy to characterize. Clearly, each $\lambda_{i}^{\boldsymbol{P}},i\in\mathcal{I}_{2}$ is an eigenvalue of $\widetilde{\boldsymbol{D}}$, and it corresponds to an eigenvector $\boldsymbol{e}_{j(i)}$, where $j(i)\geq 3$ is the row index of $\lambda_{i}^{\boldsymbol{P}}$ in $\widetilde{\boldsymbol{D}}$.

The eigenvalues and eigenvectors of $\boldsymbol{S}$ can also be precisely characterized. Due to its shape, $\boldsymbol{S}$ is sometimes referred to in the literature as an arrowhead matrix . It can be shown (see for instance [pp. 94 – 97]) that $\lambda_{1}^{\boldsymbol{S}}$ is the unique point within the interval $\lambda>\max\left\{\lambda_{i}\mathrel{\mathop{\mathchar 58\relax}}i\in\mathcal{I}_{1}\right\}$ to satisfy the equation

where $R(\lambda)$ is the function defined in (38). (Alternatively, we can use the Laplace expansion to explicitly derive the characteristic polynomial of $\boldsymbol{S}$ as

Then, by following similar arguments as those used in the proof of Lemma 1, we can reach the characterization (80) about $\lambda_{1}^{\boldsymbol{S}}$.) Furthermore, let $\boldsymbol{x}_{1}^{\boldsymbol{S}}$ be a unit-norm eigenvector of $\widetilde{\boldsymbol{D}}$ associated with $\lambda_{1}^{\boldsymbol{S}}$. It is easily checked that

where $\boldsymbol{y}=(\lambda_{1}^{\boldsymbol{S}}\boldsymbol{I}-\operatorname{diag}\left\{\lambda_{i}^{\boldsymbol{P}}\right\}_{i\in\mathcal{I}_{1}})^{-1}\widetilde{\boldsymbol{q}}$ and $\boldsymbol{0}_{r}$ is a row vector of $r$ zeroes with $r$ being the cardinality of $\mathcal{I}_{2}$. It follows that

where $R^{\prime}(\lambda)$ denotes the derivative of the function $R(\lambda)$.

To show the claim of the proposition, we consider the following three cases.

Case 1: $\lambda_{1}^{\boldsymbol{S}}>\max\left\{\lambda_{i}^{\boldsymbol{P}}\mathrel{\mathop{\mathchar 58\relax}}i\in\mathcal{I}_{2}\right\}$. We choose $\mu^{\ast}=-1/R(\lambda_{1}^{\boldsymbol{S}})$, and thus $\lambda_{1}^{\boldsymbol{S}}=R^{-1}(-1/\mu^{\ast})$. It follows from Lemma 1 that

where the second equality is due to the fact that

Using the identity (80) for $\lambda_{1}^{\boldsymbol{S}}$, we can also verify that $\mu^{\ast}$ indeed satisfies the equation (44). (Its uniqueness is always guaranteed; see Remark 6 at the end of Section III-B.) The unit-norm leading eigenvector of $\widetilde{\boldsymbol{D}}$ in this case is the vector $\boldsymbol{x}_{1}^{\boldsymbol{S}}$ defined in (81). Since $L(\mu)=R^{-1}(-1/\mu)$ in a neighborhood of $\mu^{\ast}$, the function $L(\mu)$ is differentiable at $\mu^{\ast}$ and

Substituting (84) into (82) leads to (46).

Case 2: $\lambda_{1}^{\boldsymbol{S}}<\max\left\{\lambda_{i}^{\boldsymbol{P}}\mathrel{\mathop{\mathchar 58\relax}}i\in\mathcal{I}_{2}\right\}$, in which case $\lambda_{1}^{\widetilde{D}}=\max\left\{\lambda_{i}^{\boldsymbol{P}}\mathrel{\mathop{\mathchar 58\relax}}i\in\mathcal{I}_{2}\right\}=\lambda_{1}^{\boldsymbol{P}}$, where the last equality is due to (83). The corresponding leading eigenvector has nonzero elements only in its last $r$ entries, where $r$ is the cardinality of $\mathcal{I}_{2}$. Thus,

We set $\mu^{\ast}=(\lambda_{1}^{\boldsymbol{P}}-a)^{-1}$. (Note that we are guaranteed to have $\mu^{\ast}>0$. This can be verified by observing that $\lambda_{1}^{\boldsymbol{P}}>\lambda_{1}^{\boldsymbol{S}}=a-R(\lambda_{1}^{\boldsymbol{S}})>a$, where the equality is due to (80) and the last inequality follows from the fact that $R(\lambda)<0$.) Since $R(\lambda)$ is a strictly increasing function, we have

where the equality comes from (80). It then follows from Lemma 1 that $L(\mu^{\ast})=\lambda_{1}^{\boldsymbol{P}}=\lambda_{1}^{\widetilde{\boldsymbol{D}}}$ and, moreover, $\mu^{\ast}$ satisfies the equation (44).

To characterize the eigenvector, we note that $L(\mu)\equiv\lambda_{1}^{\boldsymbol{P}}$ in a neighborhood of $\mu^{\ast}$. We then have $L^{\prime}(\mu^{\ast})=0$, which, together with (85), leads to (46).