Subexponential-Time Algorithms for Sparse PCA
Yunzi Ding Dmitriy Kunisky Alexander S. Wein Afonso S. Bandeira
math.ST cs.CC cs.DS stat.ML stat.TH


Role of coherence in the plasmonic control of molecular absorption
    Stefano Corni
    December 30, 2023
==================================================================

empty


We study the computational cost of recovering a unit-norm sparse principal component x ∈^n planted in a random matrix, in either the Wigner or Wishart spiked model (observing either W + λ xx^⊤ with W drawn from the Gaussian orthogonal ensemble, or N independent samples from (0, I_n + β xx^⊤), respectively).
Prior work has shown that when the signal-to-noise ratio (λ or β√(N/n), respectively) is a small constant and the fraction of nonzero entries in the planted vector is x_0 / n = ρ, it is possible to recover x in polynomial time if ρ≲ 1/√(n). While it is possible to recover x in exponential time under the weaker condition ρ≪ 1, it is believed that polynomial-time recovery is impossible unless ρ≲ 1/√(n). We investigate the precise amount of time required for recovery in the “possible but hard” regime 1/√(n)≪ρ≪ 1 by exploring the power of subexponential-time algorithms, i.e., algorithms running in time exp(n^δ) for some constant δ∈ (0,1). For any 1/√(n)≪ρ≪ 1, we give a recovery algorithm with runtime roughly exp(ρ^2 n), demonstrating a smooth tradeoff between sparsity and runtime.
Our family of algorithms interpolates smoothly between two existing algorithms: the polynomial-time diagonal thresholding algorithm and the exp(ρ n)-time
exhaustive search algorithm.
Furthermore, by analyzing the low-degree likelihood ratio, we give rigorous evidence suggesting that the tradeoff achieved by our algorithms is optimal.








§ INTRODUCTION




 §.§ Spiked Matrix Models


Since the foundational work of Johnstone <cit.>, spiked random matrix ensembles have been widely studied throughout random matrix theory, statistics, and theoretical data science. These models describe a deformation of one of several canonical random matrix distributions by a rank-one perturbation or “spike,” intended to capture a signal corrupted by noise. Spectral properties of these spiked models have received much attention in random matrix theory <cit.>, leading to a theoretical understanding of methods based on principal component analysis (PCA) for recovering the direction of the rank-one spike <cit.>. Spiked matrix models have also found more specific applications to problems such as community detection in graphs (see, e.g., <cit.>, or  <cit.> for surveys) and synchronization over groups (see, e.g., <cit.>).

We will study two classical variants of the spiked matrix model: the Wigner and Wishart models.
The models differ in how noise is applied to the signal vector.
In either case, let x ∈^n be the signal vector (or “spike”).
We will either have x deterministic with x_2 = 1, or x ∈^n random for each n with x_2 → 1 in probability as n →∞.



  * Spiked Wigner Model. 
Let λ > 0. Observe Y = W + λ xx^⊤, where W ∈^n × n is drawn from the Gaussian orthogonal ensemble (n), i.e., W is symmetric with entries distributed independently as W_ii∼(0,2/n) for all 1 ≤ i ≤ n, and W_ij = W_ji∼(0,1/n) for all 1 ≤ i < j ≤ n.


  * Spiked Wishart Model. 
Let β > 0 and N ∈.
Observe N samples y^(1),y^(2),…,y^(N)∈^n drawn independently from (0,I_n + β xx^⊤).
The ratio of dimension to number of samples is denoted γ n/N. We will focus on the high-dimensional regime where γ converges to a constant as n →∞. We let Y denote the sample covariance matrix Y = 1/N∑_i=1^N y^(i)y^(i)^⊤.


Each of these planted models has a corresponding null model, given by sampling from the planted model with either λ = 0 (Wigner) or β = 0 (Wishart). 

We are interested in the computational feasibility of the following two statistical tasks, to be performed given a realization of the data (either Y or {y^(1),…,y^(N)}) drawn from either the null or planted distribution.


  * Detection. Perform a simple hypothesis test between the planted model and null model. We say that strong detection is achieved by a statistical test if both the type-I and type-II errors tend to 0 as n →∞.

  * Recovery. Estimate the spike x given data drawn from the planted model. We say that a unit-norm estimator x∈^n achieves weak recovery if ⟨x,x ⟩^2 remains bounded away from zero with probability tending to 1 as n →∞.[We will also consider stronger notions of recovery: strong recovery is ⟨x̂, x ⟩^2 → 1 as n →∞ and exact recovery is x̂ = x with probability 1-o(1).] (Note that we cannot hope to distinguish between the planted models with signals x and -x.)

For high-dimensional inference problems such as the spiked Wigner and Wishart models, these two tasks typically share the same computational profile: with a given computational time budget, strong detection and weak recovery are possible in the same regions of parameter space.



 §.§ Principal Component Analysis


Simple algorithms for both detection and recovery are given by principal component analysis (PCA) of the matrix Y.
For detection, one computes and thresholds the maximum eigenvalue λ_max(Y) of Y, while for recovery one estimates x using the leading eigenvector v_max(Y).
Both the spiked Wishart and Wigner models are known to exhibit a sharp transition in their top eigenvalue as the model parameters vary. For the Wishart model, the celebrated “BBP transition” of Baik, Ben Arous, and Péché <cit.> states that the maximum eigenvalue of the sample covariance matrix Y emerges from the Marchenko–Pastur-distributed bulk if and only if β^2 > γ. Similarly, in the Wigner model, the maximum eigenvalue of Y emerges from the semicircular bulk if and only if λ > 1 <cit.>. More formally, the following statements hold.



Consider the spiked Wigner model Y = W+λ xx^⊤ with x = 1 and λ > 0 fixed. Then as n→∞,


  * if λ≤ 1, λ_max(Y) → 2 almost surely, and v_max(Y),x^2 → 0 almost surely (where v_max denotes the leading eigenvector);

  * if λ > 1, λ_max(Y) →λ + λ^-1 > 2 almost surely, and v_max(Y),x ^2 → 1-λ^-2 > 0 almost surely.





Let Y denote the sample covariance matrix in the spiked Wishart model with x = 1, β > 0 fixed, and N = N(n) such that γ n/N converges to a constant γ̅> 0 as n →∞. Then as n →∞,


  * if β^2 ≤γ̅, λ_max(Y) → (1+√(γ̅))^2 almost surely, and v_max(Y),x^2 → 0 almost surely (where v_max denotes the leading eigenvector);

  * if β^2 > γ̅, λ_max(Y) → (1+β)(1+γ̅/β) > (1+√(γ̅))^2 almost surely, and v_max(Y),x ^2 → (1-γ̅/β^2)(1+γ̅/β) > 0 almost surely.




We define the signal-to-noise ratio (SNR)  λ̂ by

    λ̂{[                     λ  in the Wigner model,;              β / √(γ) in the Wishart model. ].

Theorems <ref> and <ref> then characterize the performance of PCA in terms of λ̂. Namely, thresholding the largest eigenvalue of Y succeeds at strong detection when λ̂> 1 and fails when λ̂≤ 1; similarly, the top eigenvector succeeds at weak recovery when λ̂> 1 and fails when λ̂≤ 1. For some distributions of x, including the spherical prior (x drawn uniformly from the unit sphere) and the Rademacher prior (each entry x_i drawn i.i.d. from Unif(± 1/√(n))), it is known that the PCA threshold is optimal, in the sense that strong detection and weak recovery are statistically impossible (for any test or estimator, regardless of computational cost) when λ̂< 1 <cit.>.



 §.§ Sparse PCA


Sparse PCA, a direction initiated by Johnstone and Lu <cit.>, seeks to improve performance when the planted vector is known to be sparse in a given basis. While various sparsity assumptions have been considered, a simple and illustrative one is to take x drawn from the sparse Rademacher prior, denoted 𝒳_n^ρ, in which each entry x_i is distributed independently as

    x_i = {1/√(ρ n)     with probability ρ/2
    
    -1/√(ρ n)     with probability ρ/2
    
    0      with probability 1-ρ
    .

for some known sparsity parameter ρ∈ (0,1], which may depend on n.[We analyze our algorithms for a more general set of assumptions on x; see Definition <ref>.] The normalization ensures x→ 1 in probability as n →∞.

Consider the Wishart model (we will see that the Wigner model shares essentially the same behavior) in the regime λ̂= Θ(1) with λ̂< 1 (so that ordinary PCA fails at weak recovery). The simple diagonal thresholding algorithm proposed by <cit.> estimates the support of x by identifying the largest diagonal entries of Y. Under the condition[We use A ≲ B to denote A ≤ CB for some constant C, and use A ≪ B to denote A ≤ B/polylog(n).] ρ≲ 1/√(n log n), this has been shown <cit.> to achieve exact support recovery, i.e., it exactly recovers the support of x with probability tending to 1 as n →∞ (and once the support is known, it is straightforward to recover x). The more sophisticated covariance thresholding algorithm proposed by <cit.> has been shown <cit.> to achieve exact support recovery when ρ≲ 1/√(n).

On the other hand, given unlimited computational power, an exhaustive search over all possible support sets of size ρ n achieves exact support recovery under the much weaker assumption ρ≲ 1/log(n) <cit.>. Similarly, strong detection and weak recovery are statistically possible even when ρ is a sufficiently small constant (depending on β,γ) <cit.>, and the precise critical constant ρ^*(β,γ) is given by the replica formula from statistical physics (see, e.g., <cit.>, or <cit.> for a survey). However, no polynomial-time algorithm is known to succeed (for any reasonable notion of success) when ρ≫ 1/√(n), despite extensive work on algorithms for sparse PCA <cit.>. 

In fact, a growing body of theoretical evidence suggests that no polynomial-time algorithm can succeed when ρ≫ 1/√(n) <cit.>. Such evidence takes the form of reductions <cit.> from the planted clique problem (which is widely conjectured to be hard in certain regimes <cit.>), as well as lower bounds against the sum-of-squares hierarchy of convex relaxations <cit.>. Thus, we expect sparse PCA to exhibit a large “possible but hard” regime when 1/√(n)≪ρ≪ 1. Statistical-to-computational gaps of this kind are believed to occur and have been studied extensively in many other statistical inference problems, such as community detection in the stochastic block model <cit.> and tensor PCA <cit.>.



 §.§ Our Contributions


In this paper, we investigate precisely how hard the “hard” region (1/√(n)≪ρ≪ 1) is in sparse PCA. We consider subexponential-time algorithms, i.e., algorithms with runtime exp(n^δ + o(1)) for fixed δ∈ (0,1). We show a smooth tradeoff between sparsity (governed by ρ) and runtime (governed by δ). More specifically, our results (for both the Wishart and Wigner models) are as follows.


  * Algorithms. For any δ∈ (0,1), we give an algorithm with runtime exp(n^δ + o(1)) that achieves exact support recovery, provided ρ≪ n^(δ-1)/2.

  * Lower bounds. Through an analysis of the low-degree likelihood ratio (see Section <ref>), we give formal evidence suggesting that the above condition is essentially tight in the sense that no algorithm of runtime exp(n^δ + o(1)) can succeed when ρ≫ n^(δ-1)/2. (Our results are sharper than the sum-of-squares lower bounds of <cit.> in that we pin down the precise constant δ.)

Our algorithm involves exhaustive search over subsets of [n] of cardinality ℓ≈ n^δ. The case ℓ = 1 is diagonal thresholding (which is polynomial-time and succeeds when ρ≲ 1/√(n log n)) and the case ℓ = ρ n is exhaustive search over all possible spikes (which requires time exp(ρ n^1+o(1)) and succeeds when ρ≲ 1/(log n)). As ℓ varies in the range 1 ≤ℓ≤ρ n, our algorithm interpolates smoothly between these two extremes. For a given ρ in the range 1/√(n)≪ρ≪ 1, the smallest admissible choice of ℓ is roughly ρ^2 n, yielding an algorithm of runtime exp(ρ^2 n^1+o(1)). 

Our results extend to the case λ̂≪ 1, e.g., λ̂= n^-α for some constant α > 0. In this case, provided ρ≪λ̂^2 (which is information-theoretically necessary <cit.>), there is an exp(n^δ + o(1))-time algorithm if ρ≪λ̂n^(δ-1)/2, and the low-degree likelihood ratio again suggests that this is optimal. In other words, for a given ρ in the range λ̂/√(n)≪ρ≪λ̂^2, we can solve sparse PCA in time exp(λ̂^-2ρ^2 n^1+o(1)).

The analysis of our algorithm applies not just to the sparse Rademacher spike prior, but also to a weaker set of assumptions on the spike that do not require all of the nonzero entries to have the same magnitude. Our algorithm is guaranteed (with high probability) to exactly recover both the support of x and the signs of the nonzero entries of x. Once the support is known, it is straightforward to estimate x via the leading eigenvector of the appropriate submatrix.

In independent work <cit.>, a different algorithm for sparse PCA was proposed and shown to have essentially the same subexponential runtime as ours. Also, prior work <cit.> gave a subexponential-time algorithm for certifying the restricted isometry property, that is somewhat similar in spirit to our algorithm for sparse PCA.


Certain problems besides sparse PCA have a similar smooth tradeoff between subexponential runtime requirements and statistical power. These include refuting random constraint satisfaction problems <cit.> and tensor PCA <cit.>. In contrast, other problems have a sharp threshold at which they transition from being solvable in polynomial-time to (conjecturally) requiring essentially exponential time: exp(n^1-o(1)). Examples of this behavior can occur at the spectral transition at λ̂= 1 in the spiked Wishart and Wigner matrix models (see <cit.>) as well as at the Kesten–Stigum threshold in the stochastic block model (see <cit.>).





 §.§ Background on the Low-Degree Likelihood Ratio




A sequence of recent work on the sum-of-squares hierarchy <cit.> has led to the development of a remarkably simple method for predicting the amount of computation time required to solve statistical tasks. This method—which we will refer to as the low-degree method—is based on analyzing the so-called low-degree likelihood ratio, and is believed to be intimately connected to the power of sum-of-squares (although formal implications have not been established). We now give an overview of this method; see <cit.> for more details.

We will consider the problem of distinguishing two simple hypotheses _n and _n, which are probability distributions on some domain Ω_n = ^d(n) with d(n) = poly(n). The idea of the low-degree method is to explore whether there is a low-degree polynomial f_n: Ω_n → that can distinguish _n from _n.

We call _n the “null” distribution, which for us will always be i.i.d. Gaussian (see Definitions <ref> and <ref>).
_n induces an inner product on L^2 functions f: Ω_n → given by f,g _L^2(_n) = _Y ∼_n[f(Y)g(Y)], and a norm f_L^2(_n)^2 =  f,f _L^2(_n). For D ∈, let [Y]_≤ D denote the multivariate polynomials Ω_n → of degree at most D. For f: Ω_n →, let f^≤ D denote the orthogonal projection (with respect to ·,·_L^2(_n)) of f onto [Y]_≤ D. The following result then relates the distinguishing power of low-degree polynomials (in a certain L^2 sense) to the low-degree likelihood ratio.


Let  and  be probability distributions on Ω = ^d. Suppose  is absolutely continuous with respect to , so that the likelihood ratio L = d/d is defined. Then

    max_f ∈[Y]_≤ D∖{0}_Y ∼[f(Y)]/√(_Y ∼[f(Y)^2]) = L^≤ D_L^2().


(The proof is straightforward: the fraction on the left can be written as ⟨ f,L ⟩_L^2(_n) / f_L^2(_n), so the maximizer is f = L^≤ D.) 
The left-hand side of (<ref>) is a heuristic measure of how well degree-D polynomials can distinguish  from : if this quantity is O(1) as n →∞, this suggests that no degree-D polynomial can achieve strong detection (and indeed this is made formal by Theorem 4.3 of <cit.>).
The right-hand side of (<ref>) is the norm of the low-degree likelihood ratio (LDLR), which can be computed or bounded in many cases, making this heuristic a practical tool for predicting computational feasibility of hypothesis testing.

The key assumption underlying the low-degree method is that, for many natural distributions _n and _n, degree-D polynomials are as powerful as algorithms of runtime n^Θ̃(D), where Θ̃ hides factors of log n. This is captured by the following informal conjecture, which is based on <cit.>; in particular, see Hypothesis 2.1.5 of <cit.>.


    
    Suppose t: →.
    For “nice” sequences of distributions _n and _n, if L_n^≤ D(n)_L^2(_n) remains bounded as n →∞ whenever D(n) ≤ t(n) ·polylog(n), then there exists no sequence of functions f_n: Ω_n →{,} with f_n computable in time n^t(n) that strongly distinguishes _n and _n, i.e., that satisfies
    
    lim_n →∞_n[ f_n(Y) = ] = lim_n →∞_n[ f_n(Y) = ] = 1.


On a finer scale, it is conjectured <cit.> that if for some ε > 0 we have D(n) ≥log^1+ε(n) and L_n^≤ D(n)_L^2(_n) = O(1), then no polynomial-time algorithm can strongly distinguish _n from _n. In practice, it seems that the converse of Conjecture <ref> often holds as well, in the sense that if L_n^≤ D(n)_L^2(_n) = ω(1) for some D(n) = t(n) / polylog(n), then there is an n^t(n)-time distinguishing algorithm (however, see Remark <ref> for one caveat).

Calculations with the LDLR have been carried out for problems such as community detection <cit.>, planted clique <cit.>, the spiked Wishart model <cit.>, the spiked Wigner model <cit.>, and tensor PCA <cit.> (tensor PCA exhibits a subexponential-time tradeoff similar to sparse PCA; see <cit.>). In all of the above cases, the low-degree predictions coincide with widely-conjectured statistical-versus-computational tradeoffs.

Various leading algorithmic approaches can be approximated by low-degree polynomials and are thus ruled out by low-degree lower bounds of the form L_n^≤ D(n)_L^2(_n) = O(1). These approaches include a general class of spectral methods (see Theorem 4.4 of <cit.>) as well as the algorithms that we present in this paper (see Remark <ref>). The low-degree predictions are also conjectured to coincide with the power of the sum-of-squares hierarchy and are in particular connected to the pseudo-calibration approach <cit.>; see <cit.>. We refer the reader to Section 4 of <cit.> for further discussion of the implications (both formal and conjectural) of low-degree lower bounds.

Conjecture <ref> is informal in the sense that we have not specified the meaning of “nice" _n and _n. Roughly speaking, highly-symmetric high-dimensional problems are considered “nice” so long as _n and _n have at least a small amount of noise in order to rule out brittle high-degree algorithms such as Gaussian elimination. (In particular, we consider spiked Wigner and Wishart to be “nice.”) Conjecture 2.2.4 of <cit.> is one formal variant of the low-degree conjecture, although it uses the more refined notion of coordinate degree and so does not apply to the calculations in this paper.

We remark that if L_n_L^2(_n) = O(1) (the D = ∞) case, then it is statistically impossible to strongly distinguish _n and _n; this is a commonly-used second moment method (see e.g., <cit.>) of which Conjecture <ref> is a computationally-bounded analogue.

In this paper we give tight computational lower bounds for sparse PCA, conditional on Conjecture <ref>. Alternatively, one can view the results of this paper as a “stress test” for Conjecture <ref>: we show that Conjecture <ref> predicts a certain statistical-versus-computational tradeoff and this indeed matches the best algorithms that we know.





  
Organization.
The remainder of the paper is organized as follows.
In Section <ref>, we present our subexponential-time algorithms and our lower bounds based on the low-degree likelihood ratio.
In Section <ref>, we give proofs for the correctness of our algorithms. In Section <ref>, we give proofs for our analysis of the low-degree likelihood ratio.



  
Notation.
We use standard asymptotic notation O(·), Ω(·), Θ(·), always pertaining to the limit n →∞. We also use Õ(B) to mean O(B ·polylog(n)) and Ω̃(B) to mean Ω(B / polylog(n)). Also recall that f(n) = o(g(n)) means f(n)/g(n) → 0 as n →∞ and f(n) = ω(g(n)) means f(n)/g(n) →∞ as n →∞. An event occurs with high probability if it occurs with probability 1-o(1). We sometimes use the shorthand A ≲ B to mean A ≤ CB for an absolute constant C, and the shorthand A ≪ B to mean A ≤ B/polylog(n).




§ MAIN RESULTS



In the analysis of our algorithms, we consider the spiked Wishart and Wigner models with signal x satisfying the following properties.

For ρ∈ (0,1] and A ≥ 1, a vector x∈ℝ^n is called (ρ,A)-sparse if


  * x_2 = 1 and x_0 = ρ n, and

  * for any i∈ supp(x), 1/A√(ρ n)≤ |x_i| ≤A/√(ρ n).


Here we have used the standard notations supp(x) = {i ∈ [n]  :  x_i  0} and x_0 = | supp(x)|. We assume that ρ (which may depend on n) is chosen so that ρ n is an integer.


A lower bound on |x_i| is essential for exact support recovery, since we cannot hope to distinguish tiny nonzero entries of x from zero entries. The upper bound on |x_i| is a technical condition that is likely not essential, and is only used for recovery in the Wishart model (Theorem <ref>).



In our calculations of the low-degree likelihood ratio, we instead assume the signal x is drawn from the sparse Rademacher distribution, defined as follows.

The sparse Rademacher prior _n^ρ with sparsity ρ∈ (0,1] is the distribution on ℝ^n whereby x∼_n^ρ has i.i.d. entries distributed as

    x_i = {[      +1 / √(ρ n) with probability           ρ / 2,;      -1 / √(ρ n) with probability           ρ / 2,;                0 with probability           1 - ρ. ].

Note that x∼_n^ρ has x_2 → 1 in probability as n →∞.





 §.§ The Wishart Model


We first present our results for the Wishart model. Our algorithms and results for the Wigner model are essentially identical and can be found in Section <ref>.



The spiked Wishart model with parameters n, N ∈_+, β≥ 0, and planted signal x ∈^n is defined as follows.


  * Under _n = _n,N,β, we observe N independent samples y^(1),…,y^(N)∼𝒩(0,I_n + β x x^⊤).

  * Under _n = _n,N, we observe N independent samples y^(1), …, y^(N)∼𝒩(0,I_n).

We will sometimes specify a prior 𝒳_n for x, in which case _n first draws x ∼𝒳_n and then draws y^(1),…,y^(N) as above.




  
Detection. We first consider the detection problem, where the goal is to determine whether the given data {y^(i)} was drawn from _n or _n.









The detection algorithm is motivated by the fact that v^⊤ Yv = 1/N∑_i = 1^N  v,y^(i)^2. Under the planted model _n, y^(i)∼ (0,I_n+β xx^⊤) and thus v,y^(i)∼(0, ℓ+β v,x^2) for any fixed v ∈_n,ℓ; as a result, if v correctly “guesses" ℓ entries of x with correct signs (up to a global flip), then the contribution of v,x^2 to the variance of v,y^(i) will cause v^⊤ Y v to be large.



The runtime of Algorithm <ref> is dominated by exhaustive search over _n,ℓ during Step 3, when we compute T. Since |_n,ℓ| = nℓ 2^ℓ≤ (2n)^ℓ, the runtime is n^O(ℓ). If ℓ = ⌈ n^δ⌉ for a constant δ > 0, then the runtime is n^O(n^δ) = exp(n^δ + o(1)).



Consider the spiked Wishart model with a (ρ,A)-sparse signal x, and let γ = n/N. Let {y^(i)}_i = 1^N be drawn from either _n or _n, and let f_n be the output of Algorithm <ref>. Suppose

    ρ≤min(1,β/A^2)β/25A^2γ1/log n.

Let ℓ be any integer in the interval

    ℓ∈[25 A^4γ/β^2ρ^2 nlog n, min(1,β/A^2)A^2/βρ n],

which is nonempty due to (<ref>).
Then, the total failure probability of Algorithm <ref>  satisfies

    _n[f_n = ] + _n[f_n = ] ≤ 2exp(-β^2/48A^4γℓ^2/ρ^2 n) ≤ 2n^-25ℓ/48,

where the last inequality follows from (<ref>).



Since the runtime is n^O(ℓ), for the best possible runtime we should choose ℓ as small as possible, i.e.,

    ℓ = ⌈25A^4γ/β^2ρ^2 nlog n ⌉.




We are primarily interested in the regime n →∞ with γ = Θ(1), A = Θ(1), ρ = n^-τ for a constant τ∈ (0,1), and either β = n^-α for a constant α > 0, or β = Θ(1) with λ̂β/√(γ) < 1 (in which case α 0). In this case, the requirement (<ref>) reads ρ≤Ω(λ̂^2/log n) (or, in other words, τ > 2 α), which is information-theoretically necessary up to log factors <cit.>. Choosing ℓ as in Remark <ref> yields an algorithm of runtime n^O(1 + λ̂^-2ρ^2 n log n) = poly(n) + exp(n^2α - 2τ + 1 + o(1)).



For ⊆ [n], let Y_ denote the corresponding principal submatrix of Y (i.e., restrict to the rows and columns whose indices lie in ). An alternative detection algorithm would be to threshold the test statistic

    T' max_∈[n]ℓλ_max(Y_),

i.e., the largest eigenvalue of any ℓ×ℓ principal submatrix. One can obtain similar guarantees for this algorithm as for Algorithm <ref>.




  
Recovery. We now turn to the problem of exactly recovering the support and signs of x, given data drawn from _n. The goal is to output a vector x̅∈{-1,0,1}^n such that (x̅) = ±(x) where (x)_i = (x_i) and

    (x_i) = {[            1  if x_i > 0,;           -1 if  x_i < 0,;            0 if  x_i = 0. ].

Note that we can only hope to recover (x) up to a global sign flip, because xx^⊤ = (-x)(-x)^⊤.









For technical reasons, we divide our N samples into two subsamples of size N̅ = ⌊ N/2 ⌋ (with one sample discarded if N is odd) and produce two independent sample covariance matrices Y' and Y”. The first step of the algorithm is similar to the detection algorithm: by exhaustive search, we find the vector v^* ∈_n,ℓ maximizing v^⊤ Y' v. In the course of proving that the algorithm succeeds, we will show that v^* has nontrivial correlation with x. The second step is to recover the support (and signs) of x by thresholding z = (Y”-I)v^*. Note that z discards (i.e., does not depend on) the columns of Y” that do not lie in (v^*); since (v^*) has substantial overlap with (x), this serves to amplify the signal.


Consider the planted spiked Wishart model _n with an arbitrary (ρ,A)-sparse signal x, and let γ = n/N. Suppose

    ρ≤min(1,β/25 A^8)A^4β/400γ1/log n.

Let ℓ be any integer in the interval

    ℓ∈[10000 A^4γ/β^2ρ^2 nlog n, min(1,β/25 A^8)25A^8/βρ n],

which is nonempty due to (<ref>).
Then the failure probability of Algorithm <ref> satisfies 

    1 - _n[ sign(x̅) = ± sign(x)]≤ 6exp(-β^2/6400A^4γℓ/ρ^2 n)≤ 6n^-3/2,

where the last inequality follows from (<ref>).

 As for detection, the runtime of Algorithm <ref> is n^O(ℓ), and we can minimize this by choosing

    ℓ = ⌈10000 A^4γ/β^2ρ^2 nlog n ⌉.



Once we obtain supp(x) using Algorithm <ref>, it is straightforward to estimate x (up to global sign flip) using the leading eigenvector of the appropriate submatrix. This step of the algorithm requires only polynomial time.

Consider the planted spiked Wishart model _n with an arbitrary (ρ,A)-sparse signal x, and let γ = n/N.
Suppose we have access (e.g., via Algorithm <ref>) to =  supp(x) ⊂ [n]. Write P_ = ∑_i∈e_i e_i^⊤, y^(i)_ = P_y^(i) and Y_ = P_YP_^⊤ = 1/N∑_i = 1^Ny^(i)_y^(i)_^⊤. Let x̃ denote the unit-norm eigenvector corresponding to the maximum eigenvalue of Y_. Then, there exists an absolute constant C>0 such that, for any ϵ∈ [2(1+β)√(γρ)/Cβ,1),

    _n [x̃,x^2 ≤ 1-ϵ] ≤ 2exp(-C^2β^2nϵ^2/4(1+β)^2γρ) ≤ 2 exp(-n).



In the regime we are interested in, n →∞ with A = O(1), β = O(1), and (<ref>) is satisfied. In this case, the conclusion of Theorem <ref> gives ⟨x̃,x ⟩^2 > 1 - o(1) with high probability.





  
Low-degree likelihood. Now, we turn to controlling the low-degree likelihood ratio (LDLR) (see Section <ref>) to provide rigorous evidence that the above algorithms are optimal. In this section we take a fully Bayesian approach, and assume that the planted signal x is drawn from the sparse Rademacher prior _n^ρ. Recall that the signal-to-noise ratio is defined as λ̂β/√(γ).

As discussed in Section <ref>, we will determine the behavior of L_n^≤ D in the limit n →∞: if L_n^≤ D = O(1), this suggests hardness for n^Ω̃(D)-time algorithms. We allow the parameters D,ρ,β,γ to depend on n, which we sometimes emphasize by writing, e.g., ρ_n. For D_n = o(n), our results suggest hardness for n^Ω̃(D)-time algorithms whenever λ̂< 1 and ρ≫λ̂√(D_n/n). This is essentially tight, matching PCA (which succeeds when λ̂> 1) and our algorithm with ℓ = D_n (which succeeds when ρ≪λ̂√(D_n/n)) (however, see Remark <ref> below for one caveat).


 Under the spiked Wishart model with spike prior = _n^ρ, suppose D_n = o(n). If one of the following holds for sufficiently large n:


  (a) lim sup_n→∞λ̂_n < 1 and 
    ρ_n≥max(1,√(1/6log(1/λ̂_n)))√(D_n/n), or


  (b) lim sup_n→∞λ̂_n < 1/√(3) and

    ρ_n≥λ̂_n√(D_n/n),


then, as n→∞, L_n,N,β,^≤ D = O(1).


The following result on divergence of the LDLR serves as a sanity check: we show that L_n^≤ D indeed diverges in the regime where we know that a n^Ω̃(D)-time algorithm exists.


Under the spiked Wishart model with spike prior = _n^ρ, suppose D_n = ω(1) and D_n = o(n). If one of the following holds:


  (a) lim inf_n→∞λ̂_n > 1, or

  (b) lim sup_n→∞λ̂_n < 1, |logλ̂_n| = o(√(D_n)) and for sufficiently large n,

    ρ_n < Cλ̂_nlog^-2(1/λ̂_n)√(D_n/n)

where C is an absolute constant,

then, as n→∞, L_n,N,β,^≤ D = ω(1).



There is one regime where the above results give some unexpected behavior. Recall first that optimal Bayesian inference for sparse PCA can be performed in time n^O(ρ n) by computing the likelihood ratio. Thus if L_n^≤ D = O(1) for some D_n ≫ρ n, this suggests that the problem is information-theoretically impossible; from our results above, there are regimes where this occurs (and indeed the problem is information-theoretically impossible), yet L_n^≤ D = ω(1) for some larger D_n (which incorrectly suggests that there should be an algorithm). This is analogous to a phenomenon where the second moment of the (non-low-degree) likelihood ratio L_n can sometimes diverge even when strong detection is impossible (see, e.g. <cit.>). Luckily, this issue never occurs for us in the regime of interest D_n ≪ρ n, and therefore does not prevent our results from being tight. Note also that none of these observations contradict Conjecture <ref>.





 §.§ The Wigner Model



We now state our algorithms and results for the Wigner model. These are very similar to the Wishart case, so we omit some of the discussion.



The spiked Wigner model with parameters n∈_+, λ≥ 0, and planted signal x ∈^n is defined as follows.


  * Under _n = _n,λ, we observe the matrix Y = W + λ xx^⊤, where W∼(n).

  * Under _n, we observe the matrix Y ∼(n).












As in the Wishart case (see Remark <ref>), the runtime is n^O(ℓ). The same holds for Algorithm <ref> below.



Consider the spiked Wigner model with an arbitrary (ρ,A)-sparse signal x. Let Y be drawn from either _n or _n, and let f_n be the output of Algorithm <ref>. Suppose

    ρ≤λ^2/36A^41/log n.

Let ℓ be any integer in the interval

    ℓ∈[36 A^4/λ^2ρ^2 nlog n, ρ n],

which in nonempty due to (<ref>).
Then the total failure probability of Algorithm <ref> satisfies

    _n[f_n = ] + _n[f_n = ] ≤ 2exp(-λ^2/32A^4ℓ^2/ρ^2 n) ≤ 2n^-9 ℓ/8,

where the last inequality follows from (<ref>).



Since our lower bounds are against the class of low-degree algorithms, it is natural to ask whether our algorithms fall into this class. While our test statistic T is not a polynomial function of Y, we can instead take as a proxy the degree-2k polynomial P(Y) = ∑_v ∈_n,ℓ (v^⊤ Y v)^2k for some choice of k. Our analysis can be adapted to show that P can be used to solve strong detection under essentially the same conditions as Theorem <ref>, provided k ≳ℓlog n. Note that (up to log factors) this matches the correspondence between runtime and degree in Conjecture <ref>.










For technical reasons, our first step is to fictitiously “split” the data into two independent copies Y' and Y”. Note that

    Y' = λ/√(2)xx^⊤+ W+W̃/√(2)  and   Y” = λ/√(2)xx^⊤+ W-W̃/√(2).

Since W' W+W̃/√(2) and W”W-W̃/√(2) are independent (n) matrices, Y' and Y” are distributed as independent observations drawn from _n with the same planted signal x and with effective signal-to-noise ratio λ̅ = λ/√(2). 


Consider the planted spiked Wishart model _n with an arbitrary (ρ,A)-sparse signal x. Suppose

    ρ≤λ^2/338A^41/log n.

Let ℓ be any integer in the interval

    ℓ∈[338 A^4/λ^2ρ^2 nlog n, ρ n],

which is nonempty due to (<ref>). Then the failure probability of Algorithm <ref> satisfies 

    1-_n[ supp(x̅) =  supp(x),  sign(x̅) = ± sign(x))≤ 4exp(-λ^2/288A^4ℓ/ρ^2 n)≤ 4n^-169/144,

where the last inequality follows from (<ref>).


As in the Wishart case, once we have recovered the support, there is a standard polynomial-time spectral method to estimate x.

Consider the planted spiked Wigner model _n with an arbitrary (ρ,A)-sparse signal x. Suppose we have access (e.g., via Algorithm <ref>) to =  supp(x) ⊂ [n]. Write P_ = ∑_i∈e_i e_i^⊤ and Y_ = P_YP_^⊤. Let x̃ denote the unit-norm eigenvector corresponding to the maximum eigenvalue of Y_. Then for any ϵ∈ (4√(2ρ)/λ,1), 

    _n[x̃,x^2 ≤ 1-ϵ]
        ≤ 4exp[-n/16(λϵ-4√(2ρ))^2].



In the regime we are interested in, n →∞ with (<ref>) satisfied, so that √(ρ)/λ→ 0. In this case, the conclusion of Theorem <ref> gives ⟨x̃,x ⟩^2 > 1 - o(1) with high probability, upon choosing for example ϵ = 8√(2ρ)/λ.


We also have the following results on the behavior of the low-degree likelihood ratio.

 Under the spiked Wigner model with prior = _n^ρ, suppose D_n = o(n). If one of the following holds for sufficiently large n:


  (a) lim sup_n→∞λ_n < 1 and 
    ρ_n≥max(1,√(1/6log(1/λ_n)))√(D_n/n), or


  (b) lim sup_n→∞λ_n < 1/√(3) and

    ρ_n≥λ_n√(D_n/n),


then, as n→∞, L_n,λ,^≤ D = O(1).


 Under the spiked Wigner model with prior = _n^ρ, suppose D_n = ω(1) and D_n = o(n). If one of the following holds:


  (a) lim inf_n→∞λ_n > 1, or

  (b) lim sup_n→∞λ_n < 1, |logλ_n| = o(√(D_n)) and for sufficiently large n,

    ρ_n < Cλ_nlog^-2(1/λ_n)√(D_n/n)

where C is an absolute constant,

then, as n→∞, L_n,λ,^≤ D = ω(1).




§ PROOFS FOR SUBEXPONENTIAL-TIME ALGORITHMS





 §.§ The Wishart Model


Under ℚ_n, for any fixed v∈_n,ℓ we have v^⊤ y^(i)∼(0,ℓ) for i∈ [N] and

    v^⊤ Yv = 1/N∑_i = 1^N (v^⊤ y^(i))^2 ℓ/Nχ_N^2,

where χ_N^2 is a chi-squared random variable with N degrees of freedom, i.e., the sum of the squares of N standard gaussians.
Using Corollary <ref>, we union bound over v∈_n,ℓ for any t∈ (0,1/2):

    ℚ_n[T≥ℓ(1+t)]    ≤ | _n,ℓ|[χ_N^2 ≥N(1+t)]
       ≤nℓ2^ℓ·exp(-Nt^2/3)
       ≤exp(ℓlog(2n) -Nt^2/3).

Under the condition 

    Nt^2 ≥ 4 ℓlog (2n)

we have

    ℚ_n[T≥ℓ(1+t)] ≤exp(-1/12Nt^2).

Meanwhile, under ℙ_n, when v = v̅ correctly guesses ℓ entries and their signs in the support of x (which requires ℓ≤ρ n), for any i∈ [N] we have

    v̅^⊤ y^(i)∼(0,v̅^⊤ (I_n+β xx^⊤)v̅) = (0,ℓ+βv̅,x^2).

Therefore,

    v̅^⊤ Yv̅ = 1/N∑_i = 1^N (v̅^⊤ y^(i))^2 1/N(ℓ+βv̅,x^2)χ_N^2

where v̅,x^2 ≥ℓ^2/A^2ρ n.
As a result, by Corollary <ref>,

    ℙ_n[T<ℓ(1+t)]    ≤ℙ_n[v̅^⊤ Yv̅ < ℓ(1+t)]
       = [1/N(ℓ+βv̅,x^2)χ_N^2 < ℓ(1+t)]
       = [χ_N^2 < N(1-βℓ-A^2ρ n t/βℓ +A^2ρ n)]
       ≤exp(-N/3(βℓ-A^2ρ n t/βℓ +A^2ρ n)^2),

the last inequality requiring

    0 ≤βℓ-A^2ρ nt/βℓ +A^2ρ n≤1/2.

To satisfy t∈ (0,1/2), (<ref>) and (<ref>) at the same time, we choose 

    t = βℓ/2A^2ρ n.

Under the condition

    βℓ/A^2 n≤ρ≤β/5A^2√(γ)√(ℓ/nlog n),

which is equivalent to the interval for ℓ given in (<ref>), thresholding the statistic T at ℓ(1+t) succeeds at distinguishing _n and _n with total error probability 

    ℚ_n[T≥ℓ(1+t)]
    +ℙ_n[T<ℓ(1+t)] ≤ 2exp(-β^2/48A^4γℓ^2/ρ^2 n),

which completes the proof.



First, we give a high-probability lower bound on v^*,x. From the analysis of the detection algorithm, we know that under the condition (<ref>),

    1-2exp(-β^2/48A^4γℓ^2/ρ^2 n) 
       ≤_n [ℓ(1+βℓ/2A^2ρ n) ≤v^*^⊤ Y' v^*] 
       ≤_n[ℓ(1+βℓ/2A^2ρ n) ≤1/N̅(ℓ+β v^*,x^2)χ_N̅^2,  v^*,x^2 < ℓ^2/3A^2ρ n] 
       +_n[ℓ(1+βℓ/2A^2ρ n) ≤1/N̅(ℓ+β v^*,x^2)χ_N̅^2,  v^*,x^2 ≥ℓ^2/3A^2ρ n] 
       ≤_n[ℓ(1+βℓ/2A^2ρ n) ≤1/N̅(ℓ+βℓ^2/3A^2ρ n)χ_N̅^2] + _n[ v^*,x^2 ≥ℓ^2/3A^2ρ n],

where

    _n[ℓ(1+βℓ/2A^2ρ n) ≤1/N̅(ℓ+βℓ^2/3A^2ρ n)χ_N̅^2]    =
    _n[χ_N̅^2 ≥N̅(1+βℓ/2βℓ +6A^2ρ n)]
       ≤exp(-N̅/3(βℓ/2βℓ+6A^2ρ n)^2)
       ≤exp(-β^2/384A^4γℓ^2/ρ^2 n),

hence we have the lower bound

    _n[ v^*,x^2 ≥ℓ^2/3A^2ρ n]    ≥ 1-2exp(-β^2/48A^4γℓ^2/ρ^2 n) -  exp(-β^2/384A^4γℓ^2/ρ^2 n)
       ≥ 1-3exp(-β^2/384A^4γℓ^2/ρ^2 n).

We now fix v^* satisfying the above lower bound on v^*,x ^2. From this point onward, we will only use the second copy Y” of our data; note that, crucially, Y” is independent from v^*. To simplify the notation, we will write y^(1),…,y^(N̅) instead of y^(N̅+1),…,y^(2N̅) for the samples used to form Y”. We now adopt an equivalent representation of the observations: y^(i) = u^(i) + √(β) w^(i) x, where u^(i)∼ (0,I_n) and w^(i)∼ (0,1) are independent random gaussian vectors and scalars, respectively. Substituting this into z = (Y”-I)v^* yields
 
    z_j = 1/N̅∑_i = 1^N̅(a_ij+b_ij+c_ij +d_ij +e_ij)

where, for i ∈ [N̅] and j∈ [n],

    a_ij   = (w^(i))^2 β x_j v^*,x
    
    b_ij   = ((u^(i)_j)^2-1)v^*_j
    
    c_ij   = ∑_k≠ ju^(i)_j u^(i)_k v^*_k
    
    d_ij   = u^(i)_j w^(i)√(β) v^*,x
    
    e_ij   =  u^(i), v^*  w^(i)√(β) x_j,

with 𝔼(a_ij) = β x_j v^*,x and 𝔼(b_ij) = 𝔼(c_ij) = 𝔼(d_ij) = 𝔼(e_ij) = 0. We will show separate union bounds for these five contributions to z_j. In the following, we fix the constant μ = 1/20.



  
Union bound for a_ij.
For all j∈ supp(x),

    βℓ/√(3) A^2ρ n≤β1/A√(ρ n)ℓ/√(3)A√(ρ n)≤ |β x_j v^*,x| ≤βA/√(ρ n)Aℓ/√(ρ n) = A^2βℓ/ρ n,

so by Corollary <ref>,

    log_n[|1/N̅∑_i = 1^N̅ a_ij-β x_j v^*,x| > μβℓ/A^2ρ n]    = 
    log[|1/N̅χ_N̅^2 -1| > μβℓ/A^2ρ n·1/|β x_j v^*,x|]
       ≤log[|1/N̅χ_N̅^2 -1| > μ/A^4]
       ≤ -μ^2 N̅/3A^8.

Therefore, we may union bound over j∈ supp(x):

    _n[|1/N̅∑_i = 1^N̅ a_ij-β x_j v^*,x| ≤μβℓ/A^2ρ n,  for all  j∈ [n]]
       ≥ 1-∑_j∈ supp(x)_n[|1/N̅∑_i = 1^N̅ a_ij-β x_j v^*,x| > μβℓ/A^2ρ n]
       ≥ 1-exp(log n-μ^2 N̅/3A^8)
       ≥ 1-exp(-μ^2 n/7A^8γ).




  
Union bound for b_ij.
We have b_ij nonzero only when j∈ supp(v^*). For such j, by Corollary <ref>,

    log_n[|1/N̅∑_i = 1^N̅ b_ij| > μβℓ/A^2ρ n]    = 
    log[|1/N̅χ_N̅^2 -1| > μβℓ/A^2ρ n]
       ≤ -N̅/3(μβℓ/A^2ρ n)^2.

 Therefore,

    _n[|1/N̅∑_i = 1^N̅ b_ij| ≤μβℓ/A^2ρ n,  for all  j∈ [n]]    ≥ 1-∑_j∈ supp(v^*)_n[|1/N̅∑_i = 1^N̅ b_ij| > μβℓ/A^2ρ n]
       ≥ 1-exp(logℓ -N̅/3(μβℓ/A^2ρ n)^2)
       ≥ 1-exp(-μ^2β^2/12A^4γℓ^2/ρ^2 n),

under the condition

    ρ≤1/√(12)μβ/A^2√(γ)√(ℓ/nlogℓ).




  
Union bound for c_ij and d_ij.
In the following, let u, u' denote independent samples from (0,I_N̅). Note that 
    ũ_j^(i)∑_k≠ j u^(i)_k v^*_k ∼(0,ℓ̃_j),  where ℓ̃_j = ℓ-{j∈ supp(v^*)},

and ũ_j^(i) is independent from u^(i)_j. 
Therefore,

    ∑_i = 1^N̅ c_ij = ∑_i = 1^N̅ u^(i)_j ũ_j^(i)√(ℓ̃_j) u,u'.

Therefore, by Lemma <ref>, for the c_ij we have

    _n[|1/N̅∑_i = 1^N̅ c_ij| ≤μβℓ/A^2ρ n,  for all  j∈ [n]]    ≥ 1-∑_j = 1^n[| u,u'| > μβ√(ℓ)N̅/A^2ρ n]
       ≥ 1-2nexp(-1/4N̅(μβ√(ℓ)N̅/A^2ρ n)^2)
       ≥ 1-exp(-μ^2β^2/16A^4γℓ/ρ^2 n).

The last inequality holds under the condition

    μβ√(ℓ)/A^2ρ n≤1/2, μ^2β^2/16A^4γℓ/ρ^2 n≥log (2n),
 which follows from 
    2μβ√(ℓ)/A^2 n≤ρ≤μβ/5A^2√(γ)√(ℓ/nlog n).


Meanwhile,

    ∑_i = 1^N̅ d_ij√(β) v^*,x u,u'.

Therefore, as for the c_ij, for the d_ij we have

    _n[|1/N̅∑_i = 1^N̅ d_ij| ≤μβℓ/A^2ρ n,  for all  j∈ [n]] 
       ≥ 1-∑_j = 1^n[| u,u'| > μ√(β)N̅/A^3√(ρ n)]
       ≥ 1-2nexp(-1/4N̅(μ√(β)N̅/A^3√(ρ n))^2)
       ≥ 1-exp(-μ^2β/16A^6γ1/ρ).

The last inequality holds under the condition

    μ√(β)/A^3√(ρ n)≤1/2, μ^2β/16A^6γ1/ρ≥log(2n),
 which follows from 
    4μ^2β/A^6 n≤ρ≤μ^2β/17A^6γlog n.




  
Union bound for e_ij.
We have

    ∑_i = 1^N̅ e_ij x_j√(βℓ) u,u' ,

which is only nonzero for j∈ supp(x). Therefore,

    _n[|1/N̅∑_i = 1^N̅ e_ij| ≤μβℓ/A^2ρ n, for all  j∈ [n]]    ≥ 1-∑_j∈ supp(x)[| u,u'| > μ√(βℓ)N̅/A^2ρ n |x_j|]
       ≥ 1-∑_j∈ supp(x)[| u,u'| > μ√(βℓ)N̅/A^4√(ρ n)]
       ≥ 1-2nexp(-1/4N̅(μ√(βℓ)N̅/A^4√(ρ n))^2)
       ≥ 1-exp(-μ^2β/16A^8γℓ/ρ).

The last inequality holds under the condition

    μ√(βℓ)/A^4√(ρ n)≤1/2, μ^2β/16A^8γℓ/ρ≥log (2n),

which follows from

    4μ^2βℓ/A^8 n≤ρ≤μ^2βℓ/17A^8γlog n.




  
Final steps.
Now, combining all of the union bounds and conditions from (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>), assuming that β, γ = Θ(1) and that ω(1) ≤ℓ(n) ≤ o(n/log n), under the condition

    max(1,β/25A^8)ℓ/n≤ρ≤β/100A^2√(γ)√(ℓ/nlog n)

which is equivalent to the regime for ℓ given in (<ref>) that we are considering, we have

    _n[for some  j, |z_j-β x_j v^*,x| > βℓ/4A^2ρ n] + _n[ v^*,x^2 ≤ℓ^2/3A^2ρ n] 
       ≤_n[for some  j, |1/N̅∑_i = 1^N̅ a_ij-β x_j v^*,x| > μβℓ/A^2ρ n] 
        +_n[for some  j, |1/N̅∑_i = 1^N̅ b_ij| > μβℓ/A^2ρ n]
        +_n[for some  j, |1/N̅∑_i = 1^N̅ c_ij| > μβℓ/A^2ρ n] 
        +_n[for some  j, |1/N̅∑_i = 1^N̅ d_ij| > μβℓ/A^2ρ n]
       +_n[for some  j, |1/N̅∑_i = 1^N̅ e_ij| > μβℓ/A^2ρ n] 
        +_n[ v^*,x^2 ≤ℓ^2/3A^2ρ n]
       ≤ 5exp(-β^2/6400A^4γℓ/ρ^2 n) +exp(-β^2/384A^4γℓ^2/ρ^2 n) 
       ≤ 6exp(-β^2/6400A^4γℓ/ρ^2 n).

Since |β x_j v^*,x| ≥βℓ/√(3)A^2ρ n for j∈ supp(x) and |β x_j v^*,x| = 0 for j∉ supp(x), we conclude that, with probability at least 1-6exp(-β^2/6400A^4γℓ/ρ^2 n), for every j∈ [n], 

    j∈ supp(x)  if and only if  |z_j|≥βℓ/2√(3)A^2ρ n,

and

    sign(z_j) =  sign(x_j v^*,x),

completing the proof.




By a result in the analysis of covariance matrix estimation for subgaussian distributions (<cit.>, Remark 5.51), there exists an absolute constant C>0 such that, for any δ∈ [√(ργ)/C,1), the following holds with probability at least
1-2exp(-C^2δ^2 N/ρ):

    Y_-(P_+β xx^⊤)≤δP_+β xx^⊤ = δ(1+β).

Whenever this is true, by the definition of spectral norm we have

    δ(1+β)    ≥x̃^⊤[Y_-(P_+β xx^⊤)]x̃
       = Y_-(1+βx̃,x ^2)
       ≥ (1-δ)(1+β)-(1+βx̃,x ^2),

which is equivalent to x̃,x^2 ≥ 1-ϵ upon taking δ = βϵ/2(1+β). Thus, for any ϵ∈ (2(1+β)√(ργ)/Cβ,1),

    [x̃,x^2 ≤ 1-ϵ] ≤ 2exp(-C^2β^2nϵ^2/4(1+β)^2γρ),

which completes the proof.




 §.§ The Wigner Model


For simplicity we denote t = λℓ^2/2A^2ρ n. Under ℙ_n, when v̅ correctly guesses ℓ entries in the support of x with correct signs (which requires ℓ≤ρ n),

    v̅^⊤ Y v̅ = v̅^⊤ Wv̅ + λv̅,x^2,

where v̅^⊤ Wv̅∼𝒩(0,ℓ^2/n).
Note that 

    λv̅,x^2 ≥λℓ^2/A^2ρ n=2t.

Therefore, a standard Gaussian tail bound gives

    _n[T < t]    ≤_n[v̅^⊤ Yv̅ < t] 
       ≤[𝒩(0,ℓ^2/n) > t] 
       ≤exp(-n/2ℓ^2(λℓ^2/2A^2ρ n)^2) 
       = exp(-λ^2/8A^4ℓ^2/ρ^2 n).

Under ℚ_n, for each fixed v ∈ℐ_n,ℓ, we have

    v^⊤ Y v ∼𝒩(0,2ℓ^2/n).

By the same tail bound,

    _n[v^⊤ Y v ≥ t] ≤exp(-nt^2/4ℓ^2).

Now, by a union bound over v ∈ℐ_n,ℓ,

    _n[T ≥ t]    ≤ |ℐ_n,ℓ| exp(-nt^2/4ℓ^2)
       = nℓ 2^ℓexp(-nt^2/4ℓ^2)
       ≤exp(ℓlog(2n) -nt^2/4ℓ^2).

Under the condition

    nt^2/8ℓ^2≥ℓlog(2n) ⇐ ρ < λ/6A^2√(ℓ/nlog n),

which is equivalent to the interval for ℓ given in (<ref>), we have

    _n[T ≥ t] ≤exp(-nt^2/8ℓ^2) = exp(-λ^2/32A^4ℓ^2/ρ^2 n).

Therefore, by thresholding T at t, under the condition 

    ℓ/n≤ρ≤λ/6A^2√(ℓ/n log n),

we can distinguish _n and _n with total failure probability at most

    _n[T < t] + _n[T ≥ t] ≤exp(-λ^2/8A^4ℓ^2/ρ^2 n)+exp(-λ^2/32A^4ℓ^2/ρ^2 n)≤ 2exp(-λ^2/32A^4ℓ^2/ρ^2 n),

completing the proof.




First, we show that v^* has significant overlap with the support of x. From the analysis of the detection algorithm, provided (<ref>) holds, with probability at least 1-2exp(-λ̅^2/32A^4ℓ^2/ρ^2 n) we have

    λ̅ℓ^2/2A^2ρ n≤v^*^⊤ Y' v^* = λ̅⟨ v^*,x ⟩^2 + v^*^⊤ W' v^*.

where v^*^⊤ W' v^* ∼(0,2ℓ^2/n).
Therefore, for n sufficiently large, 

    _n[⟨ v^*,x ⟩^2 ≥ℓ^2/4A^2ρ n]    ≥(1-2exp(-λ̅^2/32A^4ℓ^2/ρ^2 n)) (1- [(0,2ℓ^2/n) ≥λ̅ℓ^2/4A^2ρ n])
       ≥ 1-2exp(-λ̅^2/32A^4ℓ^2/ρ^2 n)-exp(-λ̅^2/64A^4ℓ^2/ρ^2 n)
       ≥ 1-3exp(-λ̅^2/64A^4ℓ^2/ρ^2 n).

We now fix v^* satisfying the above lower bound on ⟨ v^*,x ⟩^2. From this point onward, we will only use the second copy Y” of our data; it is important here that Y” is independent from v^*. We will that x is successfully recovered by thresholding the entries of z = Y”v^*. Entrywise, we have

    z_i = λ̅ x_i ⟨ v^*,x ⟩ + e_i^⊤ W”v^*.

For all i ∈supp(x),

    |λ̅ x_i ⟨ v^*,x ⟩| ≥λ̅1/A√(ρ n)·ℓ/2A√(ρ n) = λ̅ℓ/2A^2ρ n.

For simplicity we denote s = λ̅ℓ/2A^2ρ n and μ = 1/3. Note that for all i∈ [n], e_i^⊤ W” v^* ∼𝒩(0,v^2/n) = 𝒩(0,ℓ/n) and therefore

    _n[|e_i^⊤ W” v^*| ≥μ s] ≤ 2 exp(-nμ^2 s^2/2ℓ).

By a union bound over all i ∈ [n],

    _n[|e_i^⊤ W” v^*| ≤μ s  for all  i]    ≥ 1 - 2n exp(-nμ^2 s^2/2ℓ) 
       ≥ 1 - exp(log (2n) -nμ^2 s^2/2ℓ)
       ≥ 1 - exp(-nμ^2 s^2/4ℓ) 
       = 1-exp(-λ̅^2/144A^4ℓ/ρ^2 n)

under the condition

    nμ^2 s^2/4ℓ≥log(2n) ⇐ ρ≤λ̅/13A^2√(ℓ/n log n) = λ/13√(2)A^2√(ℓ/n log n),

which, combined with (<ref>), is equivalent membership in the interval for ℓ that we are considering per (<ref>).
Therefore, with probability at least

    1-3exp(-λ̅^2/64A^4ℓ^2/ρ^2 n)-exp(-λ̅^2/144A^4ℓ/ρ^2 n) ≥ 1-4exp(-λ^2/288A^4ℓ/ρ^2 n)

for all j∈ [n],

    j∈ supp(x)  if and only if  |z_j|≥s/2

and

    sign(z_j) =  sign(x_j  v^*,x).

Thus, we find that thresholding the entries of z at s/2 successfully recovers the support and signs of x, completing the proof.



Since Y_x̃ = λ_max(Y_)x̃, we must have supp(x̃)⊂. Denote W_ = P_WP_^⊤ and W̅_ the ℓ×ℓ submatrix of W_ with rows and columns indexed by  (the only nonzero rows and columns). Now, the variational description of the leading eigenvector yields

    x̃^⊤ W_x̃+λx̃,x^2 = x̃^⊤ Y_x̃≥ x^⊤ Y_ x
    = x^⊤ W_x +λ.

Therefore,

    x̃,x^2 ≥ 1- 1/λ(x̃^⊤ W_x̃-x^⊤ W_ x)
    ≥ 1-1/λ(λ_max(W̅_)+λ_max(-W̅_)).

Note that W̅_ has the same law as (G̅+G̅^⊤) / √(2n), where G̅ is an ρ n×ρ n matrix whose entries are independent standard normal random variables. Now, for any ϵ > 4√(2ρ)/λ, we have √(2n)λϵ/4 > 2√(ρ n), a standard singular value estimate for Gaussian matrices (see <cit.>, Corollary 5.35) gives

    _n[x̅,x^2 ≤ 1-ϵ]
           ≤[λ_max(W̅_)+λ_max(-W̅_) ≥λϵ]
       ≤ 2[λ_max(W̅_) ≥λϵ/2]
       ≤ 2[σ_max(G̅) ≥√(2n)λϵ/4]
       ≤ 4exp[-n/16(λϵ-4√(2ρ))^2],

which conclues the proof.




§ PROOFS FOR LOW-DEGREE LIKELIHOOD RATIO BOUNDS





 §.§ Low-Degree Likelihood Ratio for Spiked Models


We begin by giving expressions for the norm of the low-degree likelihood ratio (LDLR) for the spiked Wigner and Wishart models. These expressions are derived in <cit.> and <cit.>, respectively.


    Let L_n,λ,^≤ D denote the degree-D likelihood ratio for the spiked Wigner model with parameters n,λ and spike prior . Then,
    
    L_n,λ,^≤ D^2 = _v^(1), v^(2)∼_n[∑_d = 0^D 1/d!(n/2λ^2 v^(1), v^(2)^2)^d]

where v^(1),v^(2) are drawn independently from _n.



Let L_n,N,β,^≤ D denote the degree-D likelihood ratio for the spiked Wishart model with parameters n,N,β and spike prior . Define

    φ_N(x)     (1-4x)^-N/2
    φ_N,k(x)    ∑_d = 0^k x^d∑_d_1,…,d_N
    ∑ d_i = d∏_i = 1^N2d_id_i,

so that φ_N,k(x) is the Taylor series of φ_N around x = 0 truncated to degree k. Then,

    L_n,N,β,^≤ D_2^2    = _v^(1),v^(2)∼_n[φ_N,⌊ D/2 ⌋(β^2  v^(1),v^(2)^2/4)]
       = _v^(1),v^(2)∼_n∑_d = 0^⌊ D/2 ⌋(∑_d_1,…,d_N
    ∑ d_i = d∏_i = 1^N2d_id_i)(β^2 v^(1),v^(2)^2/4)^d,

where v^(1),v^(2) are drawn independently from _n.

We consider a signal x drawn from the sparse Rademacher prior, _n = _n^ρ. The goal of this section is to prove upper and lower bounds on the LDLR expressions in (<ref>) and (<ref>) as n→∞, for certain regimes of the parameters (λ,ρ for the Wigner model and β,γ,ρ for the Wishart model).
These bounds are obtained in several steps. First, we treat the moment terms

    A_d  (nρ)^2d_v^(1),v^(2)∼_n^ρ v^(1),v^(2)^2d

from (<ref>) and (<ref>) in Section <ref>, with upper bounds given in Lemmas <ref> and <ref> and a lower bound given in Lemma <ref>. We then give a precise estimate in Lemma <ref> of the coefficient 
    ∑_d_1,…,d_N
    ∑ d_i = d∏_i = 1^N2d_id_i

in the LDLR (<ref>) of the Wishart model. Finally, by combining the above bounds, we show regimes of parameters under which the LDLR either remains bounded or diverges as n →∞. This yields the proofs of Theorems <ref> and <ref> for the Wishart model, and Theorems <ref> and <ref> for the Wigner model.



 §.§ Introduction and Estimates of A_d

In this section, we carry out combinatorial estimates of the moments A_d defined in (<ref>), which appear in the LDLR expressions (<ref>) and (<ref>). We give upper bounds (Lemmas <ref> and <ref>) and a lower bound (Lemma <ref>) on these moments.

For independent v^(1),v^(2) drawn from the sparse Rademacher prior _n^ρ, v^(1),v^(2) has the same distribution as S_n,ρ = 1/n∑_i = 1^n s_i,ρ for i.i.d. s_i,ρ with

    s_i,ρ = {[           +1 / ρ with probability         ρ^2 / 2,;           -1 / ρ with probability         ρ^2 / 2,;                0 with probability         1 - ρ^2, ].

and kth moment (for k>0) given by

    _s_i, ρ^k = {[             0  for  k  odd,;       ρ^2 - k for  k  even. ].

Therefore, the moments of v^(1),v^(2) have the combinatorial description

    _v^(1),v^(2)∼_n^ρ v^(1),v^(2)^2d   = n^-2d𝔼S_n,ρ^2d
       = n^-2d∑_i_1,…,i_2d∈ [n]𝔼s_i_1,ρs_i_2,ρ… s_i_2d,ρ
       = n^-2d∑_a_1, …, a_n ≥ 0 
    ∑ a_i = d2d2a_1 ⋯ 2a_n 𝔼∏_a_j > 0s_j,ρ^2a_j
       = n^-2d∑_a_1, …, a_n ≥ 0 
    ∑ a_i = d2d2a_1 ⋯ 2a_n ∏_a_j > 0ρ^2-2a_j
       = n^-2d∑_a_1, …, a_n ≥ 0 
    ∑ a_i = d2d2a_1 ⋯ 2a_n ρ^2|{ i: a_i > 0 }|-2d.


Recall, from (<ref>), that

    A_d = (nρ)^2d_v^(1),v^(2)∼_n^ρ v^(1),v^(2)^2d = ∑_a_1, …, a_n ≥ 0 
    ∑ a_i = d2d2a_1 ⋯ 2a_n ρ^2|{i: a_i > 0}|.

Fix d≤ D_n, and let 1≤ k≤ d be the number of positive numbers among the {a_i}. Suppose d = w k+r, where 0≤ r< k. For positive integers b_1,b_2,…,b_k such that ∑ b_i = d, we claim that

    2d2b_1 ⋯ 2b_k≤2d2w ⋯ 2w  2(w+1) ⋯ 2(w+1) M(k),

and that equality holds if and only if {b_i} consists of r copies of (w+1) and k-r copies of w. This follows from the simple fact that, for any 1≤ i,j≤ k such that b_i≥ b_j+2, we have (2b_i)!(2b_j)! > (2(b_i-1))!(2(b_j+1))!, and therefore the left-hand side of the above inequality becomes strictly larger as we “unbalance” the b_i by replacing b_i and b_j with b_i -1 and b_j+1. As a result, the left-hand side is maximized if and only if the maximum and minimum of {b_1,b_2,…,b_k} differ by at most 1. Now, since the total number of positive integer solutions to ∑_i = 1^k b_i = d is d-1k-1, we have

    A_d  ≤∑_k = 1^d nkd-1k-1ρ^2kM(k).


Before proceeding to bounds on the A_d, we introduce the following result, which will be useful in several estimates in this section.
 
Suppose D_n = o(n). Then, for sufficiently large n,

    n(n-1)⋯ (n-D_n+1)/n^D_n >1/2 e^-D_n^2/n.



By Stirling's formula, as n→∞ and n-D_n→∞ (which is ensured by D_n = o(n)),

    n(n-1)⋯ (n-D_n+1)/n^D_n = n!/(n-D_n)! n^D_n   ∼√(2π n)(n/e)^n/√(2π (n-D_n))(n-D_n/e)^n-D_nn^D_n
       ∼1/e^D_n(1+D_n/n-D_n)^n-D_n.

Here the relation ∼ means that the quotient of the quantities on either side tends to 1 as n→∞. Since D_n = o(n), for large enough n such that D_n < n/3, we have

    log[1/e^D_n(1+D_n/n-D_n)^n-D_n]    = -D_n+(n-D_n)log(1+D_n/n-D_n)
       ≥ -D_n+(n-D_n)(D_n/n-D_n -D_n^2/2(n-D_n)^2)
       ≥ -D_n^2/n,

and the lemma follows.



In the setting of the spiked Wishart or Wigner model with sparse Rademacher prior _n^ρ, suppose D_n = o(n). If for some μ > 0 we have

    ρ≥max(1,√(1/6μ))√(D_n/n),

then for sufficiently large n and for any 1≤ d≤ D_n, A_d defined by (<ref>) satisfies

    A_d ≤∑_k = 1^d G(k) ≤ 2d e^μ d+d^2/n G(d) = 2d e^μ d+d^2/nnd(2d)!/2^dρ^2d,

where

    G(k) nkd-1k-1ρ^2kM(k).



Fix 1≤ d≤ D_n. Recall that the first inequality in (<ref>) is a restatement of (<ref>). By a simple comparison argument, we observe that M(k) is monotone increasing with respect to k. For any 1≤ k < d/2, we have

    G(k+1)/G(k) = nk+1d-1kρ^2k+2M(k+1)/nkd-1k-1ρ^2kM(k)≥(n-k)(d-k)/k(k+1)ρ^2 ≥(n-d)(d/2)/k(k+1)d/n > 1

if ρ≥√(D_n/n). Therefore, 
    ∑_k <d/2G(k) ≤⌊d/2⌋∑_k≥d/2G(k).
 Meanwhile, for d/2≤ k < d, by (<ref>),

    M(k) = max_b_1, …, b_k > 0 
    ∑ b_i = d2d2b_1 ⋯ 2b_k = 2d2 ⋯ 2  4 ⋯ 4 = (2d)!/24^d-k2^2k-d,

since the maximum is attained when {b_i} has (d-k) occurrences of 2 and (2k-d) occurrences of 1. As a result, with the help of Lemma <ref>,

    G(k)/G(d)   = nkd-1k-1ρ^2kM(k)/ndρ^2dM(d)
       ≤1/(6ρ^2)^d-k·nkdk/nd
       ≤1/(6ρ^2)^d-k·d!/(n-k)(n-k-1)… (n-d+1)k!·d!/k!(d-k)!
       ≤1/(6ρ^2)^d-k·2e^d^2/n/n^d-k·(d!/k!)^21/(d-k)!
       ≤ 2e^d^2/n(d^2/6ρ^2 n)^d-k·1/(d-k)!.

Thus for ρ≥√(D/6μ n)≥√(d/6μ n)
we have d^2/6ρ^2 n≤μ d and

    ∑_k = 1^d G(k) ≤(⌊ d/2⌋+1)∑_k≥d/2G(k) ≤ 2de^d^2/n∑_k≥d/2(μ d)^d-k/(d-k)!G(d)≤ 2d e^μ d+d^2/n G(d),

completing the proof.




In the setting of the spiked Wishart or Wigner model with sparse Rademacher prior _n^ρ, suppose D_n = o(n). If for some μ < 1/√(3) we have

    ρ≥μ√(D_n/n),

then for sufficiently large n and for any 11≤ d≤ D_n, 

    A_d ≤∑_k = 1^d G(k) ≲√(d)e^d^2/n(11e/30)^d/2μ^-2d G(d) = √(d)e^d^2/n(11e/30)^d/2μ^-2dnd(2d)!/2^dρ^2d,

where A_d is defined in (<ref>) and G(k) is defined as in Lemma <ref>.


As in the proof of Lemma <ref>, for sufficiently large n and for any 1≤ d≤ D_n and any 1≤ k <d/2, we have

    G(k+1)/G(k)≥(n-k)(d-k)/k(k+1)ρ^2 ≥(n-k)(d-k)/k(k+1)(μ^2d/n) ≥μ^2.

Therefore,

    ∑_k < d/2G(k) ≤ (μ^-2+μ^-4+…+μ^-2(⌈ d/2 ⌉ -1))G(⌈ d/2 ⌉) ≲μ^-2⌈ d/2 ⌉∑_k≥⌈ d/2 ⌉G(k).

Meanwhile, for d/2≤ k < d,

    G(k)/G(d)≤ 2e^d^2/n(d^2/6ρ^2 n)^d-k·1/(d-k)!≤ 2e^d^2/n(d/6μ^2)^d-k·1/(d-k)!.

Summing these quantities, we find (the last inequality requiring d≥ 11)

    ∑_k≥ d/2G(k)/G(d)   ≤ 2e^d^2/n∑_k ≤ d/2(d/6μ^2)^k·1/k!
       ≲ de^d^2/n(d/6μ^2)^⌊ d/2 ⌋·1/⌊ d/2 ⌋!
       ≲√(d)e^d^2/n(ed/6μ^2⌊ d/2 ⌋)^⌊ d/2 ⌋
       ≲√(d)e^d^2/n(11e/30)^⌊ d/2 ⌋μ^-2⌊ d/2 ⌋.

Here we have used the fact that (d / 6μ^2)^k / k! is monotone increasing for 1≤ k ≤d/2, since d / 6μ^2 > d / 2. Combining the two cases, we conclude that

    ∑_k = 1^dG(k)/G(d)≲(1+μ^-2⌈ d/2 ⌉)∑_k ≥ d/2G(k)/G(d)≲√(d)e^d^2/n(11e/30)^d/2μ^-2d,

completing the proof.




In the settings of the spiked Wishart or Wigner models with sparse Rademacher prior _n^ρ, consider a series d = d_n = o(n) with integers w = w_n satisfying w_n | d_n. Then, as d/w→∞, A_d defined by (<ref>) satisfies

    A_d ≥ (1-o(1))nd(2d)!√(w)/2^d [2(d/neρ^2)^1-1/w(w/(2w)!)^1/w]^dρ^2d.



To obtain a lower bound on A_d, we only consider the contribution to the sum from terms {a_i} with d/w-many occurrences of w and (n-d/w)-many occurrences of zero:

    ρ^-2dA_d = ∑_a_1, …, a_n ≥ 0 
    ∑ a_i = d2d2a_1 ⋯ 2a_n ρ^2|{i: a_i > 0}|-2d≥nd/w(2d)!/[(2w)!]^d/wρ^2d/w-2d T_n,d(w).

Now, we calculate the ratio

    T_n,d(w)/T_n,d(1)   = d!(n-d)!/(d/w)!(n-d/w)!2^d/[(2w)!]^d/wρ^2d/w-2d
       = d!/(d/w)!/(n-d/w)(n-d/w-1)⋯ (n-d+1)[2/[(2w)!]^1/wρ^2(1-1/w)]^d.

By Stirling's formula, for w fixed and d sufficiently large,

    d!/(d/w)! = (1+o(1))√(w)(d/e)^d-d/ww^d/w.

Meanwhile,

    (n-d/w)(n-d/w-1)⋯ (n-d+1) ≤ n^d-d/w.

Plugging into (<ref>) we get

    T_n,d(w)/T_n,d(1)≥ (1-o(1))√(w) [2(d/neρ^2)^1-1/w(w/(2w)!)^1/w]^d,

completing the proof.




 §.§ The Wishart Model

In this section, we first carry out an estimate on the extra coefficient occurring in the Wishart LDLR (<ref>) (Lemma <ref>), then use the bounds on A_d (Lemmas <ref>, <ref> and <ref>) to prove the upper bound (Theorem <ref>) and the lower bound (Theorem <ref>) on (<ref>).

Suppose D_n = o(N). There exist absolute constants c_1,c_2 > 0 such that, for sufficiently large N, for any 1≤ d≤ D_n,

    (2N)^d/d!≤∑_d_1,…,d_N
    ∑ d_i = d∏_i = 1^N2d_id_i≤ c_1 d^3/2e^c_2d^2/N(2N)^d/d!.



For the lower bound, note that for any d_i ≥ 2,

    2d_id_i = 2d_i(2d_i-1)/d_i^22(d_i-1)d_i-1≥212(d_i-1)d_i-1,

so for any d_1,…,d_N ≥ 0 such that ∑ d_i = d,

    ∏_i = 1^N2d_id_i≥ 2^d.

Summing over all of the {d_i}, we find

    ∑_d_1,…,d_N
    ∑ d_i = d∏_i = 1^N2d_id_i≥N+d-1d2^d ≥(2N)^d/d!.

For the upper bound, we separately consider those {d_i}'s with exactly k positive entries for each k = 1,2,…,d:

    ∑_d_1,…,d_N
    ∑ d_i = d∏_i = 1^N2d_id_i   = ∑_k = 1^d Nk∑_c_1,…,c_k > 0
    ∑ c_i = d∏_i = 1^k2c_ic_i
       ≤∑_k = 1^d Nkd-1k-1max_c_1,…,c_k > 0
    ∑ c_i = d∏_i = 1^k2c_ic_i.

Given any positive integers c_1,…,c_k such that ∑ c_i = d, if there are two entries c_j≥ c_ℓ≥ 2, consider c̃_j = c_j+1, c̃_ℓ = c_ℓ-1 and c̃_i = c_i for all i≠ j,ℓ. We have the comparison

    ∏_i = 1^k2c̃_ic̃_i/∏_i = 1^k2c_ic_i = 4-2/c_j+1/4-2/c_l+1 > 1.

Therefore, for fixed k, the product ∏_i = 1^k2c_ic_i is maximized when {c_i} is composed of (k-1) occurrences of 1 and one occurrence of (d-k+1). As a result,

    ∑_d_1,…,d_N
    ∑ d_i = d∏_i = 1^N2d_id_i≤∑_k = 1^d Nkd-1k-1· 2^k-12(d-k+1)d-k+1∑_k = 1^d S(k).

Since

    S(k) = Nkd-1k-1· 2^k-12(d-k+1)d-k+1≤N^k/k!dk· 2^k2(d-k)d-k

and Stirling's formula gives

    k!    ≳(k/e)^k,
    dk   ≲d^d/k^k (d-k)^d-k, 
    2(d-k)d-k   ≲ 4^d-k, 
     d!    ≲√(d)(d/e)^d,

substituting into (<ref>) and denoting k = (1-η)d, we find

    S(k) ((2N)^d/d!)^-1   ≲(Ne/k)^k·d^d/k^k (d-k)^d-k· 2^2d-k·(1/√(d)(2Ne/d)^d)^-1
        = √(d)(2d/Neη)^η d1/(1-η)^2(1-η)d.

Now, for η∈ (0,1], denote 

    h(η) (2d/Neη)^η1/(1-η)^2(1-η).

Then,

    ∑_d_1,…,d_N
    ∑ d_i = d∏_i = 1^N2d_id_i≤∑_k = 1^d S(k) ≤ dmax_1≤ k≤ dS(k) ≲ d^3/2(2N)^d/d![sup_η∈ (0,1]h(η)]^d.

The last step is to evaluate h(η). Note that

    d/d η [log h(η)] = log(2d/N)-log(η)+2log(1-η).

Since 2d/N≤2D_n/N = o(1) as n→∞, for large N the unique maximizer of h has the form

    η^* = η^*(N,d) = (2-o(1))d/N,

and consequently

    sup_η∈(0,1]h(η) = h(η^*) = (2/e(2-o(1)))^η^*(1-η^*)^-2(1-η^*)≤ e^4η^* (1-η^*)≤ e^c_2d/N.

Substituting into (<ref>) then completes the proof.


 Suppose (<ref>) holds. Let μ = -logλ̂. In the setting of Lemma <ref>, note that

    lim_n→∞D_n/N   = 0, 
    lim inf_n→∞(-1/2logλ̂_n) = -1/2log(lim sup_n→∞λ̂_n)    > 0.

For n large enough so that (c_2+1/γ)D_n/N < -1/2logλ̂_n, applying Lemma <ref> in the expression of (<ref>) yields

    L_n,N,β,^≤ D_2^2    = _v^(1),v^(2)∼_n^ρ∑_d = 0^⌊ D/2 ⌋[∑_d_1,…,d_N
    ∑ d_i = d∏_i = 1^N2d_id_i][β^2 v^(1),v^(2)^2/4]^d
       ≲∑_d = 1^⌊ D/2 ⌋[d^3/2(2N)^d/d! e^c_2 d^2/N][(β^2/4)^d (nρ)^-2d· 2d e^μ d+d^2/nnd(2d)!/2^dρ^2d]
       ≤∑_d = 1^⌊ D/2 ⌋ 2d^5/2e^c_2 dD/N(2d)!/4^d(d!)^2[e^μ+D/n(N/n) β^2]^d
       ≲∑_d = 1^∞ d^2(e^(c_2+1/γ)D/Ne^μλ̂^2)^d 
       ≤∑_d = 1^∞ d^2(λ̂^1/2)^d 
       = O(1),

where the last equality is by the assumption that lim sup_n→∞λ̂_n < 1.


 If (<ref>) holds, then λ̂_n≤ 1/√(3) for sufficiently large n. In the setting of Lemma <ref>, suppose (c_2+1/γ)D/N < 0.001. Then, substituting the estimates in Lemma <ref> (taking μ = λ̂) into (<ref>) gives

    L_n,N,β,^≤ D_2^2    ≲∑_d = 11^⌊ D/2 ⌋[d^3/2(2N)^d/d! e^c_2d^2/N][(β^2/4)^d (nρ)^-2d·√(d)e^d^2/n(11e/30)^d/2λ̂^-2dnd(2d)!/2^dρ^2d]
       ≲∑_d = 11^⌊ D/2 ⌋ d^2e^(c_2+1/γ)dD/N(2d)!/4^d(d!)^2[β^2(N/n)λ̂^-2]^d (11e/30)^d/2
       ≲∑_d = 11^∞ d^3/2(e^0.001√(11e/30))^d 
       = O(1),

completing the proof.



In (<ref>), only counting the terms {a_i} with d occurrences of 1 and (n-d) occurrences of zero yields

    _v^(1),v^(2)∼_n^ρ v^(1),v^(2)^2d≥ n^-2dnd(2d)!/2^d.

From Lemma <ref> we have that when D_n = o(n), for sufficiently large n and for any 1≤ d≤ D_n,

    nd = n(n-1)…(n-d+1)/n^d·n^d/d!≥1/2e^-d^2/nn^d/d!≥1/2e^-dD_n/nn^d/d!.

Substituting Lemma <ref>, (<ref>) and (<ref>) into (<ref>) yields, for sufficiently large n,

    L_n,N,β,^≤ D_2^2    ≥∑_d = 1^D_n[(2N)^d/d!][β^2d/4^d n^-2dnd(2d)!/2^d]
       ≳∑_d = 1^D_n(2d)!/4^d(d!)^2(β^2/γ)^d e^-dD_n/n
       ≳∑_d = 1^D_n1/√(d)(λ̂_n^2 e^-D_n/n)^d
       ≥∑_d = 1^D_n1/√(d)
       = ω(1),

since D_n = ω(1), lim inf_n→∞λ̂_n > 1 and e^-D_n/n→ 1.



Suppose ω(1) ≤ D_n ≤ o(n). If there exists a series of positive integers w_n = o(√(D_n)) such that

    lim inf_n→∞ 2λ̂_n^2 (D_n/2neρ_n^2)^1-1/w_n(w_n/(2w_n)!)^1/w_n > 1,

then L_n,N,β,^≤ D_n_2^2 →∞ as n→∞.



 If (<ref>) holds, we can choose an ϵ > 0 such that, for sufficiently large n, 

    2λ̂_n^2 (D_n/2neρ_n^2)^1-1/w_n(w_n/(2w_n)!)^1/w_n > 1+ϵ.

Let n satisfy the above inequality. Pick μ∈ (0,1) such that

    μ^1-1/w_n(1+ϵ) > 1,

which implies

    2 (μ D_n/2neρ_n^2)^1-1/w_n(w_n/(2w_n)!)^1/w_n > γ/β^2.

In the sum (<ref>), we only consider those d ∈ (μ D_n/2,⌊ D_n/2 ⌋) that are multiples of w_n. By Lemma <ref>, Lemma <ref>, and (<ref>),

    _v^(1), v^(2)∼_n^ρ(∑_d_1,…,d_N
    ∑ d_i = d∏_i = 1^N2d_id_i)(β^2 v^(1),v^(2)^2/4)^d 
       ≳(2N)^d/d![(β^2/4)^d(nρ)^-2dnd(2d)!√(w_n)/2^dρ^2d·(γ/β^2)^d]
       ≳(2d)!/4^d(d!)^2(Nγ/n)^d 
       ≳1/√(d).

Therefore

    L_n,N,β,^≤ D_n_2^2    ≳∑_μ D_n/2 < d < ⌊ D_n/2 ⌋
     w_n | d1/√(d)
       ≳1/√(w_n)(√( D_n/2w_n)-√(μ D_n/2w_n))
       = 1-√(μ)/√(2)√(D_n)/w_n
       = ω(1),

completing the proof.



For sufficiently large n, in Lemma <ref> we choose

    w_n = ⌈log(1/λ̂_n) ⌉ = ⌈log(√(γ)/β) ⌉

which is o(√(D_n)). Reorganizing the terms, the condition (<ref>) is satisfied if for sufficiently large n,

    ρ_n < 0.991/√(2e)(w_n· 2^w_n/(2w_n)!)^1/(w_n-1)√(D_n/n)λ̂_n^w_n/(w_n-1).

Notice that

    1/√(2e)(w_n· 2^w_n/(2w_n)!)^1/(w_n-1)   = Θ(w_n^-2) 
       = Θ(log^-2 (1/λ̂_n)), 
    λ̂_n^w_n/(w_n-1)   = λ̂_n·λ̂_n^1/(⌈log(1/λ̂_n) ⌉-1)
       = Θ(λ̂_n).

Therefore, there exists an absolute constant C such that, if

    ρ_n < C√(D_n/n)λ̂_nlog^-2(1/λ̂_n),

then (<ref>) is satisfied and the divergence of L_n,N,β,^≤ D_n_2^2 follows from Lemma <ref>.




 §.§ The Wigner Model

In this section, we use the bounds on A_d (Lemmas <ref>, <ref> and <ref>) to prove the upper bound (Theorem <ref>) and the lower bound (Theorem <ref>) on the Wigner LDLR (<ref>).
  We only work with those n for which (<ref>) holds. Let μ = -logλ. Note that

    lim_n→∞D_n/N   = 0, 
    lim inf_n→∞(-1/2logλ̂_n) = -1/2log(lim sup_n→∞λ̂_n)    > 0.

For large enough n that D_n/n < -1/2logλ_n, applying Lemma <ref> in the expression of (<ref>) yields

    L_n,λ,^≤ D_n^2    = _v^(1),v^(2)∼_n^ρ[∑_d = 0^D_n1/d!(n/2λ^2)^d  v^(1),v^(2)^2d]
       ≲∑_d = 1^D_n1/d!(n/2λ^2)^d(nρ)^-2d· 2d e^μ d+d^2/nnd(2d)!/2^dρ^2d
       ≲∑_d = 1^D_nd(2d)!/4^d(d!)^2(e^μ+D_n/nλ^2)^d
       ≲∑_d = 1^D_n1/√(d)(λ^1/2)^d 
       = O(1),

where the last equation is by the assumption lim sup_n→∞λ_n < 1.


 By the assumption (<ref>), for sufficiently large n we have λ_n≤ 1/√(3). Now, Theorem <ref>(b) immediately follows from Lemma <ref> (taking μ = λ), since for large enough n that D_n/n < 0.001, we have

    L_n,λ,^≤ D_n^2 
       ≲∑_d = 11^D_n1/d!(n/2λ^2)^d(nρ)^-2d·√(d)e^d^2/n(11e/30)^d/2λ^-2dnd(2d)!/2^dρ^2d
       ≲∑_d = 11^∞√(d)(2d)!/4^d(d!)^2e^d^2/n(11e/30)^d/2
       ≲∑_d = 11^∞(e^D_n/n√(11e/30))^d 
       ≲∑_d = 11^∞(e^0.001√(11e/30))^d 
       =O(1),

completing the proof.



Substituting (<ref>) and (<ref>) into (<ref>) yields

    L_n,λ,^≤ D_n_2^2    ≥∑_d = 1^D_n1/d!(n/2λ^2)^d· n^-2dnd(2d)!/2^d
       ≳∑_d = 1^D_n(2d)!/4^d(d!)^2λ^2de^-dD_n/n
       ≳∑_d = 1^D_n1/√(d)(λ^2 e^-D_n/n)^d 
       ≥∑_d = 1^D_n1/√(d)
       = ω(1),

since D_n = ω(1), lim inf_n→∞λ_n > 1 and e^-D_n/n→ 1.



Suppose ω(1) ≤ D_n ≤ o(n). If there exists a series of positive integers w_n = o(√(D_n)) such that

    lim inf_n→∞ 2λ_n^2 (D_n/neρ_n^2)^1-1/w_n(w_n/(2w_n)!)^1/w_n > 1

then L_n,λ,^≤ D_n_2^2 →∞ as n→∞.



If (<ref>) holds, we can choose an ϵ > 0 such that for sufficiently large n,

    2λ_n^2 (D_n/neρ_n^2)^1-1/w_n(w_n/(2w_n)!)^1/w_n > 1+ϵ.

Let n satisfy the above inequality. Pick μ∈ (0,1) such that

    μ^1-1/w_n(1+ϵ) > 1.

In the sum (<ref>) we only consider those d > μ D_n that are multiples of w_n. For each of them, Lemma <ref> gives

    1/d!(n/2λ_n^2)^d 𝔼 v^(1),v^(2)^2d   ≳1/d!(n/2λ_n^2)^d· n^-2dnd(2d)!/2^d [2(d/neρ^2)^1-1/w_n(w_n/(2w_n)!)^1/w]^d
       ≥(2d)!/4^d(d!)^2·n(n-1)⋯ (n-d+1)/n^d·[2λ^2(μ D_n/neρ^2)^1-1/w_n(w_n/(2w_n)!)^1/w_n]^d
       ≳1/√(d).

Therefore,

    L_n,λ,^≤ D_n_2^2    ≳∑_μ D_n < d < D_n 
     w_n | d1/√(d)
       ≳1/√(w_n)(√(D_n/w_n)-√(μ D_n/w_n))
       = (1-√(μ))√(D_n)/w_n
       = ω(1),

completing the proof.


 For sufficiently large n, in Lemma <ref> we choose the positive integer 

    w_n = ⌈log(1/λ_n) ⌉,

which is o(√(D_n)). The divergence of L_n,λ,^≤ D_n_2^2 follows from the condition (<ref>), which is implied by the following sufficient condition: for sufficiently large n,

    ρ_n < 0.991/√(e)(w_n· 2^w_n/(2w_n)!)^1/(w-1)√(D_n/n)λ_n^w_n/(w_n-1).

Similar to the proof of Lemma <ref>, notice that

    1/√(e)(w_n· 2^w_n/(2w_n)!)^1/w-1 = Θ(w_n^-2) = Θ(log^-2 (1/λ_n));   λ_n^w_n/(w_n-1) = λ_n·λ_n^1/(⌈log(1/λ_n) ⌉-1) = Θ(λ_n),

Thus there exists an absolute constant C such that, if

    ρ_n < C√(D_n/n)λ_nlog^-2(1/λ_n),

then (<ref>) is satisfied and the divergence of L_n,λ,^≤ D_n_2^2 follows from Lemma <ref>.




§ ACKNOWLEDGMENTS

tocsectionAcknowledgments
We thank Samuel B. Hopkins, Philippe Rigollet, and Eliran Subag for helpful discussions.


tocsectionReferences
alpha






§ CHERNOFF BOUNDS

In this section, we present two Chernoff-type concentration inequalities used in our proofs.

Let u^(1),u^(2)∈ℝ^N be independent samples from (0,I_N). Then, for any 0 < t ≤ N/2,

    [| u^(1),u^(2)| ≥ t] ≤ 2exp( -t^2/4N).



Since by symmetry u^(1),u^(2) and - u^(1),u^(2) have the same distribution, it suffices to bound [ u^(1),u^(2)≥ t] for 0 < t≤ N/2. By Markov's inequality on the moment generating function, for any μ > 0,

    [ u^(1),u^(2)≥ t] ≤ e^μ u^(1),u^(2)/e^μ t = e^-μ t( e^μ x_1x_2)^N,
 
where x_1,x_2 are independent samples from (0,1). We compute

    e^μ x_1x_2 = 1/2π∬_^2exp(-x^2/2+μ xy -y^2/2)dx dy
    = (1-μ^2)^-1/2.

Take μ = t/N ∈ (0,1/2]. Note that 1-z ≥ e^-3z/2 on z∈ (0,1/4], and so 1-μ^2 ≥ e^-3μ^2/2. Hence

    [ u^(1),u^(2)≥ t] ≤exp(-t^2/N)·exp(3Nμ^2/4) = exp(-t^2/4N),

and the result follows.


The following result may be found in <cit.>.

For all 0 < z < 1,

    1/klog[χ_k^2 ≤ zk] ≤1/2(1-z+log z).

Similarly, for all z > 1,

    1/klog[χ_k^2 ≥ zk] ≤1/2(1-z+log z).




For all 0 < t ≤ 1/2,

    1/klog[|χ_k^2 -k| ≥ kt] ≤ -t^2/3.



It is easy to check that for t∈ (0,1/2],

    t+log(1-t)    ≤ -t^2/3
     -t+log(1+t)    ≤ -t^2/3.

Therefore, by Lemma <ref>,

    1/klog[|χ_k^2 -k| ≥ kt]    = 1/klog[χ_k^2 ≥ (1+t)k] +1/klog[χ_k^2 ≤ (1-t)k]
       ≤1/2(-t+log(1+t))+1/2(t+log(1-t)) 
       ≤  -t^2/3,

completing the proof.



