Non-negative Principal Component Analysis: Message Passing Algorithms and Sharp Asymptotics

Introduction

This approach is known to be consistent in low dimension. Let $\mathbf{v_{0}}$ be the solution of problem (1). If $n/p\to\infty$, then $\|\mathbf{v}_{1}-\mathbf{v_{0}}\|_{2}\to 0$ in probability [And63]. On the other hand, it is well understood that consistency can break dramatically in the high-dimensional regime $n=O(p)$. This phenomenon is crisply captured by the spiked covariance model [JL04, JL09], that postulates

where $\mathbf{v_{0}}$ has unit norm, $\mathbf{z}_{1},\mathbf{z}_{2},\dots\mathbf{z}_{p}$ are i.i.d. $p$-dimensional standard normal vectors $\mathbf{z}_{i}\sim{\sf N}(0,{\rm I}_{p}/n)$, and $\mathbf{u_{0}}=(u_{0,1},\dots,u_{0,n})^{{\sf T}}$ is a unit-norm vectorThe definition of [JL04] assumes $u_{i}\sim_{\text{i.i.d.}}{\sf N}(0,1/n)$ but for our purposes it is more convenient to consider $\mathbf{u}$ as a given deterministic vector. Equivalently, we can condition on $\mathbf{u}_{0}$.. The above model can also be written as

The spectral properties of the random matrix $\mathbf{X}$ defined by the Spiked Model have been studied in detail across statistics, signal processing and probability theory [BBAP05, BS06, BS06, Pau07, FP09, BGN12, CDMF12]. In the limit $n,p\to\infty$ with $p/n\to\alpha\in(0,\infty)$, the leading eigenvector $\mathbf{v}_{1}$ undergoes a phase transition:

In other words, Classical PCA contains information about the signal $\mathbf{v_{0}}$ if and only if the signal-to-noise ratio is above the threshold $\sqrt{\alpha}$. Below that threshold, the principal component is asymptotically orthogonal to the signal.

The failure of PCA has motivated significant effort aimed at developing better estimation methods. A recurring idea is to use additional structural information about the principal eigenvector $\mathbf{v_{0}}$, such as its sparsity [JL04, ZHT06] or its distribution (within a Bayesian framework) [Bis99, LU09]. Here we focus on the simplest type of structural information, namely we assume $\mathbf{v_{0}}$ is known to be non-negativeOf course the case in which $\mathbf{v_{0}}\in Q$ with $Q$ an arbitrary, known, orthant, can be reduced to the present one.. It is then natural to replace the Classical PCA problem with the following one (whereby we use the matrix $\mathbf{X}$ to represent the data):

Notice that this problem in non-convex and cannot be solved by standard singular value decomposition. Indeed it is in general NP-hard by reduction from maximum independent set [dKP02]. Two questions are therefore natural: given the additional complexity induced by the non-negativity constraint, does this constraint reduce the statistical error significantly? Are there efficient algorithms to solve the Non-negative PCA problem?

In this paper we answer positively to both questions within the spiked covariance model. Namely denoting by $\mathbf{v}^{+}$ the solution of the Non-negative PCA problem, we provide the following contributions:

We unveil a new phase transition phenomenon concerning $\mathbf{v}^{+}$ that is analogous to the classical one, see Eq. (3). Namely, for $\beta>\sqrt{\alpha/2}$, $\langle\mathbf{v}^{+},\mathbf{v_{0}}\rangle$ stays bounded away from , while, for $\beta<\sqrt{\alpha/2}$, there exists vectors $\mathbf{v_{0}}$ such that $\langle\mathbf{v}^{+},\mathbf{v_{0}}\rangle\to 0$ as $n,p\to\infty$.

Non-negative PCA is superior to classical PCA in this respect since $\sqrt{\alpha/2}<\sqrt{\alpha}$ strictly.

We prove an explicit formula for the asymptotic scalar product $\lim_{n\to\infty}\langle\mathbf{v}^{+},\mathbf{v_{0}}\rangle$. Non-negative PCA is superior to Classical PCA also in this respect. Namely $\langle\mathbf{v}^{+},\mathbf{v_{0}}\rangle$ is strictly larger than $|\langle\mathbf{v}_{1},\mathbf{v_{0}}\rangle|$ with high probability as $n,p\to\infty$.

Note that the non-negativity constraint breaks the rotational invariance of classical PCA (under the spiked model). As a consequence, not all spikes $\mathbf{v_{0}}$ are equally hard –or easy– to estimate. We use our theory to characterize the least favorable vectors $\mathbf{v_{0}}$.

We prove that (for any fixed $\delta>0$) a $(1-\delta)$ approximation to the non-convex optimization Non-negative PCA problem can be found efficiently with high probability with respect to the noise realization. Our algorithm has complexity of order $T_{\rm mult}\log(1/\delta)$, where $T_{\rm mult}$ is the maximum of the complexity of multiplying a vector by $\mathbf{X}$ or by $\mathbf{X}^{{\sf T}}$.

Technically, our approach has two components. We use Sudakov-Fernique inequality to upper bound the expected value of the Non-negative PCA optimization problem. We then define an iterative algorithm to solve the optimization problem, and evaluate the value achieved by the algorithm after any number $t$ of iterations. This provides a sequence of lower bounds which we prove converge to the upper bound as the number of iterations increase.

More precisely, we use an approximate message passing (AMP) algorithm of the type introduced in [DMM09, BM11]. Each iteration requires a multiplication by $\mathbf{X}$ and a multiplication by $\mathbf{X}^{{\sf T}}$ plus some lower complexity operations. While AMP is not guaranteed to solve the Non-negative PCA problem for arbitrary matrices $\mathbf{X}$, we establish the following properties:

Further the limit $\lim_{n,p\to\infty}\|\mathbf{X}\mathbf{v}^{t}\|^{2}_{2}=r(t)$ exists almost surely, and $r(t)$ can be computed explicitly as a function of the empirical law of entries of $\mathbf{v_{0}}$. Analogously, the asymptotic correlation $\lim_{n,p\to\infty}\langle\mathbf{v}^{t},\mathbf{v_{0}}\rangle=s(t)$ can be computed explicitly.

Denoting by $r_{*}$ the upper bound on the value of the optimization Non-negative PCA problem implied by Sudakov-Fernique inequality, we prove that $r(t)\geq(1-\delta)r_{*}$ for all $t\geq t_{0}(\delta)$ for some dimension-independent $t_{0}(\delta)$. This implies that Sudakov-Fernique inequality is asymptotically tight in the high-dimensional limit.

The asymptotic correlation converges to a limit as the number of iteration tends to infinity $s_{*}=\lim_{t\to\infty}s(t)$ (the convergence is, again, exponentially fast). Further, if we add the constraint $|\langle\mathbf{v},\mathbf{v_{0}}\rangle-s_{*}|\geq\delta$ to the Non-negative PCA optimization problem, Sudakov-Fernique’s upper bound on the resulting value is asymptotically smaller than $r_{*}$ for any $\delta>0$.

This implies that $\lim_{n\to\infty}\langle\mathbf{v}^{+},\mathbf{v_{0}}\rangle=s_{*}$.

Finally, we generalize our analysis to the case of symmetric matrices, namely assuming that data consist of a $n\times n$ symmetric matrix $\mathbf{X}$:

with $\mathbf{v_{0}}\geq 0$, $\|\mathbf{v_{0}}\|_{2}=1$. Here $\mathbf{Z}=\mathbf{Z}^{{\sf T}}$ is a noise matrix such that $(\mathbf{Z}_{ij})_{i\leq j}$ are independent with $\mathbf{Z}_{ij}\sim{\sf N}(0,1/n)$ for $i<j$ and $\mathbf{Z}_{ii}\sim{\sf N}(0,2/n)$.

In this case we study the analogue of the Non-negative PCA problem, namely

The non-negativity constraint on principal components arises naturally in many situations: we briefly discuss a few related areas. Let us emphasize that the theoretical understanding of the methods discussed below is much more limited than for Classical PCA.

where $k$ indexes such gene groups (or ‘layers’), and ${\boldsymbol{\rho}}^{(k)}$, ${\boldsymbol{\kappa}}^{(k)}$ indicate the level of participation of different samples or different genes in group $k$. These authors assume ${\boldsymbol{\rho}}^{(k)}_{i}{\boldsymbol{\kappa}}^{(k)}_{i}\in\{0,1\}$, but it is natural to relax this condition allowing for partial participation in group $k$, i.e. ${\boldsymbol{\rho}}^{(k)}_{i}{\boldsymbol{\kappa}}^{(k)}_{i}\in$, By a change of normalization, this constraint can be simplified to ${\boldsymbol{\rho}}^{(k)}_{i}{\boldsymbol{\kappa}}^{(k)}_{i}\geq 0$. Note a few differences with respect to our work:

We study a model with only one non-negative component. While Eq. (4) corresponds to a model with multiple $K\geq 1$ components, in practice several authors fit one ‘layer’ at a time, hence effectively reducing the problem to a single-component case.

Extending our analysis to the multiple component case will be the object of future work.

The non-negativity constraint is imposed in the model (4) on both components. This is a relatively straightforward modification of our setting.

Several studies (e.g. [LO02]) fit models of the form Eq. (4) using greedy optimization methods. Their conclusions are based on the unproven belief that these methods approximately solve the optimization problem. Our results (establishing convergence, with high probability, of an iterative method) provide some mathematical justification for this approach.

Neural signal processing. Neurons’ activity can be recorded through thin implanted electrodes. The resulting signal is a superposition of localized effects of single neurons (spikes). In order to reconstruct the single neuron activity, it is necessary to assign each spike to a specific neuron that created it, a process known as ‘spike sorting’ [Lew98, QNBS04, QP09]. Once spikes are aligned, the resulting data can be viewed as a matrix $\mathbf{X}=(\mathbf{X}_{ij})_{i\in[n],j\in[p]}$, where $i$ indexes the spikes and $j$ time (or a transform domain, e.g. wavelet domain).

In this context, principal component analysis is often used to project each row of $\mathbf{X}$ (i.e. each recorded spike) in a low dimensional space, or decomposing it as a sum of single neurons activity, see e.g. [BYS01, ZWZ+04, PMMP07]. Clustering may be carried out after dimensionality reduction. Note that each spike is a sum of single neuron activity with non-negative coefficients. In other words, the $i$-th row of $\mathbf{X}$ reads

where $\mathbf{v_{0}}^{(1)}$, …$\mathbf{v_{0}}^{(K)}$ are the signatures of $K$ neurons and $u_{0,ik}$ are non-negative coefficients.

Again, this corresponds to a multiple component version of the problem we study here. To the best of our knowledge, the non-negativity constraint has not been exploited in this context.

Non-negative matrix factorization. Initially introduced in the context of chemometrics [PT94, Paa97], non-negative matrix factorization attracted considerable interest because of its applications in computer vision and topic modeling. In particular, Lee and Seung [LS99] demonstrated empirically that non-negative matrix factorization successfully identifies parts of images, or topics in documents’ corpora.

A mathematical model to understand these findings was put forward in [DS03] and most recently studied, for instance, in [AGKM12]. Note that these results only apply under a no-noise or very-weak noise conditions, but for multiple components. Further, the aim is to approximate the original data matrix, rather than estimating the principal components.

In this sense, non-negative matrix factorization is the farther among all related areas to the scope of our work.

Approximate Message Passing. Approximate Message Passing algorithms proved successful as a fast first-order method for compressed sensing reconstruction [DMM09]. Their definition is inspired by ideas from statistical mechanics and coding theory [TAP77, MPV87, RU08], see also [Mon12] for further background. One attractive feature of AMP algorithms is that their high-dimensional asymptotics can be characterized exactly and in close form, through ‘state-evolution’ [BM11, JM13, BLM12]. Several applications and generalizations were developed by Rangan [Ran11], Schniter [VS11] and collaborators.

In particular Schniter and Cevher [SC11, PSC13] apply AMP the problem of reconstructing a vector from bilinear noisy observations, a problem that is mathematically equivalent to the one explored here. These authors consider however more complex Bayesian models, and evaluate performances through empirical simulations, while we characterize a fundamental threshold phenomenon in a worst case setting. Similar ideas were applied in [KMZ13] to the problem of dictionary learning, and in [VSM13] to hyperspectral imaging. Finally, Kabashima and collaborators [KKM+14] study low-rank matrix reconstruction using a similar approach, but focus on the case in which the rank scales linearly with the matrix dimensions.

2 Organization of the paper

In Section 2 we present formally our results, both for symmetric matrices and rectangular matrices. As mentioned above, the proof is obtained by establishing an upper bound on the value of the Non-negative PCA optimization problem using Sudakov-Fernique inequality, and a lower bound by analyzing an AMP algorithm. The upper bound is outlined in Section 3. Section 4 introduces formally AMP and its analysis, hence establishing the desired lower bound as well as the convergence properties of this algorithm. Section 5 presents a numerical illustration of the phase transition phenomenon, and of the behavior of our algorithm. Finally, Section 6 contains proofs, with some technical details deferred to the appendices.

Main results

In this section we present formally our results. For the sake of clarity , we consider first the case of symmetric (Wigner) matrices, and then the case of rectangular (or sample covariance, Wishart) matrices. Indeed formulæ for symmetric matrices are somewhat simpler. Before doing that, it is convenient to introduce some definitions. (For basic notations, we invite the reader to consult Section 2.4.)

Our results concern sequences of matrices $\mathbf{X}$ with diverging dimensions $n,p$, and are expressed in terms of the asymptotic empirical distribution of the entries of $\mathbf{v_{0}}$. This is formalized through the following definition.

converges weakly to $\mu$ and the second moment of $\mu_{\mathbf{x}(n)}$ converges as well or, equivalently, $\|\mathbf{x}(n)\|_{2}^{2}/n\to\int x^{2}\mu({\rm d}x)$.

With an abuse of terminology, we will say that $\{\mathbf{x}(n)\}_{n\geq 0}$ converges in empirical distribution to $X$ if $X$ is a random variable with law $\mu$.

Given a random variable $X$, we let $\mu_{X}$ denote its law. We next define a few functions of such a law.

Using ${\sf F}_{V}$ and ${\sf G}_{V}$ define the following ‘Rayleigh functions’

For $\beta\geq 0$, we also define ${\sf T}_{V}(\beta)$ as the unique non-negative solution of $x=\beta{\sf F}_{V}(x)$ and ${\sf S}_{V}(\beta,\alpha)$ as the unique non-negative solution of $x^{2}(1+\beta{\sf F}_{V}(x/\sqrt{\alpha})^{2})=\beta^{2}{\sf F}_{V}(x/\sqrt{\alpha})^{2}$.

Note that the above functions depend on the random variable ${V}$ only through its law $\mu_{{V}}$, but we prefer the notation –say– ${\sf F}_{{V}}$ to the more indirect ${\sf F}_{\mu_{V}}$. Existence and well-definedness of ${\sf T}_{V}$ and ${\sf S}_{V}$ are proved in Lemma 6.3 below. Further in Lemma 6.5 we prove that the functions ${\sf R}_{V}^{\text{sym}}$ respectively ${\sf R}_{V}^{\text{rec}}(\cdot,\alpha)$ have a unique maximum reached respectively at ${\sf T}_{V}(\beta)$ and at ${\sf S}_{V}(\beta,\alpha)$.

In the following, we will often state that an event holds almost surely as the dimensions of the random matrix $\mathbf{X}$ tend to infinity. It is understood that such statements hold with respect to the law of a sequence $\{\mathbf{X}_{n}\}_{n\geq 1}$ of independent random matrices distributed according to the Spiked Model or the Symmetric Spiked Model.

2 Symmetric matrices

This model has been studied in probability theory under the name of ‘low rank deformation of a Wigner matrix’. The following is a simplified version of the main theorem in [CDMF09].

Let $\mathbf{X}=\beta\mathbf{v_{0}}\mathbf{v_{0}}^{{\sf T}}+\mathbf{Z}$ be a rank-one deformation of the Gaussian symmetric matrix $\mathbf{Z}$, with $\mathbf{Z}_{ij}\sim{\sf N}(0,1/n)$ independent for $i<j$, and $\|\mathbf{v_{0}}\|_{2}=1$. Then we have, almost surely

Numerous refinements exist on this basic result, see for instance [CDMF09, Péc09, BGN11, BGGM11, CDMF+11, BGGM12, KY13, PRS13].

Our analysis provides a version of this theorem that holds for non-negative PCA, and is intriguingly similar to the original one. Its proof can be found in Appendix B.

Then (with the notation introduced in Definition 2.2), we have almost surely

The statement in Theorem 2 is dependent on the empirical distribution of the entries of $\mathbf{v_{0}}$. It is of special interest to characterize the least favorable situation, i.e. the distribution corresponding to the smallest scalar product $\langle\mathbf{v}^{+},\mathbf{v_{0}}\rangle$. This has two motivations: $(i)$ to guarantee the minimum value of $\langle\mathbf{v}^{+},\mathbf{v}_{0}\rangle$ achieved by a solution $\mathbf{v}^{+}$ of the optimization problem Symmetric non-negative PCA and $(ii)$ to describe the least favorable signal $\mathbf{v_{0}}$.

The worst-case scenario is realized for a particularly simple distribution, namely 2-atoms distribution, with an atom at . However, unlike in classical denoising [DJ94], the worst case mixture is not obtained by setting all the allowed coordinates to non-zero. In the following Theorem we are interested in the worst case among ${\bar{\varepsilon}}$-sparse signals, or equivalently in vector sequences $\{\mathbf{v_{0}}(n)\}_{n\geq 0}$ such that $\lim_{n\to\infty}\|\mathbf{v_{0}}(n)\|_{0}/n\leq{\bar{\varepsilon}}$, or $V\in{\cal P}_{{\bar{\varepsilon}}}$ since sparse signals are naturally interesting for applications.

Consider the Symmetric Spiked Model with the Symmetric non-negative PCA estimator.

If $\beta\leq 1/\sqrt{2}$, then there exists a sequence of vectors $\{\mathbf{v}_{0}(n)\}_{n\geq 0}$ such that $\lim_{n\to\infty}\|\mathbf{v_{0}}(n)\|_{0}/n=0$ and, almost surely,

For any $\beta>1/\sqrt{2}$, there exists $\varepsilon_{*}(\beta,{\bar{\varepsilon}})\in(0,{\bar{\varepsilon}}]$ such that the following is true. Let ${V}_{*}$ be the random variable with law

Then for any sequence of vectors $\{\mathbf{v_{0}}(n)\}_{n\geq 0}$ such that $\|\mathbf{v_{0}}(n)\|_{0}\leq n{\bar{\varepsilon}}$ we have, almost surely,

Equality holds if $\mathbf{v_{0}}(n)$ is the vector with $n\varepsilon_{*}$ non-zero entries, all equal to $1/\sqrt{n\varepsilon_{*}}$.

We defer this proof to Section 6.5. The worst case mixture $\varepsilon_{\#}(\beta)$ as well as the function ${\sf F}_{V_{*}}({\sf T}_{V_{*}}(\beta))$ can be expressed explicitly in terms of the Gaussian distribution function, see Section 6.5.

3 Rectangular matrices

We develop a very similar theory for the case of rectangular matrices. Our first result characterizes the value of the Non-negative PCA problem, and the estimation error, in analogy with Theorem 2. The proof can be found in Appendix B.

Let $\mathbf{X}=\sqrt{\beta}\mathbf{u_{0}}\mathbf{v_{0}}^{{\sf T}}+\mathbf{Z}$ be a rank-one deformation of the Gaussian matrix $\mathbf{Z}$ with $\mathbf{Z}_{ij}\sim{\sf N}(0,1/n)$ independent, and $\|\mathbf{u_{0}}\|_{2}=\|\mathbf{v_{0}}\|_{2}=1$. Further let $\sigma^{+}=\sigma^{+}(\mathbf{X})$ be the expected value of the Non-negative PCA problem, and $\mathbf{v}^{+}=\mathbf{v}^{+}(\mathbf{X})$ be any of the optimizers.

Then (with the notation introduced in Definition 2.2), we have almost surely

Finally, in the same fashion as Theorem 3, we can characterize the worst case signals $\mathbf{v_{0}}$.

Consider the Spiked Model, with the Non-negative PCA estimator.

If $\beta\leq\sqrt{\alpha/2}$, then there exists a sequence of vectors $\{\mathbf{v}_{0}(p)\}_{p\geq 1}$ such that $\lim_{p\to\infty}\|\mathbf{v_{0}}(p)\|_{0}/p=0$ and, almost surely,

For any $\beta>\sqrt{\alpha/2}$, there exists $\varepsilon_{rec,*}(\alpha,\beta,{\bar{\varepsilon}})\in(0,{\bar{\varepsilon}}]$ such that the following is true. Let ${V}_{*}$ be the random variable with law $(1-\varepsilon_{rec,*})\delta_{0}+\varepsilon_{rec,*}\,\delta_{1/\sqrt{\varepsilon_{rec,*}}}$. Then for any sequence of vectors $\{\mathbf{v}_{0}(p)\}_{p\geq 1}$, $\|\mathbf{v}_{0}(p)\|_{0}\leq p{\bar{\varepsilon}}$, we have. almost surely,

Equality holds if $\mathbf{v_{0}}(p)$ is the vector with $p\varepsilon_{*}$ non-zero entries, all equal to $1/\sqrt{p\varepsilon_{*}}$.

For the proof we refer to Section 6.5 which also contains explicit expressions to compute $\varepsilon_{{}^{\text{rec}},\#}$.

4 Additional notations

Upper bounds on non-negative PCA values

As mentioned above, Theorems 2 and 4 are proved in two steps. We establish an upper bound on the value of the optimization problem by using Sudakov-Fernique inequality and prove that the bound is tight by analyzing an iterative algorithm that solves the optimization problem.

The first statement concerns the Symmetric Spiked Model.

Consider the Symmetric Spiked Model, and let $\mathbf{v}^{+}=\mathbf{v}^{+}(\mathbf{X})$ be the Symmetric non-negative PCA estimator, with $\lambda^{+}=\lambda^{+}(\mathbf{X})$ the value of the corresponding optimization problem.

Then, under the assumptions of Theorem 2, we have

The second statement concern the (non-symmetric) Spiked Model.

Consider the Spiked Model and let $\mathbf{v}^{+}=\mathbf{v}^{+}(\mathbf{X})$ be the Non-negative PCA estimator, with $\sigma^{+}=\sigma^{+}(\mathbf{X})$ the value of the corresponding optimization problem.

Then, under the assumptions of Theorem 4, we have

The proof of Lemma 3.2 can be found in Section 6.2. The proof for the case of symmetric matrices, cf. Lemma 3.1, is completely analogous and we omit it.

While the above upper bounds are stated in asymptotic form, the proofs in Section 6 imply non-asymptotic upper bounds. Roughly speaking, the above upper bounds hold non-asymptotically up to an additive correction of order $1/\sqrt{n}$.

Approximate message passing algorithm

We use an algorithmic approach to prove a lower bound that matches the upper bound in Lemmas 3.1, 3.2. The algorithm is close in spirit to the usual power method that computes the leading eigenvector of a symmetric matrix $\mathbf{X}$ by iterating

Our approach differs substantially from this line of work. We develop an approximate message passing (AMP) algorithm that builds on ideas from statistical physics and graphical models [DMM09, Mon12]. Remarkably, exact high-dimensional asymptotics for these algorithms have been characterized in some generality using a method known as state evolution [BM11, BLM12]. We establish the desired lower bounds by applying this theory to our problem.

As before, we will start by considering the case of symmetric matrices and then move to rectangular matrices.

(The factor $\sqrt{n}$ is introduced here for future convenience.)

If we neglect the memory term $-{\sf b}_{t}\,f(\mathbf{v}^{t-1})$, the algorithm AMP-sym is extremely simple: It alternates between a power iteration, and an orthogonal projection onto the constraint set $\{\mathbf{v}~{}:~{}\mathbf{v}\geq 0~{},~{}\|\mathbf{v}\|\leq 1\}$. As proved in [BM11, BLM12] the memory term (‘Onsager term’) plays a crucial role in allowing for an exact high-dimensional characterization.

Note that $\mathbf{v}^{t}$ does not satisfy –in general– the positivity constraint. Indeed it is not the algorithm estimate of $\mathbf{v_{0}}$. After any number $t$ of iteration we construct the estimate

1.2 Asymptotic analysis

State evolution [DMM09, BM11, JM13, BLM12] is a mathematical technique that provides an exact distributional characterization of a class of algorithms that includes AMP-sym, under suitable probabilistic models for the matrix $\mathbf{X}$. In the present case, we will assume the Symmetric Spiked Model, with $\sqrt{n}\mathbf{v_{0}}$ converging in empirical distribution to a random variable ${V}$.

Informally, state evolution predicts that as $n\to\infty$, for any fixed $t\geq 1$, the state vector $\mathbf{v}^{t}$ is approximately normal with mean $\sqrt{n}\tau_{t}\,\,\mathbf{v_{0}}$ and covariance ${\rm I}_{n\times n}$. In other words, it can be viewed as a noisy version of the signal $\mathbf{v_{0}}$:

with ${G}\sim{\sf N}(0,1)$ independent of ${V}$. A formal statement is given below.

Consider the Symmetric Spiked Model, and assume that $\left\{\sqrt{n}\mathbf{v_{0}}(n)\right\}_{n\geq 0}$ converges in empirical distribution to a random variable $V$. Further, let $\{\tau_{t}\}_{t\geq 1}$ be defined by the state evolution recursion (35).

The proof of this result is a direct application of the results of [BM11, JM13] and can be found in Appendix A.1.

A second important result that follows from state evolution is that the sequence $\{\mathbf{v}^{t}\}_{t\geq 0}$ converges in the following asymptotic sense.

The proof of this statement is deferred to Appendix A.2.

As $t\to\infty$, $\tau_{t}\to\tau$, with $\tau$ the unique positive solution of the fixed point equation $\tau=\beta{\sf F}_{{V}}(\tau)$. By using the above two propositions, we then obtain the following lower bound, whose proof can be found in Section 6.3.

Consider the Symmetric Spiked Model, and assume that $\left\{\sqrt{n}\mathbf{v_{0}}(n)\right\}_{n\geq 0}$ converges in empirical distribution to a random variable $V$. Further, let $\{\widehat{\mathbf{v}}^{t}\}_{t\geq 0}$ be the AMP iterates as defined by AMP-sym and Eq. (33). Finally, let $\tau$ be the unique positive solution of the fixed point equation $\tau=\beta{\sf F}_{{V}}(\tau)$ (equivalently $\tau={\sf T}_{{V}}(\beta)$).

This provides the necessary lower bound that complements the upper bound based on Sudakov-Fernique inequality, cf. Section 3.

2 Rectangular matrices

After any number $t$ of iteration we construct the estimates

These satisfy the normalization and positivity constraints and are used as estimates of $\mathbf{u_{0}}$, $\mathbf{v_{0}}$.

2.2 Asymptotic analysis

We consider the high dimensional setup where $n\to\infty$, and $p=p(n)\to\infty$ with converging aspect ratio $p/n\to\alpha\in(0,1)$. We assume that $\ \{\sqrt{n}~{}\mathbf{u_{0}}(n)\}_{n\geq 0}$ converges in empirical distribution to $U$ and $\{\sqrt{p}~{}\mathbf{v_{0}}(p)\}_{p\geq 0}$ converges in empirical distribution to $V$.

The high dimensional asymptotics of $\mathbf{u}^{t}$, $\mathbf{v}^{t}$ is characterized –as in the symmetric case– through state evolution. We introduce the real-valued state evolution sequences $\{\vartheta_{t}\}_{t\geq 0}$ and $\{\mu_{t}\}_{t\geq 1}$ through the following recursion for $t\geq 0$

Consider the Spiked Model and assume that $\ \{\sqrt{n}~{}\mathbf{u_{0}}(n)\}_{n\geq 0}$ converges in empirical distribution to a random variable $U$ and $\{\sqrt{p}~{}\mathbf{v_{0}}(p)\}_{p\geq 0}$ converges in empirical distribution to a random variable $V$. Further, let $\{\mu_{t}\}_{t\geq 0}$, $\{\vartheta_{t}\}_{t\geq 1}$ be defined by the state evolution recursion SE-rec.

The proof is very similar to the one of Proposition 4.1 and is again a direct application of the results of [BM11, JM13]. We omit it to avoid redundancy.

We also have an analogous of Proposition 4.2.

We omit the proof, as it is very similar to the one of Proposition 4.4.

In the limit $t\to\infty$ (and assuming $\varepsilon>0$), the sequence defined in SE-rec converges to a nonzero fixed point $(\mu,\vartheta)$ satisfying the fixed point equations

We will prove that these equations admit a unique positive solution.

Considering $t\to\infty$ (after $n\to\infty$) we can thus prove the following.

Consider the Spiked Model and assume that $\ \{\sqrt{n}~{}\mathbf{u_{0}}(n)\}_{n\geq 0}$ converges in empirical distribution to a random variable $U$ and $\{\sqrt{p}~{}\mathbf{v_{0}}(p)\}_{p\geq 0}$ converges in empirical distribution to a random variable $V$. Further, let $\{\widehat{\mathbf{u}}^{t}\}_{t\geq 0}$, $\{\widehat{\mathbf{v}}^{t}\}_{t\geq 0}$ be the AMP estimates as defined by AMP-rec and Eq. (41) Finally, let $(\mu,\vartheta)$ be the only positive solution of the fixed point equations (44).

The proof of this theorem can be found in Section 6.3.

3 Computational complexity

As a direct consequence of the characterization of AMP established in Propositions 4.1 and 4.3, we can upper bound the number of iterations needed for Algorithms AMP-rec and AMP-sym to converge. We point out that the cost of each step of the AMP algorithms is dominated by a matrix vector multiplication. This operation can easily be parallelized and performed efficiently.

To be definite, we state the next result in the case of symmetric matrices. A completely analogous statement holds for rectangular matrices.

For any law $\mu_{V}\in{\cal P}$ and any $\delta>0$ there exists a constant $t_{0}(V,\delta)<\infty$ such that the following holds true. Under the assumptions of Proposition 4.1, let $\{\widehat{\mathbf{v}}^{t}\}_{t\geq 0}$ be the sequence of estimates produced by AMP. Then, for all fixed $t\geq t_{0}$ we have

The proof of this statement follows immediately from Theorem 2 and 6. A more careful treatment of error terms in the latter can be used to show that –indeed– $t_{0}(V,\delta)\leq C(V)\log(1/\delta)$ for some finite constant $C(V)$.

Notice that the computational cost of AMP is dominated by the one of matrix vector multiplications, call it $T_{\text{mult}}$. The above discussion indicates that the average-case complexity of the algorithms AMP-rec and AMP-sym is $O(T_{\text{mult}}\log 1/\delta)$.

Numerical illustration

We carried out numerical simulations on synthetic data generated following Symmetric Spiked Model. We use a signal $\mathbf{v_{0}}$ that takes two values:

where $S\subseteq[n]$ is of size $|S|=n\varepsilon$. It is immediate to see that the sequence $\{\sqrt{n}\mathbf{v_{0}}(n)\}_{n\geq 0}$ converges in empirical distribution to a random variable with distribution

In other words $\mu_{V}$ is the 2-points mixture $\mu_{V}=(1-\varepsilon)\delta_{0}+\varepsilon\,\delta_{1/\sqrt{\varepsilon}}$.

The predictions of Theorem 2 are stated in terms of the function ${\sf F}_{V}\equiv{\sf F}_{\varepsilon}$ that is rather explicit in this case. We have

We implemented the algorithm AMP-sym, and report in Figure 1 the results of numerical simulations with $n=10\,000$, sparsity level $\varepsilon\in\{0.001,0.1,0.8\}$, and signal-to-noise ratio $\beta\in\{0.05,0.10,\dots,1.5\}$. In each case we run AMP for $t=50$ iterations and plot the empirical average of $\langle\widehat{\mathbf{v}}^{t},\mathbf{v_{0}}\rangle$ over $32$ instances. The algorithm convergence is fast and –for our purposes– this value of $t$ is large enough so that $\tau_{t}\approx{\sf T}_{V}(\beta)$ and $\widehat{\mathbf{v}}^{t}\approx\mathbf{v}^{+}$. (See below for further evidence of this point.)

The results agree well with the asymptotic predictions of Theorem 2, namely with the curves reporting ${\sf F}_{V}({\sf T}_{V}(\beta))$. The figure also illustrates that sparse vectors (small $\varepsilon$) correspond to the least favorable signal in small signal-to-noise ratio. The value $\beta=1/\sqrt{2}$ corresponds to the phase transition.

2 Deviation from the asymptotic behavior

Theorem 2 and Proposition 4.1 predict the value of $\langle\mathbf{v_{0}},\mathbf{v}^{+}\rangle$ and $\langle\mathbf{v_{0}},\widehat{\mathbf{v}}^{t}\rangle$ in the limit $n\to\infty$. It is natural to question the validity of the prediction for moderate values of $n$.

In order to investigate this point, we performed numerical experiments with AMP by generating instances of the problem for several values of $n$ and compared the results with the asymptotic prediction of Eq. (51). The top left-hand frame in Figure 2 is obtained with $n=50,500,5000$, $\varepsilon=0.05$ and several value of $\beta$. For each point we plot the average of $\langle\widehat{\mathbf{v}}^{t},\mathbf{v_{0}}\rangle$ after $t=60$ iterations, over 32 instances.Already at $n=500$ the agreement is good, and improving with $n$.

In the top-right plot we plot the deviation between the empirical averages of $\langle\widehat{\mathbf{v}}^{t},\widehat{\mathbf{v}}\rangle\approx\langle\widehat{\mathbf{v}}^{+},\mathbf{v_{0}}\rangle$ (over 32 instances) and the asymptotic prediction ${\sf F}_{V}({\sf T}_{V}(\beta))$. The data suggest

In the bottom frames we plot the theoretical and empirical (for $n=1000$) values of $\langle\widehat{\mathbf{v}}^{t},\mathbf{v_{0}}\rangle$ for a grid of parameters $\beta,\varepsilon$. The difference between the two has average $1\cdot 10^{-3}$ and standard deviation $3\cdot 10^{-2}$.

3 Comparison with a convex relaxation

A natural convex relaxation for the Symmetric non-negative PCA problem is the semi-definite program

It is known [BAD09] that for $n\geq 5$ the completely positive cone is strictly included in the doubly non-negative cone

Hence in general this relaxation is not tight. The solution is a symmetric non-negative matrix $\hat{\mathbf{W}}$. We extract the leading eigenvector $\mathbf{v}_{1}(\hat{\mathbf{W}})$ and use its positive part as our approximation for $\mathbf{v}^{+}$.

In simulations we use CVX [GB10] to solve SDP, and compare the result to the output of AMP stopped after $t=50$ iterations. The interior point solver of CVX forces us to consider small problems. We use $n=50$, $\beta=1/\sqrt{2}$, $\varepsilon=0.3$, and average over $50$ instances.

Proofs

Positivity is immediate from the definition. The upper bound $1$ follows Cauchy-Schwarz inequality. To prove differentiability, we write ${\sf F}_{V}(x)={\sf Y}(x)/\sqrt{{\sf D}_{V}(x)}$ with

both differentiable (by dominated convergence) since $V$ and $G$ have bounded second moments, and strictly positive. Therefore ${\sf F}_{V}$ is differentiable.

A direct calculation yields the following relations

Using the last expression (and substituting the previous ones), we see that, to prove that ${\sf F}_{V}$ is increasing, it is sufficient to prove that

which directly follows from Cauchy-Schwarz inequality, even for $x<0$, and equality can not hold as $V$ and $G$ are independent.

Since ${\sf F}_{V}$ is strictly increasing, it follows that ${\sf G}_{V}$ is strictly decreasing.

Finally, the values at $x=0$ are obtained by simple calculus. The limits as $x\to\pm\infty$ follow by applying dominated convergence both to the numerator and to the denominator of ${\sf F}_{V}(x)$ (or ${\sf G}_{V}(x)$), after dividing both by $x$. ∎

In order to prove Eq. (63), first note that

admits a unique non-negative solution for each $\alpha,\beta>0$, which we denote by ${\sf T}_{V}(\beta)$ (for Eq. (70)) and ${\sf S}_{V}(\beta,\alpha)$ (for Eq. (70)).

Let us define the function ${\sf q}:x\mapsto{\sf q}(x)={\sf F}_{V}(x)/x$. We already know (by Lemma 6.1) that ${\sf F}_{V}(0)>0$, so $\lim_{x\to 0}{\sf q}(x)=\infty$. Also, since ${\sf F}_{V}(x)\leq 1$, we have $\lim_{x\to\infty}{\sf q}(x)=0$. Further ${\sf F}_{V}$ is differentiable and hence so is ${\sf q}$ on $(0,\infty)$. It is therefore sufficient to prove that ${\sf q}$ is strictly decreasing to prove existence and uniqueness of the solution of Eq. (70).

We will prove that $z\mapsto{\sf D}_{V}(\sqrt{z})$ is concave. This implies that ${\sf q}$ is decreasing: indeed, by the last equation we have

Applying the change of variable $x=\sqrt{z}$, we get

This shows that the derivative of ${\sf q}(x)$ is strictly negative provided that $\frac{{\rm d}^{2}\phantom{z}}{{\rm d}z^{2}}{\sf D}_{V}(\sqrt{z})$ is non-positive, or $z\mapsto{\sf D}_{V}(\sqrt{z})$ concave. Indeed the second term in the last expression is strictly negative because $\frac{{\rm d}\phantom{x}}{{\rm d}x}{\sf D}_{V}(x)>0$, cf. Lemma 6.1 and Eq. (57)

This concludes the proof that Eq. (70) admits a unique positive solution.

Consider now existence and uniqueness of solutions of Eq. (71). Note that this is equivalent to proving that for every $\beta>0$ there exists a unique $x>0$ such that

We know that ${\sf F}_{V}$ is an increasing function, so $x\mapsto 1/\sqrt{1+\beta{\sf F}_{V}(x/\sqrt{\alpha})^{2}}$ is a decreasing function taking positive values. The result follows by using monotonicity of ${\sf q}$.

In order to prove Eq. (72), notice that the lower bound follows from Lemma 6.1. For the upper bound, observe that ${\sf q}^{\prime}(x)=\left({\sf F}_{V}^{\prime}(x)-{\sf q}(x)\right)/x$. By evaluating it at ${\sf T}_{V}(\beta)$, and using ${\sf q}({\sf T}_{V}(\beta))=1/\beta$, ${\sf q}^{\prime}(x)\leq 0$, we get $\beta{\sf F}_{V}^{\prime}({\sf T}_{V}(\beta))<1$. ∎

Let ${\sf T}_{V}(\beta)$ and ${\sf S}_{V}(\beta,\alpha)$ be defined as per Eq. (6.3).

Then the function $x\mapsto{\sf R}_{V}^{\text{sym}}(x)$ is strictly increasing on $(0,{\sf T}_{V}(\beta))$ and strictly decreasing on $({\sf T}_{V}(\beta),+\infty)$. Similarly $x\mapsto{\sf R}_{V}^{\text{rec}}(x,\alpha)$ is strictly increasing on $(0,{\sf S}_{V}(\beta,\alpha))$ and strictly decreasing on $({\sf S}_{V}(\beta,\alpha),+\infty)$.

Recall that, by Lemma 6.1, ${\sf F}_{V}^{\prime}(x)>0$. Further, as per the proof of Lemma 6.3, $x\mapsto{\sf q}(x)={\sf F}_{V}(x)/x$ is strictly decreasing with ${\sf q}({\sf T}_{V}(\beta))=1/\beta$. This immediately implies the claim for ${\sf R}_{V}^{\text{sym}}$.

The argument for ${\sf R}_{V}^{\text{rec}}(\,\cdot\,,\alpha)$ is completely analogous. We write the derivative of ${\sf R}_{V}^{\text{rec}}(x,\alpha)$ with respect to $x$:

The claim follows again from ${\sf F}_{V}^{\prime}>0$, and using the properties of $x\mapsto{{\sf F}_{V}(x/\sqrt{\alpha})}/\{x\sqrt{1+\beta~{}{\sf F}_{V}(x/\sqrt{\alpha})^{2}}\}$ already discussed in the proof of Lemma 6.3. ∎

We proved in Lemma 6.3 that $x\mapsto\beta{\sf F}_{V}(x)$ is monotone increasing with $\beta{\sf F}_{V}(x)>x$ if $x<{\sf T}_{V}(\beta)$ and $\beta{\sf F}_{V}(x)<x$ if $x>{\sf T}_{V}(\beta)$. It follows that $\tau_{t+1}>\tau_{t}$ if $\tau_{t}<{\sf T}_{V}(\beta)$ and $\tau_{t+1}>\tau_{t}$ if $\tau_{t}>{\sf T}_{V}(\beta)$. Hence $\lim_{t\to\infty}\tau_{t}={\sf T}_{V}(\beta)$ Convergence is exponentially fast, i.e. $|{\sf T}_{V}(\beta)-\tau_{t}|\leq c_{0}\gamma_{0}^{t}$, since, by Lemma 6.3 $\beta{\sf F}^{\prime}_{V}({\sf T}_{V}(\beta))\in(0,1)$.

This proves the first second inequality. Note that ${\sf T}_{V}(\beta)$ is the global maximum of $x\mapsto{\sf R}_{V}^{\text{sym}}(x)$ and hence, in a neighborhood of ${\sf T}_{V}(\beta)$, $|{\sf R}_{V}^{\text{sym}}({\sf T}_{V}(\beta))-{\sf R}_{V}^{\text{sym}}(\tau_{t})|\leq c_{*}(\tau_{t}-{\sf T}_{V}(\beta))^{2}$ ∎

We state without proof the analogous result for the rectangular case. The argument is exactly the same as for the symmetric case.

Our results are stated in terms of ${\sf F}_{V}$ and ${\sf G}_{V}$, and depend on the law of $V$. However, when $\varepsilon(V)\to 0$, interestingly, two different phenomena occur. First, our results can be stated independently of law of ${V}$. Second, a phase transition occurs for a specific value of the signal-to-noise ratio $\beta$. This is stated formally below using the notion of uniform convergence introduced in Definition 2.3.

where f(z)\equiv[\Phi(z)-\Phi(0)\big{]}+z^{2}\big{[}\Phi(z)-1\big{]}+z\phi(z). Note that $f(0)=0$ and $f(z)$ is bounded, whence

which yields the desired uniform convergence of ${\sf D}_{V}$.

Next recall that ${\sf F}_{V}(x)=\frac{{\rm d}\phantom{x}}{{\rm d}x}\sqrt{{\sf D}_{V}(x)}$, cf. Eq. (57) and, by Lemma 6.1, $\sqrt{{\sf D}_{V}(x)}$ is strictly convex. We hence have, for all $\delta>0$,

The claim (79) follows by taking the limit $\varepsilon\to 0$ (using Eq. (78)) followed by $\delta\to 0$. The expression of $\lim_{\varepsilon(V)\to 0}{\sf G}_{V}(x)$ follows by taking the limit on the identity $x~{}{\sf F}_{V}(x)+{\sf G}_{V}(x)=\sqrt{{\sf D}_{V}(x)}$.

In order to prove Eq. (81), let ${\sf T}_{0}(\beta)$ denote the function on the right-hand side and assume by contradiction that there exists a sequence $\varepsilon_{n}\to 0$, probability measures $\mu_{V_{n}}\in{\cal P}_{\varepsilon_{n}}$ such that $\lim_{n\to\infty}{\sf T}_{V_{n}}(\beta)>x_{*}={\sf T}_{0}(\beta)+\delta$ for some $\delta>0$. As shown in the proof of Lemma 6.3, $x\mapsto x/{\sf F}_{V}(x)$ is monotone increasing. Using the definition we have, for all $n$ large enough

Taking the limit $n\to\infty$, and using Eq. (79), we get

that yields a contradiction by the definition of ${\sf T}_{0}$. Hence $\lim\sup_{\mu_{V}\in{\cal P}_{\varepsilon}}{\sf T}_{V}(\beta)\leq{\sf T}_{0}(\beta)$. The matching lower bound is proved in the same way.

Finally, the proof of Eq. (82) follows along the same lines. ∎

2 Upper bounds: Proof of Lemma 3.2

In this section we prove Lemma 3.2. As mentioned before, the proof of Lemma 3.1 is completely analogous and omitted.

The function $\mathbf{X}\mapsto M_{\mathbf{X}}(\mu)$ is Lipschitz continuous with Lipschitz constant $1$ (namely $|M_{\mathbf{X}}(\mu)-M_{\mathbf{X}^{\prime}}(\mu)|\leq\|\mathbf{X}-\mathbf{X}^{\prime}\|_{F}$). Hence, by Gaussian isoperimetry, we have

Further we claim that $\mu\mapsto M_{\mathbf{X}}(\mu)$ is uniformly continuous for $\mu\in$. Indeed if $\mathbf{v}(\mu)=\sqrt{1-\mu^{2}}\mathbf{v}_{\perp}(\mu)+\mu\mathbf{v}_{0}$ realizes the maximum over $\mathcal{W}_{\mu}$ (with $\langle\mathbf{v}_{\perp}(\mu),\mathbf{v}_{0}\rangle=0$), we have

Let ${{\cal I}}_{n}\equiv\{0,1/n,2/n,\dots\}\cap$ be a grid. By the above uniform continuity, we have, with probability at least $1-C_{2}e^{-n/C_{2}}$,

Using Eq. (93) and union bound over ${\cal I}_{n}$, we conclude that

where the last inequality holds for all $t\leq t_{0}$ with $t_{0}$ a suitable constant. In particular, by Borel-Cantelli we have, almost surely and in expectation,

In order to upper bound $\overline{M}(\mu)$, we apply Vitale’s extension of Sudakov-Fernique inequality (see e.g. [Vit00, Theorem 1] and [Cha05, Theorem 1] for a quantitative version) to the two processes $\{{\cal X}(\mathbf{u},\mathbf{v})\}$, $\{{\cal Y}(\mathbf{u},\mathbf{v})\}$ indexed by $(\mathbf{u},\mathbf{v})\in\mathcal{W}_{\mu}$ defined as follows:

The maximum in the last expression is achieved for

We next fix $\vartheta=\vartheta_{*}(\alpha,\beta)={\sf S}_{V}(\beta,\alpha)$, which is also the unique maximizer of $x\mapsto{\sf R}_{V}^{\text{rec}}(x,\alpha)$, as shown in Lemma 6.5. Note that Eq. (108) is strictly concave in $\mu\in$, with unique maximum at $\mu_{*}=(\vartheta_{*}/\beta)(1-\vartheta_{*}^{2}/\beta)^{-1/2}$. Substituting in Eq. (108), we get

where the last equality follows from the identity $\sqrt{{\sf D}_{V}(x)}=x{\sf F}_{V}(x)+{\sf G}_{V}(x)$, and from the equation $\vartheta_{*}=\beta{\sf F}_{V}(1+\beta{\sf F}_{V})^{-1/2}$ with ${\sf F}_{V}={\sf F}_{V}(\vartheta_{*}/\sqrt{\alpha})$ that holds by definition of $\vartheta_{*}={\sf S}_{V}(\beta,\alpha)$.

From Eq. (92), (99) and (110) we finally get

Next reconsidering Eq. (108) with $\vartheta=\vartheta_{*}$, we see that since the right-hand side is strictly concave in $\mu\in$, we can strengthen Eq. (110) to

for some $c_{*}>0$. We call $H(x)=c_{*}x^{2}$.

By Eq. (92) and (99) we have, almost surely,

This implies immediately Eq. (30) with $\Delta=H^{-1}$, since (as shown above) $\vartheta_{*}={\sf S}_{V}(\beta,\alpha)$, and $\mu_{*}={\sf F}_{V}({\sf S}_{V}(\beta,\alpha)/\sqrt{\alpha})$.

3 Lower bounds: Proofs of Theorem 6 and Theorem 7

In this section we prove lower bounds on the non-negative eigenvalue (singular value) that follows from the analysis of the AMP algorithm, namely Theorem 6 for symmetric matrices and Theorem 7 for rectangular matrices. The proofs are very similar in the two cases, hence we will provide details only in the case of rectangular matrices, and limit ourselves to pointing out differences arising in the symmetric setting.

4 Proof of Theorem 7

The last step follows from triangular inequality.

By Proposition 4.3 applied to $\psi(x,y)=x^{2}$, we have, almost surely,

By applying the same proposition to $\psi(x,y)=x_{+}^{2}$ we have

where the last equality follows from Stein’s lemma [Ste72]. Using together Eq. (127) and Eq. (129), we get

Using Eq. (124) and Eq. (130) in the upper bound (123), we get

Finally, substituting this result together with Eq. (126) and (130) in Eq. (120), we obtain

where the second equality follows by applying Proposition 4.3 to $\psi(x,y)=xy$ (for the numerator) and using Eq. (124) (for the denominator). Finally, the claim (46) follows by taking $t\to\infty$, and using Lemma 6.7.

The proof of claim (47) follows by the same argument and we omit it.

The proof in the symmetric case is very similar to the one for rectangular matrices, see Theorem 7. We limit ourselves to sketching the first steps. We have, using AMP-sym,

In addition, it follows from Proposition 4.3 that

5 Minimax analysis: proof of Theorems 3 and 5

In this section we prove that the least favorable vectors $\mathbf{v_{0}}$ are –asymptotically– of the following form: $(\mathbf{v_{0}})_{i}=1/\sqrt{|S|}$ for all $i\in S$, and $(\mathbf{v_{0}})_{i}=0$ otherwise, for some support $S\subseteq[p]$. Further, we characterize the least favorable size of the support $|S|$.

The proofs proceed by analyzing the expression in Theorem 2 and applying strong duality to a certain linear program over probability distributions, that is related to the function $\mu_{V}\mapsto{\sf F}_{V}(x)$. We start with some preliminary facts and definitions in Section 6.5.1. The key step is to reduce ourselves to two points mixtures: this is achieved in Section 6.5.2. Finally, in Sections 6.5.3 and 6.5.4, we use these results to prove Theorems 3 and 5. Since the proof of Theorem 5 is completely analogous to the one of Theorem 3, we will limit ourselves to mentioning the main differences.

In particular, when $v=1/\sqrt{\varepsilon}$ (and hence the above distribution has second moment equal to $1$), we write $\mu_{\varepsilon}=\mu_{\varepsilon,1/\sqrt{\varepsilon}}$. We also write –with a slight abuse of notation– ${\sf F}_{\varepsilon}(x)$ instead of ${\sf F}_{V}(x)$ when $V\sim\mu_{\varepsilon}$. Explicitly

We will also adopt the shorthand ${\sf T}_{\varepsilon}(\beta)={\sf T}_{V}(\beta)$ when $V\sim\mu_{\varepsilon}$.

We wil next establish two calculus lemmas that are useful for the following.

Let $h_{b}(v)=(\phi(v)/v)+b\Phi(v)$ denote the left-hand side of Eq. (147). Then

If $b\leq 1$, then we conclude that $h^{\prime}_{b}(v)<0$ for all $v>0$ and hence the equation $h_{b}(v)=a$ has at most one positive solution. If –on the other hand– $b>1$, then $h^{\prime}_{b}(v)<0$ for $v<v_{*}\equiv(b-1)^{-1/2}$ and $h^{\prime}_{b}(v)>0$ for $v>v_{*}$. It follows that the equation $h_{b}(v)=a$ has at most one solution in $(0,v_{*}]$ and at most one in $(v_{*},\infty)$. ∎

By simple calculus, we get $b(0)=0$, and the derivatives

Let us further recall the inequalities (valid for $x>0$)

Therefore the left-hand side of Eq. (149) is always strictly positive. Consider the right-hand side. By consulting special values of the normal distribution, we see that $b^{\prime}(1)=\phi(1)-2(1-\Phi(1))<-0.07<0$. By a change of variables we know that $x(1-\Phi(x))/\phi(x)=\int_{0}^{\infty}\exp(-z-z^{2}/(2x^{2})){\rm d}z$ or, equivalently

Since the term in curly brackets is decreasing in $x$, and is negative at $x=1$, we have $b^{\prime}(x)<0$ for all $x\geq 1$. Therefore the right-hand side of Eq. (149) is non-positive for $x\geq 1$. This proves the claim for $x\geq 1$, and we will assume hereafter $x\in(0,1)$.

Next notice that $0\leq b^{\prime\prime\prime}(x)\leq 3\phi(0)$ for $x\in(0,1)$. Therefore, by Taylor expansion and intermediate value theorem, we get, for $x\in(0,1)$,

The right-hand side is negative for $x\in(x_{0},x_{1})$ where

In particular $x_{0}\leq 2/3$, and $x_{1}>1$. It follows that the right-hand side of Eq. (149) is non-positive for $x\geq x_{0}$.

We will therefore restrict ourselves to considering $x\in(0,x_{0})\subseteq(0,2/3)$. Note that our claim can be equivalently written as

We will next develop, for $x\in(0,x_{0})$, a lower bound on the left-hand side, to be denoted by $l(x)$, and an upper bound on the right-hand side, to be denoted by $u(x)$ and prove that $l(x)\geq u(x)$. For the left hand side note that $b^{\prime\prime\prime}(x)\geq 0$ for $x\in(0,x_{0})$ and hence, again by Taylor expansion

For the right hand side note that $\Phi(x)-(1/2)\leq\phi(0)x$. Further $x\mapsto(1-\Phi(x))/\phi(x)$ is monotone decreasing. We therefore define

and obtain the upper bound $u(x)=\phi(0)x\,(1-d_{0}x)$. Hence

It is a straightforward exercise to check that indeed $d_{0}\phi(0)>(1/2)$ thus completing the proof. ∎

5.2 Reduction to two points mixtures

The main theorem of this Section shows that ${\sf F}_{V}(x)$ is minimized by probability measures $\mu_{V}$ that are mixture of at most two point masses.

Fix $x\geq 0$. Then for any random variable $V$ with probability distribution $\mu_{V}\in{\cal P}_{{\bar{\varepsilon}}}$, we have

The proof of this theorem is presented at the end of the section. Before getting to it, we’ll introduce a related problem. Note that

where ${\cal F}(x,y)$ is the value of a constrained optimization problem:

Here it is understood that ${\cal F}(x,y)=\infty$ if this problem is unfeasible.

By a rescaling of the objective function, and letting $W=xV$, we can rewrite the problem (164) as

The corresponding value is $x{\cal F}(x,y)$. Note that each of the functions $(f(w)/g(w))$, $1/g(w)$, $w^{2}/g(w)$ is bounded and Lipschitz continuous, with a finite limit as $w\to\infty$. This implies that the value $x{\cal F}(x,y)$ is achieved by a measure $\mu_{*}$ on the completed real line $[0,\infty]$, with total mass $y^{2}$. Indeed the family of normalized distributions on $[0,\infty]$ is tight and both the objective and the constraints are continuous in the weak topology. Hereafter, we shall assume this holds. Functions on $[0,\infty)$ are extended by continuity to $+\infty$.

because otherwise we could increase the lower bound Eq. (171) by increasing $\alpha$. Further, since the right-hand side is an analytic function of $w$, the infimum in Eq. (172) is achieved on a finite set $S_{*}\in[0,\infty]$, and the minimizer $\mu_{*}$ of problem (170) has support ${\rm supp}(\mu_{*})\subseteq S_{*}$ because otherwise the infimum in the lower bound (171) would not be achieved.

By Lemma 6.9 this has at most two solutions $w_{1},w_{2}$. If on the other hand $\gamma_{*}=1/2$, then the above equation reduces to $\Phi(w)=2\beta_{*}$ which has at most one solution. In both cases, at most one solution –call it $w_{*}$– is a local minimum of $h$.

We conclude that $S_{*}\subseteq\{0,w_{*},\infty\}$, and therefore the value of the problem (170) is achieved by a measure of the form

The three constraints imply the following relations

The proof is completed by the change of variables $p_{1}=g(w_{*})\,\varepsilon$, $p_{0}=({\bar{\varepsilon}}-\varepsilon)/2$, $p_{2}=qx^{2}$, $w_{*}=v_{*}x$, $p_{2}=qx^{2}$. With these substitutions Eq. (178) yields (166), and Eq. (179) yields (165). ∎

We are now in position to prove Theorem 8, that is the main result in this section.

By Lemma 6.11 and Eq. (163), we have, for any $\mu_{V}\in{\cal F}_{{\bar{\varepsilon}}}$,

where the infimum is over $q\in$, $\varepsilon\in(0,{\bar{\varepsilon}}]$, $v_{*}=\sqrt{(1-q)/\varepsilon}$. Our claim is equivalent to saying that the infimum on the right hand side is achieved when $q=0$.

Since $x>0$ is given, we will can regard the right-hand side as a function of $w=v_{*}x$ and $\varepsilon$, and substitute $qx^{2}=x^{2}-\varepsilon w^{2}$. We then define the function

which needs to be optimized over $\varepsilon\in(0,{\bar{\varepsilon}}]$, and $w\in[0,x/\sqrt{\varepsilon}]$. Our claim is equivalent to saying that the minimum cannot be in the interior of this domain.

Since $G$ is analytic in the mentioned domain, a minimum in the interior must satisfy $\partial_{w}G(w,\varepsilon)=\partial_{\varepsilon}G(w,\varepsilon)=0$. Simple calculus shows that these two conditions are equivalent –respectively– to:

Taking the ratio of these equations, we obtain the necessary condition

We conclude with a Corollary of Theorem 8. (Figure 4 provides an illustration of the argument used in the proof.)

Fix $\beta\in(0,\infty)$. Then for any random variable $V$ with probability distribution $\mu_{V}\in{\cal P}_{{\bar{\varepsilon}}}$, we have

Further, for any $\beta>1/\sqrt{2}$, the infimum on the right-hand side is achieved at some $\varepsilon_{*}\in(0,{\bar{\varepsilon}}]$.

Assume the claim (188) does not hold. Then there exists $\mu_{V}\in{\cal P}_{{\bar{\varepsilon}}}$ such that ${\sf T}_{V}(\beta)<{\sf T}_{\varepsilon}(\beta)$ for all $\varepsilon\in(0,{\bar{\varepsilon}}]$. Now, on the one hand, by definition we have

We therefore reached a contradiction, which proves the claim (188).

In order to prove that the infimum is achieved at some $\varepsilon_{*}\in(0,{\bar{\varepsilon}}]$, note that $\varepsilon\mapsto{\sf T}_{\varepsilon}(\beta)$ is clearly continuous and, by Lemma 6.8,

It is therefore sufficient to show that $\varepsilon\to{\sf T}_{\varepsilon}(\beta)$ is decreasing for $\varepsilon$ small enough. By an argument similar to the above, this follows if we show that $\varepsilon\mapsto{\sf F}_{\varepsilon}(x)$ is decreasing for $x={\sf T}_{0}(\beta)=\sqrt{\beta^{2}-(1/2)}$ and $\varepsilon$ small enough. Indeed using the definition (145) and recalling that $\Phi(w)=1-O(\phi(w))$ as $w\to\infty$, we get, for every fixed $x>0$

which is of course decreasing in $\varepsilon$ for $\varepsilon\in[0,c(x)]$ with $c(x)>0$. ∎