Detection of the number of principal components by extended AIC-type method
Jianwei Hu Jingfei Zhang Ji Zhu



Extended AIC-type method


,


Hu, Zhang and Zhu
Central China Normal Universityt1, University of Miamit2
 and University of Michigant3
Department of Statistics

Central China Normal University

Wuhan, 430079

China

e1


Department of Management Science

University of Miami

Coral Gables, FL 33124

USA

e2


Department of statistics

University of Michigan

Ann Arbor, MI 48109

USA

e3



Estimating the number of principal components is one of the fundamental problems in many scientific fields such as signal processing (or the spiked covariance model).
In this paper, we first demonstrate that, for fixed p, any penalty term of the form k'(p-k'/2+1/2)C_n may lead to an asymptotically consistent estimator under the condition that C_n→∞ and C_n/n→0. Compared with the condition in <cit.>, i.e., C_n/loglog n→∞ and C_n/n→0, this condition is significantly weakened. Then we propose to select the number of signals k by the iterated logarithm penalty (ILP).
We also extend our results to the case n,p→∞, with p/n→ c>0. In this case, for increasing k, we first investigate the limiting laws for the leading eigenvalues of the sample covariance matrix _n under the condition that λ_k>1+√(c), which extend the results in <cit.> and <cit.>. This includes the case λ_k→∞. At low SNR, since the AIC tends to underestimate the number of signals k,  the AIC should be re-defined in this case. As a natural extension of the AIC for fixed p, we propose the extended AIC (EAIC), i.e., the AIC-type method with tuning parameter γ=φ(c)=1/2+√(1/c)-log(1+√(c))/c, and demonstrate that the EAIC-type method, i.e., the AIC-type method with tuning parameter γ>φ(c), can select the number of signals k consistently. As a result, in the following two cases, (1) p fixed, n→∞,
(2) n,p→∞ with p/n→ 0, if the AIC is defined as the degeneration of the EAIC in the case n,p→∞ with p/n→ c>0, i.e., γ=lim_c→ 0+φ(c)=1, then we have essentially demonstrated that, the AIC tends to overestimate k. To achieve the consistency of the AIC-type method in the above two cases, γ>1 is required. On the other hand, in the case n,p→∞, with p/n→ c>0, we have φ(c)<1. Then in this case, we have actually explained why the AIC tends to underestimate k. Moreover, we show that the EAIC-type method is essentially tuning-free and outperforms the well-known KN estimator proposed in <cit.> and the BCF estimator proposed in <cit.>. Numerical studies indicate that the proposed method works well.


[class=MSC]
[Primary ]62E20
[; secondary ]60G05



Consistency
minimum description length
signal-to-noise ratio
spiked covariance model
Tracy-Widom distribution






§ INTRODUCTION



Detection of the number of principal components from a noisy data is one of the fundamental problems in many scientific fields. It is often the starting point for the signal parameter estimation problem such as signal processing <cit.>, wireless communications <cit.>, array processing <cit.> and finance <cit.> to list a few and has primary importance <cit.>. The most common approach to solving this problem is by using information theoretic criteria, and in particular minimum description length (MDL) <cit.>, Bayesian information criterion (BIC) and Akaike information criterion (AIC) <cit.>.

Consider the model

    (t)=𝐬(t)+𝐧(t)=(_1^⊤(t),_2^⊤(t))^⊤

where k<p, =((Ψ_1),⋯,(Ψ_k)), 𝐬(t)=(𝐬_1(t),⋯,𝐬_k(t))^⊤ and 𝐧(t)=(𝐧_1(t),⋯,𝐧_p(t))^⊤. In (<ref>), 𝐧(t) is the noise vector distributed independent of 𝐬(t) as multivariate normal with mean vector 0 and covariance matrix σ^2 _p. 𝐬(t) is distributed as multivariate normal with mean 0 and nonsingular matrix Ω and (Ψ_i): p× 1 is a vector of functions of the elements of unknown vector Ψ_i associated with i-th signal. Then, the covariance matrix Σ of (t) is given by

    Σ=Ω^⊤+σ^2 _p=([ _11 _12; _21 _22;     ]).

We assume that (t_1),⋯, (t_n) are independent observations on (t). Let λ_1≥⋯≥λ_k>λ_k+1=⋯=λ_p=λ=σ^2 denote the eigenvalues of . Note that (<ref>) is also called the spiked covariance model in <cit.>.

Using the well-known spectral representation theorem from linear algebra, we can express Σ as

    =∑_i=1^k(λ_i-λ)_i_i^⊤+λ_p,

where _i is the eigenvector of  with eigenvalue λ_i. Denoting by θ the parameter vector of the model, it follows that

    θ^⊤=(θ_1^⊤,θ_2^⊤),

where θ_1^⊤=(λ_1,⋯,λ_k,λ) and θ_2^⊤=(_1^⊤,⋯,_k^⊤).


With this  parameterization we now proceed to the derivation of the information theoretic criteria for the detection problem. Since the observation are regarded as statistically independent Gaussian random vectors with zero mean, their joint probability density is given by

    f(|θ)=∏_i=1^nf((t_i)|θ)=∏_i=1^n1/(2π)^p/2||^1/2exp-1/2^⊤(t_i)^-1(t_i),

where ^⊤=(^⊤(t_1),⋯,^⊤(t_n)).

Taking the logarithm, the log-likelihood function is given by

    log f(|θ)   =   ∑_i=1^nlog f((t_i)|θ)
       =   log∏_i=1^n1/(2π)^p/2||^1/2exp-1/2^⊤(t_i)^-1(t_i)
       =   -1/2nlog||-1/2n(^-1_n)-1/2nplog(2π).

where _n=1/n∑_i=1^n(t_i)(t_i)^⊤=1/n^⊤ is the sample covariance matrix.

Let ={λ_1,⋯,λ_p} and =^⊤, where  is orthogonal and its column _i is the eigenvector of   with eigenvalue λ_i. Also,  denote by d_1≥⋯≥ d_p the eigenvalues of _n. Let ={d_1,⋯,d_p} and

    _n=^⊤≜([ _11 _12; _21 _22;     ])=1/n([ _1^⊤_1 _1^⊤_2; _2^⊤_1 _2^⊤_2;        ]),

where  is orthogonal and its column _j is the eigenvector of _n  with eigenvalues d_j.


In the signal processing literature, the AIC estimates k by maximizing ℓ(k'):

    k̂=max_k'ℓ(k'),

where

    ℓ(k')=-1/2n(∑_i=1^k'log d_i+(p-k')logλ̂_k')-k'(p-k'/2+1/2).

Here λ̂_k'=1/(p-k')∑_i=k'+1^pd_i.


The BIC (or the MDL) estimates k by maximizing ℓ(k'):

    k̂=max_k'ℓ^B(k'),

where

    ℓ^B(k')=-1/2n(∑_i=1^k'log d_i+(p-k')logλ̂_k')-1/2k'(p-k'/2+1/2)log n.


Note that the AIC and the BIC are slightly different from that of <cit.>, where the number of unknown parameters was k'(p-k'/2) and was corrected by <cit.>.

The first consistency proof of the BIC was established in <cit.>. By using  Wilks' theorem,  <cit.> showed that any penalty term of the form k'(p-k'/2)C_n may lead to an asymptotically consistent estimator under the condition that C_n→∞ and C_n/n→0. Later, it was demonstrated that Wilks' theorem cannot be used in this problem in <cit.>. By using Taylor's expansion,   <cit.> gave an alternative proof and showed that, any penalty term of the form k'(p-k'/2+1/2)C_n may lead to an asymptotically consistent estimator under the condition that C_n/loglog n→∞ and C_n/n→0. Although  <cit.> proved that the BIC is consistent, the BIC fails to detect signals at low signal-to-noise ratio (SNR), hence underestimating the number of signals at small sample size. In contrast, while the AIC is able to detect low SNR signals, it has a non-negligible probability to overestimate the number of signals and thus is not consistent. To remedy the shortcoming of the AIC,  <cit.> proposed an modified AIC. The only difference between these two methods lies in that the penalty of the modified AIC is two times that of the AIC. In the presence of pure noise with no signals (i.e., k=0),  <cit.> showed that the modified AIC had a negligible overestimation probability for large n. However, for k>0,  <cit.> did not give an explanation why the modified AIC had a negligible overestimation probability. On the other hand, for small and medium n, the modified AIC tends to underestimate the number of signals k in our simulations.

 <cit.> and  <cit.> proposed a quite different method for estimating the number of signals, via a sequence of hypothesis tests, at each step testing the significance of the k-th eigenvalue as arising from a signal. The main tools used in these two papers are recent results from random matrix theory regarding both the distribution of noise eigenvalues and of signal eigenvalues in the presence of noise. In the absence of signals, the matrix n_n follows a Wishart distribution <cit.> with parameter n, p. In the joint limit n,p→∞, with p/n→ c>0, the distribution of the largest eigenvalue of _n converges to a Tracy-Widom distribution <cit.>,<cit.>,<cit.>,<cit.>,<cit.>. For fixed p, although    <cit.> proved the strong consistency of their estimator, our simulations show that when the SNR is low, for fixed p and medium n, the KN tends to underestimate the number of signals k.

Before further proceeding, we mention that a similar, if not identical problem, also appears in other literatures, likelihood ratio test statistic
<cit.>, Kac-Rice test <cit.>. For some other work in this topic, we refer to <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>  and so on.

In this paper, we first demonstrate that any penalty term of the form k'(p-k'/2+1/2)C_n may lead to an asymptotically consistent estimator under the condition that C_n→∞ and C_n/n→0. Compared with the condition in <cit.>, i.e., C_n/loglog n→∞ and C_n/n→0, this condition is significantly weakened. Then we propose to select the number of signals k by the iterated logarithm penalty (ILP). Moreover, we also extend our results to the case n,p→∞, with p/n→ c>0. In this case, for increasing k, we first investigate the limiting laws for the leading eigenvalues of the sample covariance matrix _n under the condition that λ_k>1+√(c), which extend the results in <cit.> and <cit.>. This includes the case λ_k→∞.
At low SNR, since the AIC tends to underestimate the number of signals k,  the AIC should be re-defined in this case. As a natural extension of the AIC for fixed p, we propose the extended AIC (EAIC), i.e., the AIC-type method with tuning parameter γ=φ(c)=1/2+√(1/c)-log(1+√(c))/c, and demonstrate that the EAIC-type method, i.e., the AIC-type method with tuning parameter γ>φ(c), can select the number of signals k consistently. Moreover, we show that the EAIC-type method is essentially tuning-free and outperforms the well-known KN estimator proposed in <cit.> and the BCF estimator proposed in <cit.>.


For the remainder of the paper, we proceed as follows.
In Section <ref>, we demonstrate that any penalty term of the form k'(p-k'/2+1/2)C_n may lead to an asymptotically consistent estimator under the condition that C_n→∞ and C_n/n→0.
We extend our results to the case n,p→∞, with p/n→ c>0 in Section <ref>. The numerical studies are given in Section <ref>.
Some further discussions are made in Section <ref>.
All proofs are given in Section <ref>.







§ FIXED P


In this section, we work under the classical setting where p, k and λ_i are all fixed as n →∞. We establish the consistency
of (<ref>) in the sense that it chooses the correct k with probability
tending to one when n goes to infinity.

We consider the following general criterion:

    k̂=max_k'ℓ(k'),

where

    ℓ(k')=-1/2n(∑_i=1^k'log d_i+(p-k')logλ̂_k')-k'(p-k'/2+1/2)C_n,

where C_n→∞ and C_n/n→0. Here λ̂_k'=1/(p-k')∑_i=k'+1^pd_i.

Throughout this paper, we assume that λ is known. Without loss of generality, we assume that λ=1. The signal-to-noise ratio is defined as

    SNR=λ_k-λ/λ=λ_k/λ-1=λ_k-1.


Define λ̂=1/(p-k)∑_i=k+1^pd_i. Note that for i=1,2,⋯, k, d_i and λ̂ are the maximum likelihood estimates of λ_i and λ, respectively. The following result which shows how close is the eigenvalue of Σ to its maximum likelihood estimate.

 For i=1,2,⋯, p,

    (λ_i-d_i)^2=O_p(p/n)=O_p(1/n),


    (p-k)(λ-λ̂)^2=O_p(p^2/n)=O_p(1/n).





Noting that for fixed p, it is difficult to demonstrate the consistency of (<ref>) directly (see the proof of Theorem <ref>). Similar to the discussion in <cit.> (see also <cit.>) , we consider an alternative method.





When k' is equal to k, ℓ(k') defined in (<ref>) can be re-written as:

    ℓ(k)=-1/2n∑_i=1^klog d_i-1/2n∑_i=k+1^p(d_i-1)-k(p-k/2+1/2)C_n.

where C_n→∞ and C_n/n→0.



Then we estimate k by maximizing ℓ̃(k'):

    k̂=max_k'ℓ̃(k'),

where 
    ℓ̃(k')=-1/2n∑_i=1^k'log d_i-1/2n∑_i=k'+1^p(d_i-1)-k'(p-k'/2+1/2)C_n.

where C_n→∞ and C_n/n→0.


Next, we establish the consistency of (<ref>).



Let ℓ̃(k') be the penalized likelihood function defined in (<ref>).

For k'<k,

    P(ℓ̃(k)>ℓ̃(k'))→1.

For k'>k,

    P(ℓ̃(k)>ℓ̃(k'))→1.





Especially, we propose to select the number of signals k by the iterated logarithm penalty (ILP). That is, we estimate k by maximizing ℓ(k'):

    k̂=max_k'ℓ(k'),

where

    ℓ(k')=-1/2n(∑_i=1^k'log d_i+(p-k')logλ̂_k')-γ k'(p-k'/2+1/2)loglog n,

where γ>0 is a tuning parameter.

As an alternative method, we may also estimate k by maximizing ℓ̃(k'):

    k̂=max_k'ℓ̃(k'),

where 
    ℓ̃(k')=-1/2n∑_i=1^k'log d_i-1/2n∑_i=k'+1^p(d_i-1)-γ k'(p-k'/2+1/2)loglog n,

where γ>0 is a tuning parameter.


Simulations show that (<ref>) and (<ref>) behave almost the same. For simplicity and comparisons with other methods, we consider (<ref>) in our simulations. For simplicity, we may use γ=1 in our simulations which gives good performance.

   <cit.> demonstrated that any penalty term of the form k'(p-k'/2+1/2)C_n leads to an asymptotically consistent estimator under the condition that  C_n/loglog n→∞ and C_n/n→0. Theorem <ref> shows that this condition can be significantly weakened as C_n→∞ and C_n/n→0. To show that our results are correct, we also consider different C_n's in our simulations in Section <ref>, such as ILP^1/2, which estimates k by maximizing ℓ̃(k'):

    k̂=max_k'ℓ̃(k'),

where 
    ℓ̃(k')=-1/2n∑_i=1^k'log d_i-1/2n∑_i=k'+1^p(d_i-1)-γ k'(p-k'/2+1/2)(loglog n)^1/2.















§ INFINITE P


In this section, we allow p, k and λ_k grow as n grows.



 §.§ Asymptotic behavior of the spiked sample eigenvalues


In this subsection, in the case n,p→∞, with p/n→ c>0, for increasing k, we investigate the limiting laws for the leading eigenvalues of the sample covariance matrix _n under the condition that λ_k>1+√(c), which extend the results in <cit.> and <cit.>. This includes the case λ_k→∞.

We wish to point out that the empirical spectral distribution (ESD) F(x) of  converges weakly to the limiting spectral distribution (LSD) H(x) as p→∞, that is,

    sup_x | F(x)-H(x)|=k/p,

where F(x)=1/p∑_i=1^pI_(-∞,x](λ_i), H(x)=I(1≤ x) and I_A(·) is the indicator function of A.

Note that the ESD of _n is given by

    F_n(x)=1/p∑_i=1^pI_(-∞,x](d_i).

With probability 1,

    F_n(x)d→F_c(x).

Here, for 0<c≤ 1, F_c(x) is given as

    F_c'(x)=f_c(x)={[ 1/2π xc√((b-x)(x-a)),         if  x∈ (a,b),;                    0,            otherwise, ].

where a=(1-√(c))^2 and b=(1+√(c))^2.

If c>1, F_c(x) has a point mass 1-1/c at the origin, that is,


    F_c(x)={[                   0,            if   x<0,;               1-1/c,         if   0≤ x<a,; 1-1/c+∫_a^xf_c(t)dt,         if   a≤ x<b, ].

where a and b are the same as in the case 0<c≤ 1. We remark that ∫_a^bf_c(t)dt=1 or 1/c according to 0<c≤ 1 or c>1 respectively.

By a result of  <cit.> (see also  <cit.>),

    sup_x | F_n(x)-F_c(x)|=O_p(k/n).

This implies that the ESD of _n also converges weakly to F_c(x) as n→∞ under the condition where the number of spikes, k, is fixed or k=o(n). From the MP law, we have the easy consequence that if i/p→α∈ (0,1), then d_ia.s.→μ_1-α, where μ_α is the α-th quantile of the MP law, that is F_c(μ_α)=α.


The i-th largest eigenvalue, λ_i, is said to be a distant spiked eigenvalue if ψ'(λ_i)>0 where ψ(λ_i)=λ_i+cλ_i/λ_i-1. Equivalently, λ_i is a distant spiked eigenvalue if λ_i>1+√(c).



Define x(t) as

    (t)=^1/2(t),

where (t) is distributed with mean 0 and covariance matrix _p. Thus, the covariance matrix of (t) is given by .

For fixed k, under the assumption that λ_1 is bounded,  according to <cit.>, the following result holds.


In the case n,p→∞ with p/n→ c>0, suppose that for any 1≤ i≤ n, 1≤ j≤ p, E(((t_i))_j^4)≤ C, λ_1 is bounded and λ=1.
If λ_i is a distant spiked eigenvalue, we have

    d_ia.s.⟶ψ(λ_i).



The above definition and result are a special case of a more general definition and result in <cit.> and <cit.>.

For fixed k, under the assumption that λ_k→∞ as p→∞, a new limiting result for the distant spiked eigenvalues was also established in <cit.>.



In the case n,p→∞ with p/n→ c>0, suppose that for any 1≤ i≤ n, 1≤ j≤ p, E(((t_i))_j^4)≤ C and λ=1. If λ_k→∞ as p→∞, for  any 1≤ i≤ k, we have

    d_i/ψ(λ_i)a.s.⟶1.





For fixed k, from Lemmas <ref> and <ref>, under the assumption that


  (1) at least one λ_i (1≤ i≤ k-1) satisfying λ_i→∞ as p→∞,

  (2)1+√(c)<λ_k<∞,

for any 1≤ i≤ k,

    d_i/ψ(λ_i)a.s.⟶1.

is still not established.




For increasing k, under the assumption that λ_k→∞ as p→∞,  <cit.> established the following result.


In the case n,p→∞ with p/n→ c>0, suppose that for any 1≤ i≤ n, 1≤ j≤ p, E(((t_i))_j^4)≤ C and λ_k+1≥⋯≥λ_p are bounded. If λ_k→∞ as p→∞ and k=min{o(n^1/6), o(λ_i^-1/2)}, for  any 1≤ i≤ k, we have

    d_i/λ_i-1=max{O_p(1/λ_i), O_p(k^4/n)}.



Based on Lemma <ref>,  <cit.> also established the central limit theorem for increasing k, under the assumption that λ_k→∞ as p→∞.



In the case n,p→∞ with p/n→ c>0, suppose that for any 1≤ i≤ n, 1≤ j≤ p, E(((t_i))_j^4)≤ C and λ_k+1≥⋯≥λ_p are bounded, there exists a positive constant C not depending on n such that λ_i-1/λ_i≥ C>1, i=1,2,…,k. If λ_k→∞ as p→∞ and k=min{o(n^1/6), o(λ_i^-1/2)}, for  any 1≤ i≤ k, we have

    √(n)(d_i-ψ(λ_i))/ψ(λ_i)d⟶ N(0, σ_i^2),

where σ_i^2 is defined in <cit.>.


For more work on this topic, we refer to <cit.>, <cit.>,  <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.> and <cit.>.

For increasing k, under the assumption that λ_k>1+√(c), we will extend the results in <cit.> and <cit.>.

Our first result gives the limits for the spiked eigenvalues of _n, d_1≥⋯≥ d_k, which will be used to establish the consistency of the following EAIC-type method.

Suppose that k=o(n^1/3) and p/n-c=O(√(k/n))(=o(n^-1/3)). If λ_k>1+√(c), for any 1≤ i≤ k,

    d_i/ψ(λ_i)a.s.⟶1,


where ψ(λ_i)=λ_i+cλ_i/λ_i-1.





To demonstrate this result, we need the following result which extends Lemma <ref>. This result can be of independent interest.


For any positive semi-definite matrix _n, denoting by l_1≥⋯≥ l_n the eigenvalues of _n, then there exist some universal constant C>0 and some constant ρ such that the following inequality holds

    P(||1/n_11^-1/2_1^⊤_n_1_11^-1/2-1/n(_n)||>t)
       ≤   2exp(Ck-nt^2/ρ∑_s=1^nl_s^2/n)

for all 0<t<ρ∑_s=1^nl_s^2/n/64l_1, where ||·|| is the spectral norm.




For any positive semi-definite matrix _n, denoting by l_1≥⋯≥ l_n the eigenvalues of _n, then there exist some universal constant C>0 and some constant ρ such that the following inequality holds



    P(||1/n_11^-1/2_1^⊤(+_n)_1Σ_11^-1/2-1/n(+_n)||>t)
       ≤   2exp(Ck-nt^2/ρ∑_s=1^n(1+l_s)^2/n),

for all 0<t<ρ∑_s=1^n(1+l_s)^2/n/64(1+l_1).


Our second result shows that the spiked eigenvalues of _n, d_1≥⋯≥ d_k are √(n/k)-consistent, which will be used to establish Theorem <ref>.

Compared with Theorem <ref>, we need one more condition: there exists a positive constant C not depending on n such that λ_i-1/λ_i≥ C>1, i=1,2,…,k.
This condition implies that the spiked eigenvalues are well-separated and is also used in <cit.>.



Suppose that k=o(n^1/4), p/n-c=O(√(k/n))(=o(n^-3/8)), and there exists a positive constant C not depending on n such that λ_i-1/λ_i≥ C>1, i=1,2,…,k. If λ_k>1+√(c), for any 1≤ i≤ k,

    d_i-ψ(λ_i)/λ_i=O_p(√(k/n)).



To demonstrate this result, we need the following Lemma <ref>. We give some notations.

Define

    m_1(d)=∫x/d-xdF_c(x),


    m_2(d)=∫x^2/(d-x)^2dF_c(x),


    m_3(d)=∫x/(d-x)^2dF_c(x),


    _̋n(d)=1/n_2(d -_22)^-1_2^⊤,


    K_n(d)=1/n_1^⊤(+_̋n(d))_1,


    (d)=1/n_11^-1/2_1^⊤(+_̋n(d))_1_11^-1/2-(1+cm_1(d)).

Let 
    [ _11=(λ_1,…,λ_k)^⊤, ]

where  is an orthogonal matrix.

Obviously, 
    _11^-1=(1/λ_1,…,1/λ_k)^⊤.


Also, let

    (d_i)   =   (a_st)_k× k≜^⊤_11^-1/2(d_i- K_n(d_i))_11^-1/2
       =   ^⊤(_11^-1d_i-1/n_11^-1/2_1^⊤(+_̋n(d_i))_1_11^-1/2)
       =   (d_i-(1+cm_1(d_i))λ_1/λ_1,…,d_i-(1+cm_1(d_i))λ_k/λ_k)-^⊤(d_i).


According to <cit.>, for λ_i>1+√(c),  we have m_1(ψ(λ_i))=1/λ_i-1.
Thus, for the j-th diagonal element  with index j≠ i,

    e_j≜d_i-(1+cm_1(d_i))λ_j/λ_j   =   d_i/λ_j-(1+cm_1(d_i))
       →   ψ(λ_i)/λ_j-(1+cm_1(ψ(λ_i)))
       =   (λ_i/λ_j-1)(1+c/λ_i-1)

is uniformly bounded away from 0.



Let

    ((d_i))=∏_j=1^kê_j≜∏_j=1^k(e_j+ϵ_j),

where ê_j's are eigenvalues of matrix (d_i).

Similar to the proof of (<ref>) in Theorem <ref>, we have

    ∑_i=1^k|ϵ_j|=O_p(√(k^3/n)).



For j≠ i, under the condition that k=o(n^1/3), ê_j is uniformly bounded away from 0 since the contributions from ^⊤(d_i) tend to O_p(√(k^3/n))=o_p(1).

By (<ref>), for sufficiently large n, ((d_i))=0 if and only if ê_i=0.

However, ê_i=0 does not imply that a_ii=0.


For any 1≤ i≤ k, the following lemma gives the order of the i-th diagonal element of matrix (d_i),  which is the key to establish Theorems <ref> and <ref>.




Suppose that k=o(n^1/4), p/n-c=O(√(k/n))(=o(n^-3/8)), and there exists a positive constant C not depending on n such that λ_i-1/λ_i≥ C>1, i=1,2,…,k. If λ_k>1+√(c), for any 1≤ i≤ k,

    a_ii=d_i-(1+cm_1(d_i))λ_i/λ_i-[^⊤(d_i)]_ii=O_p(k^2/n).





Our third result demonstrates that the spiked eigenvalues of _n, d_1≥⋯≥ d_k are asymptotical normal.

Suppose that k=o(n^1/4), p/n-c=o(n^-1/2), and there exists a positive constant C not depending on n such that λ_i-1/λ_i≥ C>1, i=1,2,…,k. If λ_k>1+√(c), for any 1≤ i≤ k,

    √(n)(d_i-ψ(λ_i))/λ_ia.s.⟶T(ψ(λ_i)),

where

    T(ψ(λ_i))≜(1+cm_3(ψ(λ_i))λ_i)^-1[^⊤_11^-1/2 R_n(ψ(λ_i))_11^-1/2]_ii,

and

    R_n(ψ(λ_i))=1/√(n)[_1^⊤(+_̋n(ψ(λ_i)))_1-(+_̋n(ψ(λ_i)))_11].


For any 1≤ i≤ k, the limiting distribution of T(ψ(λ_i)) is normal. For more details, we refer to <cit.> and <cit.>.

Especially, we consider the case _11 is diagonal as discussed in <cit.>. Similar to the discussion in <cit.>, for any 1≤ i≤ k, we have

    √(n)(d_i-ψ(λ_i))/λ_id⟶N(0,σ_λ_i^2),

where σ_λ_i^2=2-2c/(λ_i-1)^2.


Note that throughout this paper we assume that (t) is normal. Under the condition that for any 1≤ i≤ n, 1≤ j≤ k, E((_11^-1/2_1(t_i))_j^4)≤ C, we have the following results.



For any positive semi-definite matrix _n, denoting by l_1≥⋯≥ l_n the eigenvalues of _n, suppose that for any 1≤ i≤ n, 1≤ j≤ k, E((_11^-1/2_1(t_i))_j^4)≤ C hold,  we have

    ||1/n_11^-1/2_1^⊤_n_1_11^-1/2-1/n(_n)I||=√(1/n∑_s=1^nl_s^2)O_p(k/√(n)).






For any positive semi-definite matrix _n, denoting by l_1≥⋯≥ l_n the eigenvalues of _n, suppose that for any 1≤ i≤ n, 1≤ j≤ k, E((_11^-1/2_1(t_i))_j^4)≤ C hold,  we have

    ||1/n_11^-1/2_1^⊤(+_n)_1_11^-1/2-1/n(+_n)I||=√(1/n∑_s=1^n(1+l_s)^2)O_p(k/√(n)).





Suppose that k=o(n^1/4), p/n-c=O(k/√(n))(=o(n^-1/4)), and for any 1≤ i≤ n, 1≤ j≤ k, E((_11^-1/2_1(t_i))_j^4)≤ C. If λ_k>1+√(c), for any 1≤ i≤ k,

    d_i/ψ(λ_i)a.s.⟶1,


where ψ(λ_i)=λ_i+cλ_i/λ_i-1.


Suppose that k=o(n^1/5), p/n-c=O(k/√(n))(=o(n^-3/10)), for any 1≤ i≤ n, 1≤ j≤ k, E((_11^-1/2_1(t_i))_j^4)≤ C, and there exists a positive constant C not depending on n such that λ_i-1/λ_i≥ C>1, i=1,2,…,k. If λ_k>1+√(c), for any 1≤ i≤ k,

    d_i-ψ(λ_i)/λ_i=O_p(k/√(n)).





Suppose that k=o(n^1/5), p/n-c=O(k/√(n))(=o(n^-3/10)), for any 1≤ i≤ n, 1≤ j≤ k, E((_11^-1/2_1(t_i))_j^4)≤ C, and there exists a positive constant C not depending on n such that λ_i-1/λ_i≥ C>1, i=1,2,…,k. If λ_k>1+√(c), for any 1≤ i≤ k,

    a_ii=d_i-(1+cm_1(d_i))λ_i/λ_i-[^⊤(d_i)]_ii=O_p(k^3/n).





Suppose that k=o(n^1/6), p/n-c=o(n^-1/2), for any 1≤ i≤ n, 1≤ j≤ k, E((_11^-1/2_1(t_i))_j^4)≤ C, and there exists a positive constant C not depending on n such that λ_i-1/λ_i≥ C>1, i=1,2,…,k. If λ_k>1+√(c), for any 1≤ i≤ k,

    √(n)(d_i-ψ(λ_i))/λ_ia.s.⟶(1+cm_3(ψ(λ_i))λ_i)^-1[^⊤_11^-1/2 R_n(ψ(λ_i))Σ_11^-1/2]_ii,

where

    R_n(ψ(λ_i))=1/√(n)[_1^⊤(+_̋n(ψ(λ_i))X_1-(+_̋n(ψ(λ_i))_11].






 §.§ Extended AIC

In this subsection, we extend our results to the case n,p→∞, with p/n→ c>0. In this case, the AIC tends to underestimate the number of signals k (see Section <ref>). It is well-known that, for fixed p, the AIC tends to overestimate the number of signals k.  This implies that, in the case n,p→∞, with p/n→ c>0, the AIC does not work and should be re-defined. In view of this, we propose an extended AIC (EAIC).

The EAIC estimates k by maximizing ℓ(k'):

    k̂=max_k'ℓ(k'),

where

    ℓ(k')=-1/2n(∑_i=1^k'log d_i+(p-k')logλ̂_k')-γ k'(p-k'/2+1/2),

where γ=φ(c)=1/2+√(1/c)-log(1+√(c))/c.

 Note that when γ is equal to 1, (<ref>) is the AIC. In the following, we first show that the EAIC is a natural extension of the AIC for fixed p. In the next subsection, we further explain that, for the EAIC, why γ should be defined as φ(c)=1/2+√(1/c)-log(1+√(c))/c.

In the following two cases,


  (1)p fixed, n→∞,

  (2) n,p→∞ with p/n→ 0,
 c is defined as c=p/n. Noting that c→ 0+, by Taylor's expansion, we get

    [                            φ(c)                               =        1/2+√(1/c)-log(1+√(c))/c;                                                               = 1/2+√(1/c)-(√(c)-c/2(1+o(1)))/c;                                                               =                          1+o(1);                                                               →                              1. ]


As a result, in the following three cases,


  (1) p fixed, n→∞,

  (2) n,p→∞ with p/n→ 0,

  (3) n,p→∞ with p/n→ c, where 0<c<∞,

approximately, the AIC may be uniformly written as (<ref>).

Note that the EAIC is different from that of <cit.> which proposes the following BCF for estimating the number of principal components k.

In the case 0<c<1,  <cit.> estimates k by minimizing ℓ^1(k'):

    k̂=min_k'ℓ^1(k'),

where ℓ^1(k')=(p-k')logd̅_k'-∑_i=k'+1^plog d_i-(p-k'-1)(p-k'+2)/n, d̅_k'=∑_i=k'+1^pd_i/(p-k').

Note that,  the BCF defined in (<ref>) is essentially equivalent to the AIC and is referred as the AIC in <cit.>. In the next subsection, we will show that,  at low SNR, the BCF defined in (<ref>) tends to underestimate the number of signals k, especially in the case c→1-  (see also Section <ref>).

In the case c>1,  <cit.> replaces the p in (<ref>) by n-1
 and estimates k by minimizing ℓ^2(k'):

    k̂=min_k'ℓ^2(k'),

where ℓ^2(k')=(n-1-k')logd̅_k'-∑_i=k'+1^n-1log d_i-(n-k'-2)(n-k'+1)/p, d̅_k'=∑_i=k'+1^n-1d_i/(n-1-k').

Simulations show that, the BCF defined in (<ref>) is better than the AIC and is referred as the quasi-AIC in <cit.>. However, at low SNR, the BCF defined in (<ref>) still tends to underestimate the number of signals k, especially in the case c→1+  (see also Section <ref>).




 §.§ Consistency of the EAIC-type method


Noting that the EAIC tends to overestimate the number of signals k, we propose the EAIC-type method which chooses the correct k with probability
tending to one when n goes to infinity.

The EAIC-type method estimates k by maximizing ℓ(k'):

    k̂=max_k'ℓ(k'),

where

    ℓ(k')=-1/2n(∑_i=1^k'log d_i+(p-k')logλ̂_k')-γ k'(p-k'/2+1/2),

where γ>φ(c)=1/2+√(1/c)-log(1+√(c))/c.





To achieve the consistency of (<ref>), we  need two additional conditions,

    ψ(λ_k)-1-logψ(λ_k)>2γ c,


    λ_k>1+√(c),

where ψ(λ_k)=λ_k+cλ_k/λ_k-1.

Since there is a tuning parameter γ in (<ref>), λ_k>1+√(c) does not implies ψ(λ_k)-1-logψ(λ_k)>2γ c. However, by γ>φ(c) and ψ(λ_k)-1-logψ(λ_k)>2γ c, we have λ_k>1+√(c).

Next, we establish the consistency of (<ref>).


  In the case n,p→∞ with p/n→ c>0, let ℓ(k') be the penalized likelihood function defined in (<ref>). Suppose that k=o(n^1/3) and p/n-c=O(√(k/n))(=o(n^-1/3)), and the number of candidate models, q, satisfies q = o(p).

For k'<k, if (<ref>) and (<ref>) hold (or if γ>φ(c) and (<ref>) hold),

    P(ℓ(k)>ℓ(k'))→1.

For k'>k, if γ>φ(c),

    P(ℓ(k)>ℓ(k'))→1.

where φ(c)=1/2+√(1/c)-log(1+√(c))/c.


Note that throughout this paper, we assume that (t) is normal. Under the more general condition that for any 1≤ i≤ n, 1≤ j≤ k, E((_11^-1/2_1(t_i))_j^4)≤ C, we have the following result.

  In the case n,p→∞ with p/n→ c>0, let ℓ(k') be the penalized likelihood function defined in (<ref>). Suppose that k=o(n^1/4), p/n-c=O(k/√(n))(=o(n^-1/4)), for any 1≤ i≤ n, 1≤ j≤ k, E((_11^-1/2_1(t_i))_j^4)≤ C, and the number of candidate models, q, satisfies q = o(p).

For k'<k, if (<ref>) and (<ref>) hold (or if γ>φ(c) and (<ref>) hold),

    P(ℓ(k)>ℓ(k'))→1.

For k'>k, if γ>φ(c),

    P(ℓ(k)>ℓ(k'))→1.

where φ(c)=1/2+√(1/c)-log(1+√(c))/c.


Note that for fixed p and n→∞, γ>φ(c) becomes γ>1. It is well-known that, for fixed p, the AIC tends to overestimate the number of signals k. To achieve the consistency of the AIC-type method, the tuning parameter γ should be larger than one. That is, in the case n,p→∞ with p/n→ c>0, the EAIC can be seen as a natural extension of the AIC for fixed p. There is a similar discussion in the case n,p→∞ with p/n→ 0.


 In the following two cases,


  (1) p fixed, n→∞,

  (2) n,p→∞ with p/n→ 0,
 if the AIC is defined as the degeneration of the EAIC in the case n,p→∞ with p/n→ c>0, i.e., γ=lim_c→ 0+φ(c)=1, then we have essentially demonstrated that, the AIC tends to overestimate k. To achieve the consistency of the AIC-type method in the above two cases, γ>1 is required.




By Theorem <ref>, in the case n,p→∞ with p/n→ c>0, we have actually explained why the AIC tends to underestimate the number of signals k.
We note that γ>φ(c) is required in our theoretic analysis, to avoid overestimating the number of signals k. On the other hand, (<ref>) requires that γ cannot be too large, to avoid underestimating the number of signals k. Noting that, in the case n,p→∞ with p/n→ c>0, we always have φ(c)<1. Thus, the penalty of the AIC may be too large. As a result, it is not surprising that the AIC may tend to underestimate the number of signals k.



Since γ>φ(c) is required in our theoretic analysis, we set γ=1.1φ(c) for simulation studies which gives good performance. In this sense, the EAIC-type method is essentially tuning-free.

For fixed k, to prove the consistency of the BCF,  <cit.> gives the following two conditions,


  (1) in the case 0<c<1,

    ψ(λ_k)-1-logψ(λ_k)>2c    and


  (2) in the case c>1,

    ψ(λ_k)/c-1-log(ψ(λ_k)/c)>2/c.




Compared with (<ref>) and (<ref>), (<ref>) contains a tuning parameter γ. By simulations, we found that when c≥0.12, we have γ=1.1φ(c)<1. That is, when 0.12≤ c≤1, (<ref>) is weaker than (<ref>). This implies that, when 0.12≤ c≤1, in the case where (<ref>) holds while (<ref>) does not hold, the EAIC-type method is better than the BCF. When c>1, although it is difficult to compare (<ref>) with (<ref>), by simulations, we found that there do exist cases where (<ref>) holds while (<ref>) does not hold (see Section <ref>). This implies that, in the case  where (<ref>) holds while (<ref>) does not hold, the EAIC-type method is still better than the BCF. We also note that in the case c>1, if c is close to one, (<ref>) and (<ref>) are almost the same. Thus, in the case where c is close to one, the EAIC-type method outperforms the BCF. The reason lies in that, when c=p/n is close to one, (<ref>) is weaker than (<ref>) and (<ref>). As a result, the EAIC-type method outperforms the BCF at least in the following two cases,


  (1) in the case 0.12≤ c≤1,

  (2) in the case c>1, and c is close to one.

Moreover, note that (<ref>) is essentially equivalent to the AIC. Since when c≥0.12, we have γ=1.1φ(c)<1. That is, when (1) 0.12≤ c≤ 1; (2) c is larger than one and is close to one, the BCF tends to underestimate the number of signals k.

On the other hand, in the case where c is larger than one, simulations show that the EAIC-type method is still comparable to the BCF (see Section <ref>).



For fixed k,  <cit.>  established the consistency of the BCF. Along the way of the proposed method in this paper, it is possible to extend their results to increasing k.




Although the BCF  defined in (<ref>) and (<ref>) is tuning-free and simulation results are encouraging,  <cit.> essentially defines two different criteria in the case 0<c<1 and c>1, respectively. As a result, to achieve the consistency of the BCF, two different consistency conditions are required (i.e., (<ref>) and (<ref>)). Compared with (<ref>) and (<ref>), (<ref>) with γ>φ(c) is a natural extension of the AIC for fixed p. We also note that, both the formula (<ref>) and the consistency condition (<ref>) are more simple.



































§ EXPERIMENTS



We compare the performance of the ILP with the BIC, the AIC and the KN in a series of simulations with λ=σ^2=1. As suggested in <cit.>, the confidence level α=10^-4 was used in the KN. Our performance measure is the probability of the successful recovery of the number of signals,

    Pr(k̂=k).

We restrict our attention to candidate values for the true number of signals in the range k'∈{0,1,…,min{p-1,15}} in simulations. Each simulation in this section is repeated 200 times.

In our experiments, we define x(t) as

    (t)=^1/2(t),

where ={λ_1,⋯, λ_k,λ,⋯,λ} and (t) is distributed as multivariate normal with mean 0 and covariance matrix _p. Thus, the covariance matrix of (t) is given by .




 §.§ Fixed p


Simulation 1. In this simulaiton, we investigate how the accuracy of the ILP changes as the SNR varies.
We set p=12, k=3, γ=1, λ_1=⋯=λ_k-1=1+2SNR, λ_k=1+SNR and λ=1. We let SNR=δ√(4(p-k/2+1/2)loglog n/n), where δ increases from 1 to 2.

We compare the ILP with the BIC, the AIC, the modified AIC and the KN.  From Table <ref>, we can see that, when δ≥1.5, the ILP outperforms other methods in general.




 §.§ Infinite p

Note that the AIC-type method with γ>φ(c) is called the EAIC-type method in this paper. In simulations, we set γ=1.1φ(c).

Simulation 2. In this simulation, in the case n,p→∞ with p<n, we investigate how the accuracy of the EAIC-type method changes as the SNR varies.
We set n=500, p=200, k=10, λ_1=⋯=λ_k-1=1+2SNR, λ_k=1+SNR and λ=1. Note that when c=p/n=0.4, γ=1.1φ(c)=0.94 and √(c)=0.632. We let δ increase from 0.5 to 2.5. We also note that when δ=0.5, (<ref>) does not hold, while when 1≤δ≤2.5, (<ref>) holds. It can be seen from Table <ref> that, when
δ=0.5, the success rate of the EAIC-type method is low. On the other hand, when δ≥ 1, the success rate is high. This simulation verifies that Theorem <ref> is correct, i.e., if (<ref>) and (<ref>) hold, the EAIC-type method is consistent. It can also be seen from Table <ref> that, in the case n,p→∞ with p<n, the EAIC-type method is  better than the KN and the BCF , especially for δ≤1.

Note that when δ=1, (<ref>) holds while (<ref>) does not hold. Thus, it not surprising that in this case, the EAIC-type method is better than the BCF. On the other hand, although (<ref>) holds, ψ(λ_k)-1-logψ(λ_k)-2γ c=0.017 is close to zero. As a result, when
δ=1, the success rate of the EAIC-type method is not high.

Simulation 3. In this simulation, in the case n,p→∞ with p>n, we investigate how the accuracy of the EAIC-type method changes as the SNR varies.
We set n=200, p=500, k=10, λ_1=⋯=λ_k-1=1+2SNR, λ_k=1+SNR and λ=1. Note that when c=p/n=2.5, γ=1.1φ(c)=0.83 and √(c)=1.581. We let δ increase from 1.5 to 4.5. We also note that when δ=1.5, 2.5, (<ref>) does not hold, while when 2.68≤δ≤4.5, (<ref>) holds. It can be seen from Table <ref> that, when
δ=1.5, 2.5, the success rate of the EAIC-type method is low. On the other hand, when δ≥ 3.5, the success rate is high. This simulation verifies that Theorem <ref> is correct, i.e., if (<ref>) and (<ref>) hold, the EAIC-type method is consistent. It can also be seen from Table <ref> that, in the case n,p→∞ with p>n, the EAIC-type method is  better than the KN and the BCF, especially for δ≤3.5.

Note that when δ=2.68, (<ref>) holds while (<ref>) does not hold. Thus, it not surprising that in this case, the EAIC-type method is better than the BCF. On the other hand, although (<ref>) holds, ψ(λ_k)-1-logψ(λ_k)-2γ c=0.0085 is close to zero. As a result, when
δ=2.68, the success rate of the EAIC-type method is not high.


Simulation 4. In this simulation, in the case n,p→∞ with p=n, we investigate how the accuracy of the EAIC-type method changes as the SNR varies.
We set n=200, p=200, k=10, λ_1=⋯=λ_k-1=1+2SNR, λ_k=1+SNR and λ=1. Note that when c=p/n=1, γ=1.1φ(c)=0.89 and √(c)=1. We let δ increase from 1 to 3. We also note that when δ=1,1.5, (<ref>) does not hold, while when 2≤δ≤3, (<ref>) holds. It can be seen from Table <ref> that, when
δ=1,1.5, the success rate of the EAIC-type method is low. On the other hand, when δ≥ 2, the success rate is high. This simulation verifies that Theorem <ref> is correct, i.e., if (<ref>) and (<ref>) hold, the EAIC-type method is consistent. It can also be seen from Table <ref> that, in the case n,p→∞ with p=n, the EAIC-type method is better than the KN and the BCF, especially for δ≤2.5.

Note that when δ=2, (<ref>) holds while (<ref>) does not hold. Thus, it not surprising that in this case, the EAIC-type method significantly outperforms  the BCF.

Simulation 5. In this simulation, in the case n,p→∞, at low SNR, i.e., SNR=2√(p/n), we further compare the EAIC-type method with the BCF.
We set k=10, λ_1=⋯=λ_k-1=1+2SNR, λ_k=1+SNR and λ=1. It can be seen from Tables <ref> and  <ref> that, when
c=p/n is less than or close to one, the EAIC-type method outperforms the BCF. On the other hand, when c=p/n is larger than one, the EAIC-type method is still comparable to the BCF.




§ DISCUSSION


In this paper, we have demonstrated that any penalty term of the form k'(p-k'/2+1/2)C_n may lead to an asymptotically consistent estimator under the condition that C_n→∞ and C_n/n→0. Compared with the condition in <cit.>, i.e., C_n/loglog n→∞ and C_n/n→0, this condition has been significantly weakened.  We have also extended our results to the case n,p→∞, with p/n→ c>0. In this case, for increasing k, we have investigated the limiting law for the leading eigenvalues of the sample covariance matrix _n under the condition that λ_k>1+√(c). This includes the case λ_k→∞. At low SNR, since the AIC tends to underestimate the number of signals k,  the AIC should be re-defined in this case. As a natural extension of the AIC for fixed p, we have proposed the extended AIC (EAIC), i.e., the AIC-type method with tuning parameter γ=φ(c)=1/2+√(1/c)-log(1+√(c))/c, and demonstrated that the EAIC-type method, i.e., the AIC-type method with tuning parameter γ>φ(c), can select the number of signals k consistently. In simulations, we set γ=1.1φ(c) which gives good performance and outperforms the KN proposed in <cit.> and the BCF proposed in <cit.>. That is, the EAIC-type method is essentially tuning-free.


We have noted that in the case n,p→∞ with p/n→ c>0, if (<ref>) or (<ref>) does not hold, γ=1.1φ(c) tends to underestimate the number of of principal components k. As a result, 0<γ<1.1φ(c) may be a better choice. In this case, we may use the cross-validation method to choose the tuning parameter γ, which will be explored in  future work.




§ APPENDIX

This section presents auxiliary lemmas used in this paper and proofs for main results stated in this paper.


 §.§ Auxiliary Lemmas






Let ,∈_n, assume that ,+ are both normal, let λ_1,…,λ_n be the eigenvalues of  in some given order, and let λ̂_1,…,λ̂_n be the eigenvalues of + in some given order. There is a permutation σ(·) of the integers 1,…,n such that

    ∑_l=1^k(λ̂_σ(i)-λ_i)^2≤ ||||_2^2=(^*).






Suppose =(_1,⋯, _p)' is sub-Gaussian with constant ρ>0 and with mean 0 and covariance matrix . Let ^(1),⋯, ^(n) be n independent copies of Y. Then there exist some universal constant C>0 and some constant ρ_1 depending ρ, such that the sample covariance matrix of ^(1),⋯, ^(n), , satisfies

    P(->t)≤ 2exp(-nt^2ρ_1+Cp),

for all 0<t<ρ_1. Here · is the spectral norm.

For more work on this topic, the reader is referred to <cit.>, <cit.>, <cit.> and <cit.>.



Let  and  be p× p Hermitian matrices with eigenvalues ordered as λ_1≥⋯≥λ_p and μ_1≥⋯≥μ_p, respectively. Then

    sup_1≤ j≤ p|μ_j-λ_j|≤-,

where · is the spectral norm.






 §.§ Proofs for Main Results




By Lemmas <ref> and <ref>, for i=1,2,⋯, p,

    |λ_i-d_i|≤max_1≤ i≤ p|ϕ_i|=||_n-||=O_p(√(p/n)).

Thus, for i=1,2,⋯, p, we have

    (λ_i-d_i)^2=O_p(p/n)=O_p(1/n).

For fixed p, by Lemmas <ref> and <ref>,

    (p-k)(λ-λ̂)^2   =   (p-k)(λ-1/p-k∑_i=k+1^pd_i)^2
       =   (∑_i=k+1^p(λ-d_i))^2/(p-k)
       ≤   pmax_1≤ i≤ p(d_i-λ_i)^2
       ≤   pmax_1≤ i≤ p|ϕ_i|^2
       =   p||_n-||^2
       =   O_p(p^2/n)








Noting that (p-k)(λ-λ̂)^2=O_p(1/n), by Taylor's expansion, we get

    log L_k   =   -1/2n(∑_i=1^klog d_i+(p-k)logλ̂)
       =   -1/2n∑_i=1^klog d_i-1/2n(p-k)log((1/p-k∑_i=k+1^pd_i-λ)+λ)
       =   -1/2n∑_i=1^klog d_i-1/2n(p-k)log((1/p-k∑_i=k+1^pd_i-1)+1)
       =   -1/2n∑_i=1^klog d_i-1/2n(p-k)((1/p-k∑_i=k+1^pd_i-1)
          -1/2(1/p-k∑_i=k+1^pd_i-1)^2(1+o(1)))
       =   -1/2n∑_i=1^klog d_i-1/2n∑_i=k+1^pd_i+1/2n(p-k)+O_p(1)
       =   -1/2n∑_i=1^klog d_i-1/2n∑_i=k+1^p(d_i-1)+O_p(1).

Noting that C_n→∞, (<ref>) holds.




Let logL̃_k=-1/2n∑_i=1^klog d_i-1/2n∑_i=k+1^p(d_i-1).


Suppose k'<k. We have

    ℓ̃(k)-ℓ̃(k')=logL̃_k-logL̃_k'-γ(k-k')(p-(k+k')/2+1/2)C_n.





By Lemmas <ref> and <ref>, for i=1,2,⋯, p,

    |λ_i-d_i|≤max_1≤ i≤ p|ϕ_i|=||_n-||=O_p(√(1/n)).


Since λ_k>1, for large n,

    λ_k-1-logλ_k>ϵ,

where ϵ>0 is a sufficiently small constant.

Since

    (d_i-1-log d_i)-(λ_i-1-logλ_i)P⟶0,

for large n, we have

    P{ℓ̃(k)>ℓ̃(k')}   =   P{logL̃_k-logL̃_k'>γ(k-k')(p-(k+k')/2+1/2)C_n}
       =   P{1/2n∑_i=k'+1^k(d_i-1-log d_i)>γ(k-k')(p-(k+k')/2+1/2)C_n}
       ≥   P{1/2n(k-k')(λ_k-1-logλ_k)+o_p(n)
          >γ(k-k')(p-(k+k')/2+1/2)C_n}
       ≥   P{1/2n(λ_k-1-logλ_k)+o_p(n)>γ(p-k/2+1/2)C_n}
       →   1.



Suppose k'>k. Since d_i-1=O_p(√(1/n)) for i>k, by Taylor's expansion, we get

    logL̃_k-logL̃_k'   =   1/2n∑_i=k+1^k'(log d_i+1-d_i)
       =   1/2n∑_i=k+1^k'(log (1+(d_i-1))+1-d_i)
       =   -1/4n∑_i=k+1^k'(d_i-1)^2(1+o(1))
       =   -O_p(1).

Thus,

    P{ℓ̃(k)>ℓ̃(k')}   =   P{logL̃_k-logL̃_k'>γ(k-k')(p-(k+k')/2+1/2)C_n}
       >   P{-O_p(1)>γ(k-k')(p-(k+k')/2+1/2)C_n}
       →   1.








 Let Ł={l_1,⋯,l_n} and

    _n=Ł^⊤,

where  is orthogonal and its column _j is the eigenvector of Ł  with eigenvalues l_j.
Without loss of generality, we assume that E(_1(t))=. Then,

Let

    =_11^-1/2_1^⊤.

For 1≤ s≤ n, the s-th column of ,

    _s=(_11^-1/2_1^⊤)_s∼ N(, ),

and are independent and identically distributed.

Noting that

    (1/n_11^-1/2_1^⊤(+_n)X_1_11^-1/2)_ij   =   (1/n(_11^-1/2_1^⊤)(+Ł)(^⊤_1_11^-1/2))_ij
       =   (1/n(+Ł)^⊤)_ij
       =   1/n∑_s=1^n(1+l_s)z_isz_js,

we have

    1/n^̌⊤_11^-1/2_1^⊤(+_n)_1_11^-1/2=̌1/n∑_s=1^n(1+l_s)(∑_i,j=1^kv_iv_jz_isz_js),

where ^̌⊤=̌1. Hence,

    E(1/n_11^-1/2_1^⊤(+_n)_1_11^-1/2)_ij   =   1/n∑_s=1^n(1+l_s)E(z_isz_js)
       =   1/n∑_s=1^n(1+l_s)I_ij.

That is,

    E(1/n_11^-1/2_1^⊤(+_n)_1_11^-1/2)   =   1/n∑_s=1^n(1+l_s)=1/n(+_n).

Noting that,

    1/n^̌⊤(_11^-1/2_1^⊤(I+A_n)_1_11^-1/2-∑_s=1^n(1+l_s))   =   1/n∑_s=1^n(1+l_s)(∑_i,j=1^kv_iv_jz_isz_js)-1/n∑_s=1^n(1+l_s)^̌⊤   =   1/n∑_s=1^n(1+l_s)(∑_i,j=1^kv_iv_jz_isz_js)-1/n∑_s=1^n(1+l_s),

we have

    P(1/n^̌⊤(_11^-1/2_1^⊤(+_n)_1_11^-1/2-(+_n))>̌t)
       =   P(^̌⊤(_11^-1/2_1^⊤(+_n)_1_11^-1/2-∑_s=1^n(1+l_s))>̌nt)
       =   ∫_∏_s=1^ndF(z_1s,…,z_ks)
       ≤   ∫_exp(a∑_s=1^n(1+l_s)(∑_i,j=1^kv_iv_jz_isz_js)-a∑_s=1^n(1+l_s))/exp(ant)∏_s=1^ndF(z_1s,…,z_ks)
       ≤   ∫_exp(a∑_s=1^n(1+l_s)(∑_i,j=1^kv_iv_jz_isz_js)-a∑_s=1^n(1+l_s))/exp(ant)∏_s=1^ndF(z_1s,…,z_ks)
       =   exp(-ant-a∑_s=1^n(1+l_s))∏_s=1^nE(exp(a(1+l_s)∑_i,j=1^kv_iv_jz_isz_js))
       =   exp(-ant-a∑_s=1^n(1+l_s))∏_s=1^nE(exp(a(1+l_s)W_s)),

where 
    ={∑_s=1^n(1+l_s)(∑_i,j=1^kv_iv_jz_isz_js)-∑_s=1^n(1+l_s)>nt},


    W_s=∑_i,j=1^kv_iv_jz_isz_js=(∑_i=1^kv_iz_is)(∑_j=1^kv_jz_js)=(∑_i=1^kv_iz_is)^2∼χ^2(1).




    E(exp(a(1+l_s)W_s))   =   ∑_i=0^∞∫(a(1+l_s)w)^i/i!dF(w)
       ≤   1+a(1+l_s)+∑_i=2^∞∫(a(1+l_s)w)^i/i!dF(w)
       ≤   1+a(1+l_s)+∑_i=2^∞(a(1+l_s))^iE(w^i)/i!
       ≤   1+a(1+l_s)+∑_i=2^∞(2a(1+l_s))^i/i!Γ(1/2+i)/Γ(1/2)
       ≤   1+a(1+l_s)+∑_i=2^∞(2a(1+l_s))^i
       ≤   1+a(1+l_s)+4ρ_1a^2(1+l_s)^2,

where 2a(1+l_1)<1, ρ_1≥1/1-2a(1+l_1).

Hence, for any 0<a<1/2(1+l_1), we have

    ∏_s=1^nE(exp(a(1+l_s)W_s))   ≤   e^∑_s=1^nlog(1+a(1+l_s)+4ρ_1a^2(1+l_s)^2)
       ≤   e^a∑_s=1^n(1+l_s)+4ρ_1a^2∑_s=1^n(1+l_s)^2.


That is, for any 0<a<1/2(1+l_1), we have

    P(1/n^̌⊤(_11^-1/2X_1^⊤(+_n)_1_11^-1/2-(+_n))>̌t)
       ≤   exp(-ant+4ρ_1a^2∑_s=1^n(1+l_s)^2).



Let f(a)=-ant+4ρ_1a^2∑_s=1^n(1+l_s)^2.
Setting the derivative f'(a) to zero gives a=nt/8ρ_1∑_s=1^n(1+l_s)^2.
Thus,

    P(1/n^̌⊤(_11^-1/2_1^⊤(+_n)_1_11^-1/2-(+_n)I)>̌t)
       ≤   exp(-nt^2/16ρ_1∑_s=1^n(1+l_s)^2/n).








Furthermore, according to the proof of Lemma 3 in <cit.>, we have

    P(||1/n_11^-1/2_1^⊤(+_n)_1_11^-1/2-1/n(+_n)||>t)
       ≤   P(sup_1≤ j≤5^k1/n|_̌j^⊤(_11^-1/2_1^⊤(+_n)_1_11^-1/2-∑_s=1^n(1+l_s))_̌j|>t/4)
       ≤   ∑_j=1^5^k2P(1/n_̌j^⊤(_11^-1/2_1^⊤(+_n)_1_11^-1/2-∑_s=1^n(1+l_s))_̌j>t/4)
       ≤    2×5^kexp(-nt^2/16^2ρ_1∑_s=1^n(1+l_s)^2/n)
       =   2exp(Ck-nt^2/ρ∑_s=1^n(1+l_s)^2/n),

where _̌j^⊤_̌j=1, ρ=16^2ρ_1.

Noting that 2a(1+l_1)<1, we have

    t<4ρ_1∑_s=1^n(1+l_s)^2/n/1+l_1=ρ∑_s=1^n(1+l_s)^2/n/64(1+l_1).






By definition, each d_i solves the equation

    0=| d -_n|=| d -_22|| d - K_n(d)|,

where

    K_n(d)=1/n_1^⊤(+_̋n(d))_1,


    _̋n(d)=1/n_2(d -_22)^-1_2^⊤.

The d_i's then solve the determinant equation

    | d - K_n(d)|=0.

That is, d_i's are solutions of

    | d-1/n_1^⊤(+_̋n(d))_1|=0.

Or equivalently,

    | _11^-1d-1/n_11^-1/2_1^⊤(+_̋n(d))_1_11^-1/2|=0

Next, we demonstrate that the above determinant (<ref>) has k solutions in  the interval ((1+√(c))^2, ∞). For d>(1+√(c))^2,
by the definition of F_c(x), it is easy to see that m_1(d)=O(1) and m_2(d)=O(1).


Let β_j, j=1,…, p-k be the eigenvalues of _22=1/n_2^⊤_2. Then we have

    1/n(_̋n(d))   =1/n(1/n_2(d -_22)^-1_2^⊤)=1/n((d -_22)^-1_22)
       =1/n∑_j=1^p-kβ_j/d-β_j=p-k/n1/p-k∑_j=1^p-kβ_j/d-β_j
       =p-k/n∫x/d-xdF_n(x)=p-k/n∫x/d-xdF_c(x)(1+O_p(k/n))
       =(c+(p/n-c)-k/n)m_1(d)(1+O_p(k/n))
       =cm_1(d)(1+O_p(√(k/n))).



Moreover, we have

    1/n(_̋n^2(d))   =1/n((d -_22)^-1_22(d -_22)^-1_22)=1/n∑_j=1^p-kβ_j^2/(d-β_j)^2
       =p-k/n1/p-k∑_j=1^p-kβ_j^2/(d-β_j)^2=p-k/n∫x^2/(d-x)^2dF_n(x)
       =p-k/n∫x^2/(d-x)^2dF_c(x)(1+O_p(k/n))
       =(c+(p/n-c)-k/n)m_2(d)(1+O_p(k/n))
       =cm_2(d)(1+O_p(√(k/n))).


Combing (<ref>) and (<ref>), we have

    1/n∑_s=1^n(1+l_s)^2=1+2/n(_̋n(d))+1/n(_̋n^2(d))=1+O_p(1),

where l_1≥⋯≥ l_n are the eigenvalues of _̋n(d).

By Corollary <ref>, we have

    ||1/n_11^-1/2_1^⊤(+_̋n(d))_1_11^-1/2-1/n(+_̋n(d))I||=O_p(√(k/n)).


As a result, we have the following spectral decomposition of

    1/n_11^-1/2_1^⊤(+_̋n(d))_1_11^-1/2-1/n(+_̋n(d))
    
    =   (d)(v_1(d),…,v_k(d))^⊤(d),

where (d) is an orthogonal matrix and max_1≤ j≤ k|v_j(d)|=O_p(√(k/n)).

By (<ref>), we have

    1/n(+_̋n(d))=1+cm_1(d)(1+O_p(√(k/n))),

Combing (<ref>) and (<ref>), we have

    1/n_11^-1/2_1^⊤(+_̋n(d))_1_11^-1/2-(1+cm_1(d))
    
    =   (d)≜(d)(ṽ_1(d),…,ṽ_k(d))^⊤(d),

where max_1≤ j≤ k|ṽ_j(d)|=O_p(√(k/n)).

Note that

    _11^-1d-1/n_11^-1/2_1^⊤(+_̋n(d))_1_11^-1/2
       =   _11^-1d-(1+cm_1(d))-[1/n_11^-1/2_1^⊤(+_̋n(d))_1_11^-1/2-(1+cm_1(d))]
       =   _11^-1d-(1+cm_1(d))-(d).

and

    [ _11^-1=(1/λ_1,…,1/λ_k)^⊤. ]


Multiplying both sides of (<ref>) by ^⊤ from the left and by  from the right yields

    |^⊤(_11^-1d-1/n_11^-1/2_1^⊤(+_̋n(d))_1_11^-1/2)|
       =   |(d-(1+cm_1(d))λ_1/λ_1,…,d-(1+cm_1(d))λ_k/λ_k)-^⊤(d)|
       =   0.

Let ê_σ(i)'s be the eigenvalues of the following matrix,

    (d-(1+cm_1(d))λ_1/λ_1,…,d-(1+cm_1(d))λ_k/λ_k)-^⊤(d),

where σ(·) is a permutation  of the integers 1,…,n.


Then, d_i's are the solutions of

    ∏_i=1^kê_σ(i)=0


That is, for  any 1≤ i≤ k, d_i is the solution of

    ê_σ(i)=0.


Also, let e_i's be the eigenvalues of the following matrix,

    (d-(1+cm_1(d))λ_1/λ_1,…,d-(1+cm_1(d))λ_k/λ_k).

The following result shows how close are the eigenvalues of ê_σ(i)'s to e_i's.


Noting that

    ^⊤(d)=^⊤(d)(ṽ_1(d),…,ṽ_2(d))^⊤(d),


we have

    ((^⊤(d))^2)=((ṽ_1^2(d),…,ṽ_k^2(d)))=kO_p(k/n)=O_p(k^2/n).

By Lemma <ref>, we have

    ∑_i=1^k(ê_σ(i)-e_i)^2≤((^⊤(d))^2).

Hence,

    ∑_i=1^k|ê_σ(i)-e_i|≤√(k∑_i=1^k(ê_σ(i)-e_i)^2)=O_p(√(k^3/n)).


Combining (<ref>) and (<ref>), for  any 1≤ i≤ k, d_i tends to the solution of

    d_i-(1+cm_1(d_i))λ_i/λ_i=0.


According to <cit.>, for  any 1≤ i≤ k, we have

    d_i/ψ(λ_i)a.s.⟶1.






By definition, each d_i solves the equation

    0=| d -_n|=| d -_22|| d - K_n(d)|,

where

    K_n(d)=1/n_1^⊤(+_̋n(d))_1,


    _̋n(d)=1/n_2(d -_22)^-1_2^⊤.

The d_i's then solve the determinant equation

    | d - K_n(d)|=0.

That is

    | d_i-1/n_1^⊤(+_̋n(d_i))_1|=0.

Or equivalently,

    | _11^-1d_i-1/n_11^-1/2_1^⊤(+_̋n(d_i))_1_11^-1/2|=0

Note that

    Σ_11^-1d_i-1/n_11^-1/2_1^⊤(+_̋n(d_i))_1_11^-1/2
       =   Σ_11^-1d_i-(1+cm_1(d_i))-[1/n_11^-1/2_1^⊤(+_̋n(d_i))_1_11^-1/2-(1+cm_1(d_i))]
       =   _11^-1d_i-(1+cm_1(d_i))-(d_i),

where (d_i)=1/n_11^-1/2_1^⊤(+_̋n(d_i))_1_11^-1/2-(1+cm_1(d_i)).
Similar to the proof of (<ref>) in Theorem <ref>, we have

    (d_i)=(d_i)(ṽ_1(d_i),…,ṽ_k(d_i))^⊤(d_i),

where (d_i) is an orthogonal matrix and max_1≤ j≤ k|ṽ_j(d_i)|=O_p(√(k/n)).


Note that
 
    _11^-1=(1/λ_1,…,1/λ_k)^⊤.


Multiplying both sides of (<ref>) by ^⊤ from the left and by  from the right yields

    |^⊤(_11^-1d_i-1/n_11^-1/2_1^⊤(+_̋n(d_i))_1_11^-1/2)|
       =   |(d_i-(1+cm_1(d_i))λ_1/λ_1,…,d_i-(1+cm_1(d_i))λ_k/λ_k)-^⊤(d_i)|
       =   0.


Let

    (d_i)=(b_st)_k× k
    
    =   ^⊤(d_i)
    
    =   ^⊤(d_i)(ṽ_1(d_i),…,ṽ_k(d_i))^⊤(d_i)
    ≜   (d_i)(ṽ_1(d_i),…,ṽ_k(d_i))^⊤(d_i),

where (d_i)=^⊤(d_i)=(q_st)_k× k is an orthogonal matrix. For any 1≤ s≤ k, 1≤ t≤ k, and s≠ t, by the definitions of (d_i) and (d_i),   we have a_st=b_st.

Similar to the proof of (<ref>) in Theorem <ref>, we have

    ∑_s=1^k∑_t=1^kb_st^2=(^2(d_i))=O_p(k^2/n).


We also have

    ∑_t=1^kb_it^2=∑_t=1^k(∑_l=1^kṽ_l(d_i)q_ilq_tl)^2
    
    =   ∑_t=1^k∑_l=1^kṽ_l^2(d_i)q_il^2q_tl^2+∑_t=1^k∑_l=1^k∑_s≠ l^kṽ_l(d_i)ṽ_s(d_i)q_ilq_tlq_isq_ts
    
    =   ∑_t=1^k∑_l=1^kṽ_l^2(d_i)q_il^2q_tl^2+∑_l=1^k∑_s≠ l^kṽ_l(d_i)ṽ_s(d_i)q_ilq_is∑_t=1^kq_tlq_ts
    
    =   ∑_t=1^k∑_l=1^kṽ_l^2(d_i)q_il^2q_tl^2=∑_l=1^kṽ_l^2(d_i)q_il^2∑_t=1^kq_tl^2
    
    =   ∑_l=1^kṽ_l^2(d_i)q_il^2≤max_1≤ l≤ kṽ_l^2(d_i)∑_l=1^kq_il^2
    
    =   O_p(k/n).


Note that

    ((d_i))=∑_t=1^k(-1)^i+ta_itH_it=0,

where the (k-1)-th order determinant H_it is the cofactor of a_it, obtained by deleting the i-th row and the t-th column of the k-th order determinant ((d_i)).

Let

    H_ii=∏_j≠ iê_j'≜∏_j≠ i(e_j+ϵ_j'),

where ê_j'(j≠ i) are eigenvalues of corresponding matrix _̋ii.
Noting that _̋ii is a symmetric matrix, similar to the proof of (<ref>) in Theorem <ref>, we have

    ∑_j≠ i^k|ϵ_j'|=O_p(√(k^3/n)).


Noting that

    |∑_t≠ i1/e_tϵ_t'|   ≤   ∑_t≠ i1/|e_t||ϵ_t'|
       ≤   1/min_t≠ i|e_t|∑_t≠ i|ϵ_t'|
       ≤   1/min_t≠ i|e_t|O_p(√(k^3/n))
       =   O_p(√(k^3/n)),

and

    |∑_t≠ i∑_s≠ i,1≤ s<t1/e_te_sϵ_t'ϵ_s'|   ≤   ∑_t≠ i∑_s≠ i,1≤ s<t1/|e_te_s||ϵ_t'ϵ_s'|
       ≤   1/min_t≠ i,s≠ i,1≤ s<t|e_te_s|∑_t≠ i∑_s≠ i,1≤ s<t|ϵ_t'ϵ_s'|
       ≤   1/min_t≠ i,s≠ i,1≤ s<t|e_te_s|(∑_t≠ i|ϵ_t'|)^2
       =   1/min_t≠ i,s≠ i,1≤ s<t|e_te_s|O_p(k^3/n)
       =   O_p(k^3/n),

we have

    H_ii   =∏_j≠ i(e_j+ϵ_j')
       =∏_j≠ ie_j+∏_j≠ ie_j∑_t≠ i1/e_tϵ_t'+∏_j≠ ie_j∑_t≠ i∑_s≠ i,1≤ s<t1/e_te_sϵ_t'ϵ_s'+⋯
       =∏_j≠ ie_j+O_p(√(k^3/n))∏_j≠ ie_j+O_p(k^3/n)∏_j≠ ie_j+⋯
       =(1+O_p(√(k^3/n)))∏_j≠ ie_j.


For t≠ i,
 
    H_it=∏_j≠ i, t^kê_j”ê_t”,

where ê_j”(j≠ i, t) and ê_t” are eigenvalues of corresponding matrix _̋it.

By Gerschgorin's circle theorem, we have

    |H_it|=∏_j≠ i, t^k|ê_j'||ê_t”|≤∏_j≠ i, t^k(|e_j|+∑_r=1^k|b_jr|)|ê_t”|≤∏_j≠ i, t^k(|e_j|+∑_r=1^k|b_jr|)|∑_r≠ t^k|a_tr|.



Noting that

    ∑_s≠ i,t1/|e_s|∑_r=1^k|b_sr|   ≤   1/min_s≠ i,t|e_s|∑_s≠ i,t∑_r=1^k|b_sr|
       ≤   1/min_s≠ i,t|e_s|√(k^2∑_s≠ i,t∑_r=1^kb_sr^2)
       =   O_p(√(k^4/n)),

and

    ∑_s≠ i,t∑_u≠ i,t,1≤ u<s1/|e_se_u|∑_r=1^k|b_sr|∑_r=1^k|b_ur|
       ≤   1/min_s≠ i,t,u≠ i,t,1≤ u<s|e_se_u|∑_s≠ i,t∑_r=1^k|b_sr|∑_u≠ i,t,1≤ u<s∑_r=1^k|b_ur|
       =   1/min_s≠ i,t,u≠ i,t,1≤ u<s|e_se_u|O_p(k^4/n)
       =   O_p(k^4/n),

we have

    ∏_j≠ i, t^k(|e_j|+∑_r=1^k|b_jr|)
       =   ∏_j≠ i,t|e_j|+∏_j≠ i,t|e_j|∑_s≠ i,t1/|e_s|∑_r=1^k|b_sr|+∏_j≠ i,t|e_j|∑_s≠ i,t∑_u≠ i,t,1≤ u<s1/|e_se_u|∑_r=1^k|b_sr|∑_r=1^k|b_ur|+⋯
       =   ∏_j≠ i,t|e_j|+O_p(√(k^4/n))∏_j≠ i,t|e_j|+O_p(k^4/n)∏_j≠ i,t|e_j|+⋯
       =   (1+O_p(√(k^4/n)))∏_j≠ i,t|e_j|
       ≤   1+O_p(√(k^4/n))/min_j≠ i|e_j|∏_j≠ i|e_j|

Hence,

    |H_it|≤1+O_p(√(k^4/n))/min_j≠ i|e_j|∏_j≠ i|e_j||ê_t”|.

Combing (<ref>), (<ref>), (<ref>), and noting that for any 1≤ s≤ k, 1≤ t≤ k, s≠ t, a_st=b_st,
we have

    |∑_t≠ i^ka_itH_it|   ≤1+O_p(√(k^4/n))/min_j≠ i|e_j|∏_j≠ i^k|e_j|∑_t≠ i^k|a_itê_t”|
       ≤1+O_p(√(k^4/n))/min_j≠ i|e_j|∏_j≠ i^k|e_j|√(∑_t≠ i^ka_it^2∑_t≠ i^kê_t”^2)
       ≤1+O_p(√(k^4/n))/min_j≠ i|e_j|∏_j≠ i^k|e_j|√(∑_t≠ i^ka_it^2∑_t≠ i^k(∑_r≠ t^k|a_tr|)^2)
       ≤1+O_p(√(k^4/n))/min_j≠ i|e_j|∏_j≠ i^k|e_j|√(∑_t≠ i^kb_it^2∑_t≠ i^k(∑_r≠ t^k|b_tr|)^2)
       ≤1+O_p(√(k^4/n))/min_j≠ i|e_j|∏_j≠ i^k|e_j|√(∑_t≠ i^kb_it^2∑_t≠ i^kk∑_r≠ t^kb_tr^2)
       =1+O_p(√(k^4/n))/min_j≠ i|e_j|O_p(√(k^4/n^2))∏_j≠ i^k|e_j|.

Combing (<ref>),  (<ref>), and (<ref>), we have

    |a_ii|   =   |∑_t≠ i^ka_itH_it/H_ii|
       ≤   1+O_p(√(k^4/n))/1+O_p(√(k^3/n))1/min_j≠ i|e_j|O_p(√(k^4/n^2))
       =   1+O_p(√(k^4/n))/min_j≠ i|e_j|O_p(√(k^4/n^2))
       =   O_p(√(k^4/n^2)).



Thus,  for any 1≤ i≤ k, the i-th diagonal element of  satisfies

    d_i-(1+cm_1(d_i))λ_i/λ_i-[^⊤(d_i)]_ii=O_p(√(k^4/n^2)).







According to (<ref>) in Lemma <ref>, we have

    [^⊤(d_i)]_ii=b_ii=∑_l=1^kṽ_l(d_i)q_il^2=O_p(√(k/n)).

By Lemma <ref>, we have

    d_i-(1+cm_1(d_i))λ_i/λ_i=O_p(√(k/n)).

For  any 1≤ i≤ k, by Theorem <ref>, we have d_i>(1+√(c))^2. Let d_i=α_i+cα_i/α_i-1≜ψ(α_i).  Then, we get

    α_i=-c+d_i+1+√((c-d_i-1)^2-4d_i)/2,


    α̣_i/ḍ_i=1/2+1/2d_i-c-1/√((c-d_i-1)^2-4d_i)>0.

As a result, we have
α_i>1+√(c). According to <cit.>, m_1(d_i)=m_1(ψ(α_i))=1/α_i-1. Thus, (<ref>) becomes

    (α_i/λ_i-1)(1+c/α_i-1)=O_p(√(k/n)).

That is,

    α_i=λ_i(1+O_p(√(k/n))).

Thus, for  any 1≤ i≤ k, we have

    d_i-ψ(λ_i)/λ_i   =   α_i+cα_i/α_i-1-λ_i-cλ_i/λ_i-1/λ_i
       =   1+O_p(√(k/n))+c(1+O_p(√(k/n)))/λ_i(1+O_p(√(k/n)))-1-1-c/λ_i-1
       =   O_p(√(k/n))-cO_p(√(k/n))/(λ_i(1+O_p(√(k/n)))-1)(λ_i-1)
       =   O_p(√(k/n)).







By definition, each d_i solves the equation

    0=| d -_n|=| d -_22|| d - K_n(d)|,

where

    K_n(d)=1/n_1^⊤(+_̋n(d))_1,


    _̋n(d)=1/n_2(d -_22)^-1_2^⊤.

The d_i's then solve the determinant equation

    | d - K_n(d)|=0.

That is

    | d_i-1/n_1^⊤(+_̋n(d_i))_1|=0.


The random form K_n(d_i) can be decomposed as follows

    K_n(d_i)1/n_1^⊤(+_̋n(d_i))_1
    
    =   1/n(_1^⊤(+_̋n(d_i))_1-(+_̋n(d_i))_11)+1/n(+_̋n(d_i))_11
    
    =   1/√(n)R_n(d_i)+1/n(+_̋n(d_i))_11,

where

    R_n(d_i)=1/√(n)(_1^⊤(+_̋n(d_i))_1-(+_̋n(d_i))_11).

Let

    δ_i=√(n)(d_i-ψ(λ_i))/λ_i=O_p(√(k)).


Then the form d_i - K_n(d_i) can be decomposed as follows

    d_i - K_n(d_i)=ψ(λ_i)+1/√(n)δ_iλ_i- K_n(ψ(λ_i))-( K_n(d_i)- K_n(ψ(λ_i))).




Furthermore, using ^-1-^-1=^-1(-)^-1 and treating =(ψ(λ_i)+1/√(n)δ_iλ_i)-_22, =ψ(λ_i)-_22, we have

    _11^-1/2[ K_n(d_i)- K_n(ψ(λ_i))]_11^-1/2
    
    =   _11^-1/21/n_1^⊤1/n_2{[(ψ(λ_i)+1/√(n)δ_iλ_i)-_22]^-1
       -(ψ(λ_i)-_22)^-1}_2^⊤_1_11^-1/2
    
    =   -1/√(n)δ_iλ_i_11^-1/21/n_1^⊤1/n_2{[(ψ(λ_i)+1/√(n)δ_iλ_i)-_22]^-1
       [ψ(λ_i)-_22]^-1}_2^⊤_1_11^-1/2.

Treating

    _n=λ_i1/n_2[(ψ(λ_i)+1/√(n)δ_iλ_i)-_22]^-1(ψ(λ_i)-_22)^-1_2^⊤,



we have

    1/n(_n)
    
    =   1/n{λ_i1/n_2[(ψ(λ_i)+1/√(n)δ_iλ_i)-_22]^-1(ψ(λ_i)-_22)^-1_2^⊤}
    
    =   λ_i1/n{[(ψ(λ_i)+1/√(n)δ_iλ_i)-_22]^-1(ψ(λ_i)-_22)^-1_22}
    
    =   λ_i1/n∑_j=1^p-kβ_j/(ψ(λ_i)+1/√(n)δ_iλ_i-β_j)(ψ(λ_i)-β_j)
    
    =   λ_ip-k/n1/p-k∑_j=1^p-kβ_j/(ψ(λ_i)+1/√(n)δ_iλ_i-β_j)(ψ(λ_i)-β_j)
    
    =   λ_ip-k/n1/p-k∑_j=1^p-kβ_j/(ψ(λ_i)-β_j)^2(ψ(λ_i)-β_j)/(ψ(λ_i)+1/√(n)δ_iλ_i-β_j)
    
    =   λ_ip-k/n1/p-k∑_j=1^p-kβ_j/(ψ(λ_i)-β_j)^2(1+O_p(1/√(n)δ_i))
    
    =   λ_ip-k/n∫x/(ψ(λ_i)-x)^2dF_n(x)(1+O_p(1/√(n)δ_i))
    
    =   λ_ip-k/n∫x/(ψ(λ_i)-x)^2dF_c(x)(1+O_p(1/√(n)δ_i))(1+O_p(k/n))
    
    =   λ_i(c+(p/n-c)-k/n)m_3(ψ(λ_i))(1+O_p(1/√(n)δ_i))(1+O_p(k/n))
    
    =   λ_icm_3(ψ(λ_i))(1+O_p(√(k/n))),

where m_3(d)=∫x/(d-x)^2dF_c(x).

Noting that

    λ_im_3(ψ(λ_i))   =   λ_i∫x/(ψ(λ_i)-x)^2dF_c(x)
       =   ∫λ_i/(ψ(λ_i)-x)x/(ψ(λ_i)-x)dF_c(x)
       ≤   λ_i/ψ(λ_i)-(1+√(c))^2∫x/ψ(λ_i)-xdF_c(x)
       =   O(1),

we have

    1/n(_n)=O_p(1).



Note that,

    1/n(_n^2)
       =   1/n{λ_i1/n_2[(ψ(λ_i)+1/√(n)δ_iλ_i)-_22]^-1(ψ(λ_i)-_22)^-1_2^⊤
          λ_i1/n_2[(ψ(λ_i)+1/√(n)δ_iλ_i)-_22]^-1(ψ(λ_i)-_22)^-1_2^⊤}
       =   λ_i^21/n{[(ψ(λ_i)+1/√(n)δ_iλ_i)-_22]^-1(ψ(λ_i)-_22)^-1_22
          [(ψ(λ_i)+1/√(n)δ_iλ_i)-_22]^-1(ψ(λ_i)-_22)^-1_22}
       =   λ_i^21/n∑_j=1^p-kβ_j^2/(ψ(λ_i)+1/√(n)δ_iλ_i-β_j)^2(ψ(λ_i)-β_j)^2
       =   λ_i^2p-k/n1/p-k∑_j=1^p-kβ_j^2/(ψ(λ_i)+1/√(n)δ_iλ_i-β_j)^2(ψ(λ_i)-β_j)^2
       =   λ_i^2p-k/n∫x^2/(ψ(λ_i)-x)^4dF_n(x)(1+O_p(1/√(n)δ_i))
       =   λ_i^2p-k/n∫x^2/(ψ(λ_i)-x)^4dF_c(x)(1+O_p(1/√(n)δ_i))(1+O_p(k/n))
       =   λ_i^2(c+(p/n-c)-k/n)∫x^2/(ψ(λ_i)-x)^4dF_c(x)(1+O_p(1/√(n)δ_i))(1+O_p(k/n))
       =   λ_i^2 cm_4(ψ(λ_i))(1+O_p(√(k/n))),

where m_4(d)=∫x^2/(d-x)^4dF_c(x).

Noting that

    λ_i^2m_4(ψ(λ_i))   =   λ_i^2∫x^2/(ψ(λ_i)-x)^4dF_c(x)
       =   ∫λ_i^2/(ψ(λ_i)-x)^2x^2/(ψ(λ_i)-x)^2dF_c(x)
       ≤   λ_i^2/(ψ(λ_i)-(1+√(c))^2)^2∫x^2/(ψ(λ_i)-x)^2dF_c(x)
       =   O(1),

we have

    1/n(_n^2)=O_p(1).



Hence, we have

    1/n∑_s=1^nl_s^2=1/n(_n^2)=O_p(1).

where l_1≥⋯≥ l_n are the eigenvalues of _n.


By Lemma <ref>, we have

    ||1/n_11^-1/2_1^⊤_n_1_11^-1/2-1/n(_n)||=O_p(√(k/n)).

Combining (<ref>), (<ref>), (<ref>), and (<ref>), we have the following spectral decomposition of

    _11^-1/2[ K_n(d_i)- K_n(ψ(λ_i))]_11^-1/2
    
    =   -1/√(n)δ_iλ_icm_3(ψ(λ_i))(1+O_p(√(k/n)))
       -1/√(n)δ_i(d_i)(w_1(d_i),…,w_k(d_i))^⊤(d_i),

where (d_i) is orthogonal and max_1≤ j≤ k|w_j(d_i)|=O_p(√(k/n)).

Noting that p/n-c=o(n^-1/2), similar to the proof of (<ref>) in Theorem <ref>, we have

    1/n(+_̋n(d_i))=1+cm_1(ψ(λ_i))(1+o(n^-1/2)).


Combining (<ref>), (<ref>), (<ref>), and (<ref>), we have


    _11^-1/2(d_i - K_n(d_i))_11^-1/2
    
    =   _11^-1ψ(λ_i)+1/√(n)δ_iλ_i_11^-1-[1+cm_1(ψ(λ_i))(1+o(n^-1/2))]
       -1/√(n)_11^-1/2 R_n(ψ(λ_i))_11^-1/2+1/√(n)δ_iλ_icm_3(ψ(λ_i))(1+O_p(√(k/n)))
       +1/√(n)δ_i(d_i)(w_1(d_i),…,w_k(d_i))^⊤(d_i).



Note that

    ψ(λ_i)-(1+cm_1(ψ(λ_i)))λ_j/λ_j=(λ_i/λ_j-1)(1+c/λ_i-1),

and

    [ _11^-1=(1/λ_1,…,1/λ_k)^⊤. ]


Multiplying both sides of (<ref>) by ^⊤ from the left and by  from the right, we have

    |^⊤_11^-1/2(d_i - K_n(d_i))_11^-1/2|
       =   |((λ_i/λ_1-1)(1+c/λ_i-1),…,(λ_i/λ_k-1)(1+c/λ_i-1))
          +cm_1(ψ(λ_i))o(n^-1/2)+1/√(n)δ_iλ_i(1/λ_1,…,1/λ_k)
          -1/√(n)^⊤_11^-1/2 R_n(ψ(λ_i))_11^-1/2
          +1/√(n)δ_iλ_icm_3(ψ(λ_i))(1+O_p(√(k/n)))
          -1/√(n)δ_i^⊤(d_i)(w_1(d_i),…,w_k(d_i))^⊤(d_i)|.


Note that ^⊤(d_i) is an orthogonal matrix and max_1≤ j≤ k|w_j(d_i)|=O_p(√(k/n)).
The i-th diagonal element of

    1/√(n)δ_i^⊤(d_i)(w_1(d_i),…,w_k(d_i))^⊤(d_i)

is 1/√(n)δ_iO_p(√(k/n)). Since (λ_i/λ_i-1)(1+c/λ_i-1)=0, by Lemma <ref>, the i-th diagonal element of

    (d_i)=^⊤_11^-1/2(d_i - K_n(d_i))_11^-1/2

satisfies

    1/√(n)δ_i(1+cm_3(ψ(λ_i))λ_i)-[1/√(n)^⊤_11^-1/2 R_n(ψ(λ_i))_11^-1/2]_ii
          +1/√(n)δ_iO_p(√(k/n))+o(n^-1/2)=O_p(√(k^4/n^2)).




That is,

    1/√(n)δ_i(1+cm_3(ψ(λ_i))λ_i)-[1/√(n)^⊤_11^-1/2 R_n(ψ(λ_i))_11^-1/2]_ii
          +1/√(n)δ_iO_p(√(k/n))+o(n^-1/2)+O_p(√(k^4/n^2))=0.


Or equivalently,

    |δ_i(1+cm_3(ψ(λ_i))λ_i)-[^⊤_11^-1/2 R_n(ψ(λ_i))_11^-1/2]_ii
          +δ_iO_p(√(k/n))+o(1)+O_p(√(k^4/n))|=0.


Noting that

    δ_i=√(n)(d_i-ψ(λ_i))/λ_i=O_p(√(k)),

under the condition k=o(n^1/4), δ_i tends to a solution of

    |δ_i(1+cm_3(ψ(λ_i))λ_i)-[^⊤_11^-1/2 R_n(ψ(λ_i))Σ_11^-1/2]_ii|=0.



That is, for any 1≤ i≤ k,

    δ_ia.s.⟶(1+cm_3(ψ(λ_i))λ_i)^-1[^⊤_11^-1/2 R_n(ψ(λ_i))_11^-1/2]_ii.







Similar to the proof of Corollary <ref>, by Chebyshev's inequality, we have

    P(||1/n_11^-1/2_1^⊤(+_n)_1_11^-1/2-1/n(+_n)I||>t)
       ≤   P(max_1≤ i≤ k∑_1≤ j≤ k|1/n∑_s=1^n(1+l_s)z_isz_js-1/n∑_s=1^n(1+l_s)I_ij|>t)
       ≤   ∑_1≤ i≤ kP(∑_1≤ j≤ k∑_s=1^n(1+l_s)|z_isz_js-I_ij|>nt)
       =   ∑_1≤ i,j≤ k∑_1≤ j≤ k∑_s=1^n(1+l_s)^2E(z_is^2z_js^2)/n^2t^2
       =   ∑_1≤ i,j≤ k∑_s=1^n(1+l_s)^2E(z_is^2z_js^2)/n^2t^2
       ≤   ∑_1≤ i,j≤ k∑_s=1^n(1+l_s)^2√(E(z_is^4)E(z_js^4))/n^2t^2
       ≤   Ck^2∑_s=1^n(1+l_s)^2/n^2t^2.


Thus, we have

    ||1/n_11^-1/2_1^⊤(+_n)_1_11^-1/2-1/n(+_n)||=√(1/n∑_s=1^n(1+l_s)^2)O_p(k/√(n)).







Let log L_k'=-1/2n∑_i=1^k'log d_i-1/2n(p-k')logλ̂_k'. According to <cit.>,

    λ̂=1/(p-k)∑_i=k+1^pd_ia.s.⟶1.



Suppose k'<k. We have

    ℓ(k)-ℓ(k')=log L_k-log L_k'-γ(k-k')(p-(k+k')/2+1/2).


By Theorem <ref>, for k'<i≤ k,

    d_i/ψ(λ_i)a.s.⟶ 1.


Since ψ(λ_k)-1-logψ(λ_k)>2γ c, by Taylor's expansion, we get

    P{ℓ(k)>ℓ(k')}
       =   P{log L_k-log L_k'>γ(k-k')(p-(k+k')/2+1/2)}
       =   P{-1/2n∑_i=1^klog d_i-1/2n(p-k)logλ̂+1/2n∑_i=1^k'log d_i+1/2n(p-k')logλ̂_k'
          >γ(k-k')(p-(k+k')/2+1/2)}
       =   P{-1/2n(p-k)logλ̂+1/2n(p-k')logλ̂_k'-1/2n∑_i=k'+1^klog d_i
          >γ(k-k')(p-(k+k')/2+1/2)}
       =   P{-1/2n(p-k')logλ̂+1/2n(p-k')logλ̂_k'+1/2n(k-k')logλ̂
          -1/2n∑_i=k'+1^klog d_i>γ(k-k')(p-(k+k')/2+1/2)}
       =   P{1/2n(p-k')log(λ̂_k'/λ̂)+1/2n(k-k')logλ̂
          -1/2n∑_i=k'+1^klog d_i>γ(k-k')(p-(k+k')/2+1/2)}
       =   P{1/2n(p-k')log[(1-(k-k')/(p-k'))(∑_i=k'+1^pd_i/∑_i=k+1^pd_i)]+1/2n(k-k')logλ̂
          -1/2n∑_i=k'+1^klog d_i>γ(k-k')(p-(k+k')/2+1/2)}
       =   P{-1/2n(k-k')+1/2n(p-k')log[(1/p-k∑_i=k'+1^kd_i+1/p-k∑_i=k+1^pd_i)/(1/p-k∑_i=k+1^pd_i)]
          +o(n(k-k'))-1/2n∑_i=k'+1^klog d_i>γ(k-k')(p-(k+k')/2+1/2)}.



For the case 1/p-k∑_i=k'+1^kψ(λ_i)=o(1), noting that
1/p-k∑_i=k'+1^kd_i=o_p(1), we have

    P{ℓ(k)>ℓ(k')}
       =   P{-1/2n(k-k')+1/2n(p-k')log[(1/p-k∑_i=k'+1^kd_i+1/p-k∑_i=k+1^pd_i)/(1/p-k∑_i=k+1^pd_i)]
          +o(n(k-k'))-1/2n∑_i=k'+1^klog d_i>γ(k-k')(p-(k+k')/2+1/2)}
       =   P{-1/2n(k-k')+1/2n(p-k')[(1/p-k∑_i=k'+1^kd_i)/(1/p-k∑_i=k+1^pd_i)]+o(n(k-k'))
          -1/2n∑_i=k'+1^klog d_i>γ(k-k')(p-(k+k')/2+1/2)}
       =   P{1/2n∑_i=k'+1^k(d_i-1-log d_i)+o(n(k-k'))
          >γ(k-k')(p-(k+k')/2+1/2)}
       =   P{1/2n∑_i=k'+1^k(ψ(λ_i)-1-logψ(λ_i))>γ(k-k')(p-(k+k')/2+1/2)(1+o(1))}
       ≥   P{ψ(λ_k)-1-logψ(λ_k)>2γ(p-(k+k')/2+1/2)(1+o(1))/n}
       ≥   P{ψ(λ_k)-1-logψ(λ_k)>2γ p(1+o(1))/n}
       →   1.


For the case 1/p-k∑_i=k'+1^kψ(λ_i)/C → 1, where C is a constant, noting that
1/p-k∑_i=k'+1^kd_i/Ca.s.⟶ 1, we have

    ∏_i=k'+1^kd_i≤ (1/k-k'∑_i=k'+1^kd_i)^k-k',


    ∑_i=k'+1^klog d_i≤(k-k')log(1/k-k'∑_i=k'+1^kd_i)=O_p(klog p),

we have

    P{ℓ(k)>ℓ(k')}
       =   P{-1/2n(k-k')+1/2n(p-k')log[(1/p-k∑_i=k'+1^kd_i+1/p-k∑_i=k+1^pd_i)/(1/p-k∑_i=k+1^pd_i)]
          +o(n(k-k'))-1/2n∑_i=k'+1^klog d_i>γ(k-k')(p-(k+k')/2+1/2)}
       ≥   P{-1/2n(k-k')+1/2n(p-k')log C+o(n(k-k'))
          -1/2n∑_i=k'+1^klog d_i>γ(k-k')(p-(k+k')/2+1/2)}
       →   1.



For the case 1/p-k∑_i=k'+1^kd_i/C_n→1, where C_n→∞, noting that
1/p-k∑_i=k'+1^kd_i/C_na.s.⟶1, we have

    ∑_i=k'+1^klog d_i≤(k-k')log(1/k-k'∑_i=k'+1^kd_i)=O_p(k(log C_n+log p))

we have

    P{ℓ(k)>ℓ(k')}
       =   P{-1/2n(k-k')+1/2n(p-k')log[(1/p-k∑_i=k'+1^kd_i+1/p-k∑_i=k+1^pd_i)/(1/p-k∑_i=k+1^pd_i)]
          +o(n(k-k'))-1/2n∑_i=k'+1^klog d_i>γ(k-k')(p-(k+k')/2+1/2)}
       ≥   P{-1/2n(k-k')+1/2n(p-k')log C_n+o(n(k-k'))
          -1/2n∑_i=k'+1^klog d_i>γ(k-k')(p-(k+k')/2+1/2)}
       →   1.


Suppose k'>k. For k<i≤ k'=o(p),

    d_ia.s.⟶μ_1.

According to <cit.>, μ_1=(1+√(c))^2.

Since γ>1/2+√(1/c)-log(1+√(c))/c, by Taylor's expansion, we get

    P{ℓ(k)>ℓ(k')}
       =   P{log L_k-log L_k'>γ(k-k')(p-(k+k')/2+1/2)}
       =   P{-1/2n∑_i=1^klog d_i-1/2n(p-k)logλ̂+1/2n∑_i=1^k'log d_i+1/2n(p-k')logλ̂_k'
          >-γ(k'-k)(p-(k'+k)/2+1/2)}
       =   P{-1/2n(p-k)logλ̂+1/2n(p-k')logλ̂_k'+1/2n∑_i=k+1^k'log d_i
          >-γ(k'-k)(p-(k'+k)/2+1/2)}
       =   P{-1/2n(p-k')logλ̂+1/2n(p-k')logλ̂_k'-1/2n(k'-k)logλ̂
          +1/2n∑_i=k+1^k'log d_i>-γ(k'-k)(p-(k'+k)/2+1/2)}
       =   P{1/2n(p-k')log(λ̂_k'/λ̂)-1/2n(k'-k)logλ̂
          +1/2n∑_i=k+1^k'log d_i>-γ(k'-k)(p-(k'+k)/2+1/2)}
       =   P{1/2n(p-k')log[(1+(k'-k)/(p-k'))(∑_i=k'+1^pd_i/∑_i=k+1^pd_i)]-1/2n(k'-k)logλ̂
          +1/2n∑_i=k+1^k'log d_i>-γ(k'-k)(p-(k'+k)/2+1/2)}
       =   P{1/2n(k'-k)+1/2n(p-k')log[(1/p-k∑_i=k+1^pd_i-1/p-k∑_i=k+1^k'd_i)/(1/p-k∑_i=k+1^pd_i)]
          -o(n(k'-k))+1/2n∑_i=k+1^k'log d_i>-γ(k'-k)(p-(k'+k)/2+1/2)},




    P{ℓ(k)>ℓ(k')}
       =   P{1/2n(k'-k)-1/2n(p-k')[(1/p-k∑_i=k+1^k'd_i)/(1/p-k∑_i=k+1^pd_i)]-o(n(k-k'))
          +1/2n∑_i=k+1^k'log d_i>-γ(k'-k)(p-(k'+k)/2+1/2)}
       =   P{-1/2n∑_i=k+1^k'(d_i-1-log d_i)>-γ(k'-k)(p-(k'+k)/2+1/2)(1+o(1))}
       =   P{-1/2n((1+√(c))^2-1-2log(1+√(c)))>-γ(p-(k'+k)/2+1/2)(1+o(1))}
       =   P{γ>n/(p-(k'+k)/2+1/2)(c/2+√(c)-log(1+√(c))(1+o(1))}
       =   P{γ>p/(p-(k'+k)/2+1/2)n/p(c/2+√(c)-log(1+√(c)))(1+o(1))}
       →   1.




ims
=14.5pt





