On high-dimensional wavelet eigenanalysis
P.A. and H.W. were partially supported by the French ANR AMATIS 2011 grant #ANR-11-BS01-0011. G.D.'s long term visits to ENS de Lyon were supported by the school, the CNRS and the Carol Lavin Bernick faculty grant.
AMS Subject classification. Primary: 60G18, 60B20, 42C40. Secondary: 62H25.
Keywords and phrases: wavelets, operator self-similarity, random matrices.

Patrice Abry
Univ Lyon, ENS de Lyon,
Univ Claude Bernard,
Laboratoire de Physique,
F-69342 Lyon, France
B. Cooper Boniece
Department of Mathematics and Statistics
Washington University in St. Louis
Gustavo Didier
Mathematics Department
Tulane University
Herwig Wendt
IRIT-ENSEEIHT, CNRS (UMR 5505),
Université de Toulouse, France

In this paper, we characterize the asymptotic and large scale behavior of the eigenvalues of wavelet random matrices in high dimensions. We assume that possibly non-Gaussian, finite-variance $p$-variate measurements are made of a low-dimensional $r$-variate ($r \ll p$) fractional stochastic process with non-canonical scaling coordinates and in the presence of additive high-dimensional noise. The measurements are correlated both time-wise and between rows. We show that the $r$ largest eigenvalues of the wavelet random matrices, when appropriately rescaled, converge to scale invariant functions in the high-dimensional limit. By contrast, the remaining $p-r$ eigenvalues remain bounded.  Under additional assumptions, we show that, up to a log transformation, the $r$ largest eigenvalues of wavelet random matrices exhibit asymptotically Gaussian distributions. The results have direct consequences for statistical inference.

§ INTRODUCTION

A wavelet is a unit $L^2(\bbR)$-norm function that annihilates polynomials. For a fixed (octave) $j \in \bbN$, a wavelet random matrix is given by

\begin{equation}\label{e:W(a2^j)_intro}
{\mathbf W}(a(n)2^j) = \frac{1}{n_{a,j}}\sum^{n_{a,j}}_{k=1} D(a(n)2^j,k)D(a(n)2^j,k)^*.
\end{equation}

In (<ref>), $^*$ denotes transposition, $n_{a,j} = n/a(n)2^j$ is the number of wavelet-domain observations for a sample size $n$, and each random vector $D(a(n)2^j,k)$ is the wavelet transform of a multivariate stochastic process $Y$ at the dyadic scale $a(n)2^j$ and shift $k \in \bbZ$. The entries of $\{D(a(n)2^j,k)\}_{k \in \bbN}$ are generally correlated. The so-named wavelet eigenanalysis methodology consists in using the behavior across scales of the eigenvalues of wavelet random matrices to study the fractality of stochastic systems (Abry and Didier <cit.>). In this paper, we characterize the asymptotic and large scale behavior of the eigenvalues of wavelet random matrices in high dimensions. The (possibly non-Gaussian) underlying stochastic process $Y = \{Y(t)\}_{t \in \bbZ}$ is assumed to have the form

\begin{equation}\label{e:Y(t)}
Y(t) = P X(t) + Z(t), \quad t \in \bbZ.
\end{equation}

In (<ref>), both $Y$ and the noise term $Z = \{Z(t)\}_{t \in \bbZ}$ are (high-dimensional) $p=p(n)$-variate processes, $P=P(n)$ is a rectangular coordinates matrix and, for fixed $r \in \bbN$, $X = \{X(t)\}_{t \in \bbZ}$ is a (low-dimensional) $r$-variate fractional process. In particular, the measurements $Y$ are correlated time-wise and between rows. We show that, if the ratio $a(n)p(n)/n$ converges to a positive constant, under appropriate rescaling the $r$ largest eigenvalues of ${\mathbf W}(a(n)2^j)$ converge in probability to scale invariant functions.  By contrast,  the remaining $p(n)-r$ eigenvalues remain bounded in probability. In addition, we show that, up to a log transformation, the $r$ largest eigenvalues of ${\mathbf W}(a(n)2^j)$ exhibit asymptotically Gaussian distributions. The results bear direct consequences for statistical inference starting from high-dimensional measurements of the form of a signal-plus-noise system (<ref>), where $X$ is a latent process containing fractal information and $P$ (as well as $Z$) is unknown.

Since the 1950s, the spectral behavior of large dimensional random matrices has attracted considerable attention from the mathematical research community. In quantum mechanics, for example, random matrices are of great interest as statistical mechanical models of infinite dimensional and possibly unknown Hamiltonian operators (e.g., Mehta and Gaudin <cit.>, Dyson <cit.>, Arous and Guionnet <cit.>, Soshnikov <cit.>, Mehta <cit.>, Deift <cit.>, Anderson et al. <cit.>, Tao and Vu <cit.>, Erdős et al. <cit.>). Random matrices have also naturally emerged as one essential mathematical framework for the modern era of “Big Data" (Brody <cit.>), when hundreds to several tens of thousands of time series get recorded and stored on a daily basis. In arbitrary dimension $p$, stochastic modeling has to cope with the curse of dimensionality: as the dimension increases, the amount of data needed for traditional statistical inference typically grows exponentially (Giraud <cit.>, Yao et al. <cit.>, Wainwright <cit.>). Consequently, one is often interested in understanding the behavior of random constructs such as the spectral distribution of sample covariance matrices when the dimension $p$ is comparable to the sample size $n$ (e.g., Tao and Vu <cit.>, Xia et al. <cit.>, Paul and Aue <cit.>).

In turn, scale invariance manifests itself in a wide range of natural and social phenomena such as in climate studies (Isotta et al. <cit.>), critical phenomena (Sornette <cit.>), dendrochronology (Bai and Taqqu <cit.>), hydrology (Benson et al. <cit.>) and turbulence (Kolmogorov <cit.>). In a multidimensional setting, scaling behavior does not always appear along standard coordinate axes, and often involves multiple (scaling) relations. A $\bbR^r$-valued stochastic process $X$ is called operator self-similar (o.s.s.; Laha and Rohatgi <cit.>, Hudson and Mason <cit.>) if it exhibits the scaling property

\begin{equation}\label{e:def_ss}
\{X(ct)\}_{t\in\bbR} \stackrel {\textnormal{f.d.d.}}{=} \{c^H X(t)\}_{t\in\bbR}, \quad c>0.
\end{equation}

In (<ref>), $H$ is some (Hurst) matrix $H$ whose eigenvalues have real parts lying in the interval $(0,1]$ and $c^H := \exp\{\log(c) H\} = \sum^{\infty}_{k=0} \frac{(\log(c) H)^k}{k!}$. A canonical model for multivariate fractional systems is operator fractional Brownian motion (ofBm), namely, a Gaussian, o.s.s., stationary-increment stochastic process (Maejima and Mason <cit.>, Mason and Xiao <cit.>, Didier and Pipiras <cit.>). In particular, ofBm is the natural multivariate generalization of the classical fBm (Mandelbrot and Van Ness <cit.>).

The literature on random matrices under dependence as well as on high-dimensional stochastic processes has been expanding at a fast pace (e.g., Basu and Michailidis <cit.>, Chakrabarty et al. <cit.>, Merlevède and Peligrad <cit.>, Che <cit.>, Steland and von Sachs <cit.>, Taylor and Salhi <cit.>, Wang et al. <cit.>, Zhang and Wu <cit.>, Erdős et al. <cit.>, Merlevède et al. <cit.>, Bourguin et al. <cit.>). In applications, models of the form (<ref>) appear, for example, in neuroscience and fMRI imaging (Ciuciu et al. <cit.>, Liu et al. <cit.>) and in econometrics (e.g., Brown <cit.>, Lam and Yao <cit.>).

In the characterization of scaling properties, the use of eigenanalysis was first proposed in Meerschaert and Scheffler <cit.> for operator stable laws, and later in Becker-Kern and Pap <cit.> for o.s.s. processes in the time domain. It has also been used in the cointegration literature (e.g., Phillips and Ouliaris <cit.>, Li et al. <cit.>, Zhang et al. <cit.>). In Abry and Didier <cit.>, wavelet eigenanalysis is proposed in the construction of a general methodology for the statistical identification of the scaling (Hurst) structure of ofBm in low dimensions.

In Abry et al. <cit.> and Boniece et al. <cit.>, presented without proofs and containing preliminary computational studies, wavelet random matrices were first used in the modeling of high-dimensional systems. In this paper, we construct the mathematical foundations of wavelet eigenanalysis in high dimensions by investigating the properties of the eigenvalues of large wavelet random matrices. We assume measurements given by (<ref>), where the fractional behavior of $X$ is characterized by a scaling matrix of the Jordan form

\begin{equation}\label{e:H=PHdiag(h1,...,hn)P^(-1)H}
H = P_H \hspace{0.5mm}\textnormal{diag}(h_1,\hdots,h_r)\hspace{0.5mm}P^{-1}_{H},
\end{equation}

with real eigenvalues. The measurements $Y$ display correlation time-wise and between rows. In our main results (Section <ref>), we consider the behavior of the eigenvalues $\lambda_{1}\big({\mathbf W}(a(n)2^j)\big) \leq \hdots \leq \lambda_{p(n)}\big({\mathbf W}(a(n)2^j)\big)$ of wavelet random matrices in the three-fold limit as the sample size ($n$), dimension ($p(n)$) and scale ($a(n)$) go to infinity. We show that the $r$ largest eigenvalues of wavelet random matrices display fractal – or scaling – properties determined by (<ref>) as well as (under assumptions) asymptotically Gaussian fluctuations. In particular, such eigenvalues are explosive. By contrast, the remaining $p(n)-r$ eigenvalues do not exhibit fractality and remain bounded.

To be more precise, under very general assumptions, we establish that, for positive functions $\xi_q(\cdot)$,

\begin{equation}\label{e:eigenvalue_scaling_intro}
\lambda_{p(n)-r+q}\Big(\frac{{\mathbf W}(a(n)2^j)}{a(n)^{2h_{q}+1}} \Big) \stackrel{\bbP}\rightarrow \xi_{q}(2^j), \quad q = 1,\hdots,r,
\end{equation}

whereas $\lambda_{\ell}\big({\mathbf W}(a(n)2^j)\big)$, $\ell = 1,\hdots,p(n)-r$, are bounded in probability (Theorem <ref>; see Figure <ref> for an illustration). Moreover, under slightly stronger conditions, we show that the random vector

\begin{equation}\label{e:log-eigenvalue_normality_intro}
\sqrt{\frac{n}{a(n)2^j}}\Big(\hspace{1mm}\lambda_{p(n)-r+q}\Big(\frac{{\mathbf W}(a(n)2^j)}{a(n)^{2h_{q}+1}} \Big) - \lambda_{p(n)-r+q}\Big(\frac{\bbE{\mathbf W}(a(n)2^j)}{a(n)^{2h_{q}+1}} \Big) \hspace{1mm}\Big), \quad q = 1,\hdots, r,
\end{equation}

is asymptotically Gaussian (Theorem <ref>). In (<ref>) and (<ref>), we assume

\begin{equation}\label{e:three-fold_lim}
\lim_{n \rightarrow \infty} \frac{p(n) a(n)}{n} =: c \in [0,\infty).
\end{equation}

In particular, the convergence rate in (<ref>) also involves the scaling limit. Note that this stands in sharp contrast with traditional high-dimensional analysis of sample covariance matrices, in which one considers the ratio $\lim_{n \rightarrow \infty} p(n)/n$ and the largest eigenvalue often exhibits universality in the form of Tracy–Widom fluctuations (e.g., Bai and Silverstein <cit.>, Lee and Schnelli <cit.>).

To the best of our knowledge, this paper provides the first mathematical study of the high-dimensional properties of wavelet random matrices. Characterizing the fractality of wavelet eigenvalues requires specific arguments building upon the Courant-Fischer variational formulation of eigenvalues. In turn, characterizing their fluctuations involves handling Taylor expansions of wavelet log-eigenvalues both in the high-dimensional limit and in the presence of explosive terms, the latter due to the discrepant scaling rates in (<ref>). For the sake of clarity and mathematical generality, our assumptions are stated directly in the wavelet domain (see Section <ref>). We further provide encompassing classes of examples covered by the assumptions. Supposing $X$ is an ofBm, we describe a very broad Gaussian framework for the noise term $Z$ satisfying the assumptions used in Section <ref>. To this end, we draw upon the concentration of measure phenomenon (Ledoux <cit.>, Boucheron et al. <cit.>) to nonasymptotically control the behavior of large dimensional wavelet random matrices formed from colored noise (see Section <ref>). We also discuss some simple finite-variance and non-Gaussian instances $X$ of interest, hence illustrating the broad scope of the assumptions (see Section <ref>; for comments on non-Gaussian $Z$, see Example <ref>).

Our results have direct consequences for the empirical identification and description of fractality in high-dimensional physical systems. Non-canonical scaling coordinates $P(n)P_H$ generally mix together slow and fast scaling laws in fractional stochastic processes, leading to the so-called amplitude and dominance effects (see Abry and Didier <cit.> for a detailed discussion). These effects manifest themselves in the form of strong biases in standard, univariate-like statistical methodology when applied to measurements of multidimensional phenomena such as Internet traffic (e.g., Abry and Didier <cit.>, Section 6), cointegration (see, for instance, Kaufmann and Stern <cit.>, Schmith et al. <cit.> on climate science) and systems modeled in blind source separation problems (e.g., Comon and Jutten <cit.>; see also Example <ref>). As a consequence of the results in this paper, the multidimensional fractality of the latent process $X$ can be statistically characterized by means of the consistent estimation the scaling eigenvalues $h_1,\hdots,h_r$ (see (<ref>)) even in the presence of high-dimensional noise and of non-canonical coordinates $P(n)P_H$ (see Section <ref> for a discussion).

This paper is organized as follows. In Section <ref>, we provide the basic wavelet framework, definitions and wavelet-domain assumptions used throughout the paper. In Section <ref>, we provide the main results on the asymptotic and high-dimensional behavior of wavelet eigenvalues. In Section <ref>, we provide a Gaussian framework, as well as Gaussian and non-Gaussian examples. In Section <ref>, we lay out conclusions and discuss several open problems that this work leads to, including new aspects of the theory of wavelet random matrices, as well as consequences for statistical inference and modeling. All proofs can be found in the appendix.

[The high-dimensional asymptotic behavior of log-eigenvalues]The high-dimensional asymptotic behavior of wavelet log-eigenvalues. For $r=6$ and $h_i \in \{0.1,0.3, 0.5,0.6,0.8,0.9\}$, the plots display the asymptotic behavior of $\frac{1}{2}[(\log \Lambda_q(2^j))/j-1]$, where $\Lambda_q(2^j):=\lambda_q(\mathbf{W}(2^j))$. In all plots, $p$ and $n_j = n/2^j$ increase while their ratio remains fixed at ${p/n_j}= c = 1/2$ (left column) and $c=1/4$ (right column). The coordinates matrix $P P_H$ is randomly drawn for each $n$. The $r=6$ largest wavelet log-eigenvalues approach $h_1,\hdots,h_6$ as $n$, $p$ and $j$ grow, whereas the remaining wavelet log-eigenvalues (see (<ref>)) tend toward zero as $j$ increases in all instances. With $c$ fixed (along each column), the larger sample size $n$ (bottom row) results in near-convergence at smaller octaves $j$ (n.b.: axes have been shifted to align the curves). Similarly, as $c$ decreases (between left and right columns), the near-convergence also occurs at smaller $j$.

§ FRAMEWORK AND NOTATION

For $m \in \bbN$, let $M(m,\bbR)$ and $M(m,\bbC)$ be the spaces of $m \times m$ real- and complex-valued matrices, respectively. Also, let $M(m_1,m_2,\bbR)$ be the spaces of $m_1 \times m_2$ real-valued matrices. Let ${\mathcal S}(m,\bbC)$, ${\mathcal S}_{\geq 0}(m,\bbC)$ and ${\mathcal S}_{> 0}(m,\bbC)$, respectively, be the space of $m \times m$ Hermitian matrices and the cones of $m \times m$ positive semidefinite and positive definite Hermitian matrices. Also, let ${\mathcal S}_{\geq 0}(m,\bbR)$ and ${\mathcal S}_{> 0}(m,\bbR)$ be the their real-valued analogs. Let $O(m)$ and $U(m)$ be the orthogonal and unitary groups of $m \times m$ matrices, respectively. Throughout the manuscript, $\|M\|$ denotes the Euclidean norm of a matrix $M \in M(p,\bbR)$ in arbitrary dimension $p$, i.e.,
$\|M\| = \sqrt{\sup_{{\mathbf u} \in S^{p-1}} {\mathbf u}^* M M^* {\mathbf u}} = \sqrt{\sup_{{\mathbf u} \in S^{p-1}} {\mathbf u}^* M^* M {\mathbf u}}$. For any $M \in {\mathcal S}(p,\bbR)$,

\begin{equation}\label{e:lambda1(M)=<...=<lambdap(M)}
- \infty < \lambda_1(M) \leq \hdots \leq \lambda_q(M) \leq \hdots \leq \lambda_p(M) < \infty
\end{equation}

denotes the set of ordered eigenvalues of the matrix $M$. For $S = (s_{i_1,i_2})_{i_1,i_2=1,\dots,n}\in M(n,\mathbb{R})$, let

\begin{equation}\label{e:vec_definitions}
\textnormal{vec}_{{\mathcal S}}(S)=(s_{11},s_{21},\dots,s_{n1},s_{22},s_{32},\dots,s_{n2},\dots,s_{nn}).
\end{equation}

We use the asymptotic notation

\begin{equation}\label{e:OP(1),OmegaP(1)}
o_\bbP(1), \quad O_\bbP(1), \quad \Omega_\bbP(1),
\end{equation}

to describe sequences of random vectors or matrices whose norms vanish, are bounded above and bounded away from zero, respectively, in probability.

§.§ Wavelet analysis

Throughout the paper, we make use of a wavelet multiresolution analysis (MRA; see Mallat <cit.>, chapter 7), which decomposes $L^2(\mathbb{R})$ into a sequence of approximation (low-frequency) and detail (high-frequency) subspaces $V_j$ and $W_j$, respectively, associated with different scales of analysis $2^{j}$. In almost all mathematical statements, we assume the following conditions hold on the underlying wavelet MRA.

Assumption $(W1)$: $\psi \in L^2(\bbR)$ is a wavelet function, namely, it satisfies the relations

\begin{equation}\label{e:N_psi}
\int_{\bbR} \psi^2(t)dt = 1 , \quad \int_{\bbR} t^{p}\psi(t)dt = 0, \quad p = 0,1,\hdots, N_{\psi}-1, \quad \int_{\bbR} t^{N_{\psi}}\psi(t)dt \neq 0,
\end{equation}

for some integer (number of vanishing moments) $N_{\psi} \geq 2$.

Assumption ($W2$): the scaling and wavelet functions

\begin{equation}\label{e:supp_psi=compact}
\textnormal{$\phi\in L^1(\bbR)$ and $\psi\in L^1(\bbR)$ are compactly supported }
\end{equation}

and $\widehat{\phi}(0)=1$.

Assumption $(W3)$: there is $\alpha > 1$ such that

\begin{equation}\label{e:psihat_is_slower_than_a_power_function}
\sup_{x \in \bbR} |\widehat{\psi}(x)| (1 + |x|)^{\alpha} < \infty.
\end{equation}

Assumption ($W4$): the function

\begin{equation}\label{e:sum_k^m_phi(.-k)}
\sum_{k\in \mathbb{Z}}k^m\phi(\cdot-k)
\end{equation}

is a polynomial of degree $m$ for all $m=0,\hdots,N_{\psi}-1$.

Conditions (<ref>) and (<ref>) imply that $\widehat{\psi}(x)$ exists, is everywhere differentiable and its first $N_{\psi}-1$ derivatives are zero at $x = 0$. Condition (<ref>), in turn, implies that $\psi$ is continuous (see Mallat <cit.>, Theorem 6.1) and, hence, bounded.

Note that assumptions ($W1-W4$) are closely related to the broad wavelet framework for the analysis of $\kappa$-th order ($\kappa \in \bbN \cup \{0\}$) stationary increment stochastic processes laid out in Moulines et al. <cit.> and Roueff and Taqqu <cit.>. The Daubechies scaling and wavelet functions generally satisfy ($W1-W4$) (see Moulines et al. <cit.>, p. 1927, or Mallat <cit.>, p. 253). Usually, the parameter $\alpha$ increases to infinity as $N_{\psi}$ goes to infinity (see Moulines et al. <cit.>, p. 1927, or Cohen <cit.>, Theorem 2.10.1).

We further suppose the wavelet coefficients stem from Mallat's pyramidal algorithm (Mallat <cit.>, chapter 7). For expositional simplicity, in our description of the algorithm we use the $\bbR^p$-valued process $Y$, though analogous developments also hold for both $X$ and $Z$. Initially, suppose an infinite sequence of (generally dependent) random vectors

\begin{equation}\label{e:infinite_sample}
\{Y(\ell)\}_{\ell \in \bbZ},
\end{equation}

associated with the starting scale $2^j = 1$ (or octave $j = 0$), is available. Then, we can apply Mallat's algorithm to extract the so-named approximation $(A(2^{j+1},\cdot))$ and detail $(D(2^{j+1},\cdot))$ coefficients at coarser scales $2^{j+1}$ by means of an iterative procedure. In fact, as commonly done in the wavelet literature, we initialize the algorithm with the process

\begin{equation}\label{e:Btilde}
\bbR^p \ni \widetilde{Y}(t) := \sum_{k \in \bbZ} Y(k)\phi(t-k), \quad  t\in \bbR.
\end{equation}

By the orthogonality of the shifted scaling functions $\{\phi(\cdot - k)\}_{k \in \bbZ}$,

\begin{equation}\label{e:a(0,k)}
\bbR^p \ni  A(2^0,k)= \int_\bbR \widetilde{Y}(t)\phi(t-k)dt= Y(k), \quad k \in \bbZ
\end{equation}

(see Stoev et al. <cit.>, proof of Lemma 6.1, or Moulines et al. <cit.>, p. 160; cf. Abry and Flandrin <cit.>, p. 33). In other words, the initial sequence, at octave $j = 0$, of approximation coefficients is given by the original sequence of random vectors. To obtain approximation and detail coefficients at coarser scales, we use Mallat's iterative procedure

\begin{equation}\label{e:Mallat}
A(2^{j+1},k) = \sum_{k'\in \mathbb{Z}} u_{k'-2k}A(2^j,k'),\quad D(2^{j+1},k) =\sum_{k'\in \mathbb{Z}}v_{k'-2k} A(2^j,k'), \quad j \in \mathbb{N} \cup \{0\}, \quad k \in \mathbb{Z},
\end{equation}

where the (scalar) filter sequences $\{ u_k:=2^{-1/2}\int_\bbR \phi(t/2)\phi(t-k)dt \}_{k\in\bbZ}$, $\{v_k:=2^{-1/2}\int_\bbR\psi(t/2)\phi(t-k)dt\}_{k\in\bbZ}$ are called low- and high-pass MRA filters, respectively. Due to the assumed compactness of the supports of $\psi$ and of the associated scaling function $\phi$ (see condition (<ref>)), only a finite number of filter terms is nonzero, which is convenient for computational purposes (Daubechies <cit.>). Hereinafter, we assume without loss of generality that $\text{supp}(\phi) = \text{supp}(\psi)=[0,T]$ (cf. Moulines et al <cit.>, p. 160). Moreover, the wavelet (detail) coefficients $D(2^j,k)$ of $Y$ can be expressed as

\begin{equation}\label{e:disc2}
\bbR^p \ni D(2^j,k) = \sum_{\ell\in\bbZ} Y(\ell)h_{j,2^jk-\ell},
\end{equation}

where the filter terms are defined by

\begin{equation}\label{e:hj,l}
\bbR \ni h_{j,\ell} =2^{-j/2}\int_\bbR\phi(t+\ell)\psi(2^{-j}t)dt.
\end{equation}

If we replace (<ref>) with the realistic assumption that only a finite length series

\begin{equation}\label{e:finite_sample}
\{Y(k)\}_{k=1,\hdots,n}
\end{equation}

is available, writing $\widetilde{Y}^{(n)}(t) := \sum_{k=1}^{n}Y(k)\phi(t-k)$, we have $\widetilde{Y}^{(n)}(t) = \widetilde{Y}(t)$ for all $t\in (T,n+1)$ (cf. Moulines et al. <cit.>). Noting $D(2^j,k) = \int_\bbR \widetilde Y(t) 2^{-j/2}\psi(2^{-j}t - k)dt$ and $D^{(n)}(2^j,k) = \int_\bbR \widetilde Y^{(n)}(t)2^{-j/2}\psi(2^{-j}t - k)dt$, it follows that the finite-sample wavelet coefficients $D^{(n)}(2^j,k)$ of $\widetilde{Y}^{(n)}(t) $ are equal to $D(2^j,k)$ whenever $\textnormal{supp } \psi(2^{-j} \cdot - k) = (2^jk, 2^j (k+T)) \subseteq (T,n+1)$. In other words,

\begin{equation}\label{e:d^n=d}
D^{(n)}(2^j,k) = D(2^j,k),\qquad \forall (j,k)\in\{(j,k): 2^{-j}T\leq k\leq 2^{-j}(n+1)-T\}.
\end{equation}

Equivalently, such subset of finite-sample wavelet coefficients is not affected by the so-named border effect (cf. Craigmile et al. <cit.>, Percival and Walden <cit.>, Didier and Pipiras <cit.>). Moreover, by (<ref>) the number of such coefficients at octave $j$ is given by $n_j = \lfloor 2^{-j}(n+1-T)-T\rfloor$. Hence, $n_j \sim 2^{-j}n$ for large $n$. Thus, for notational simplicity we suppose

\begin{equation}\label{e:nj=n/2^j}
n_{j} = \frac{n}{2^j}
\end{equation}

holds exactly and only work with wavelet coefficients unaffected by the border effect.

§.§ Wavelet random matrices and assumptions

Throughout the paper, we assume observations stem from the model (<ref>). The independent “signal" $X = \{X(t)\}$ and the noise component $Z = \{Z(t)\}$ are $\bbR^r$-valued and $\bbR^p$-valued, respectively, where $r$ is fixed and $p = p(n)$.  The deterministic matrix $P = P(n)$ can be expressed as

\begin{equation}\label{e:P(n)}
M(p,r,\bbR) \ni P(n) = \Big( \mathbf{p}_{ 1}(n),\hdots,\mathbf{p}_{ r}(n) \Big), \quad \|\mathbf{p}_{ q}(n)\| = 1, \hspace{3mm} q = 1,\hdots,r.
\end{equation}

For the sake of clarity and mathematical generality, we state the appropriate conditions for the convergence in probability as well as for the asymptotic normality of wavelet log-eigenvalues directly in the wavelet domain. For $j \in \bbN$, $k \in \bbZ$ and a dyadic sequence $\{a(n)\}_{n \in \bbN}$, the random vectors

D(a(n)2^j,k)  \in \bbR^p, \quad  D_X(a(n)2^j,k)   \in \bbR^r \quad \textnormal{and} \quad D_Z(a(n)2^j,k) \in \bbR^p

denote the wavelet transform at scale $a(n)2^j$ of the stochastic processes $Y$, $X$ or $Z$, respectively. Whenever well-defined, the wavelet random matrix – or sample wavelet (co)variance – of $Y$ at scale $a(n)2^j$ is denoted by

\begin{equation}\label{e:W(a(nu))}
{\mathcal S}_{\geq 0}(p,\bbR) \ni \mathbf{W}(a(n) 2^j)) = \frac{1}{n_{a,j}}\sum^{n_{a,j}}_{k=1}D(a(n)2^j,k)D(a(n)2^j,k)^*, \quad n_{a,j} = \frac{n}{a(n)2^j}.
\end{equation}

The remaining wavelet random matrix terms $\W_X,\W_{XZ},\W_Z$ are naturally defined as

\mathbf{W}_X(a(n) 2^j)) = \frac{1}{n_{a,j}}\sum^{n_{a,j}}_{k=1}D_X(a(n)2^j,k)D_X(a(n)2^j,k)^* \in {\mathcal S}_{\geq 0}(r,\bbR),

\mathbf{W}_{XZ}(a(n) 2^j)) = \frac{1}{n_{a,j}}\sum^{n_{a,j}}_{k=1}D_X(a(n)2^j,k)D_Z(a(n)2^j,k)^* \in M(r,p,\bbR),

\begin{equation}\label{e:WZ(a(n)2^j)}
\mathbf{W}_{Z}(a(n) 2^j)) = \frac{1}{n_{a,j}}\sum^{n_{a,j}}_{k=1}D_Z(a(n)2^j,k)D_Z(a(n)2^j,k)^* \in {\mathcal S}_{\geq 0}(p,\bbR).
\end{equation}

Further define the auxiliary random matrix

\begin{equation}\label{e:B-hat_a(2^j)}
\widehat{{\mathbf B}}_a(2^j) = P_H^{-1}a(n)^{-H-(1/2)I}{\mathbf W}_X(a(n)2^j)a(n)^{-H^*-(1/2)I} (P_H^*)^{-1}\in {\mathcal S}_{\geq 0}(r,\bbR),
\end{equation}

as well as its mean ${\mathbf B}_a(2^j) := \bbE \widehat{{\mathbf B}}_a(2^j)$. The matrix $\widehat{{\mathbf B}}_a(2^j)$ should be interpreted as a version of ${\mathbf W}_X(a(n)2^j)$ after compensating for scaling and non-canonical coordinates (cf. relation (<ref>), which displays canonical scaling). In (<ref>), we assume that the scaling matrix $H$ has the Jordan form

\begin{equation}\label{e:H_is_diagonalizable}
H = P_H \textnormal{diag}(h_1,\hdots,h_r) P^{-1}_H, \quad  P_H \in GL(r,\bbR), \quad -1/2 < h_1 \leq \hdots \leq h_r <\infty.
\end{equation}

For the sake of illustration, when $X$ is an ofBm, $H$ is a Hurst matrix whose ordered eigenvalues satisfy $0 < h_1 \leq \hdots \leq h_r \leq 1$ (see Example <ref>).

We make use of the following assumptions in the main results of this paper (Section <ref>). For expository purposes, we first state the assumptions, and then provide some interpretation.

Assumption $(A1)$: The wavelet random matrix

{\mathbf W}(a(n)2^j) = P(n){\mathbf W}_X(a(n)2^j)P^*(n) + {\mathbf W}_Z(a(n)2^j)

\begin{equation}\label{e:W=PWXP*+WZ+PWXZ+WXZ*P*}
+ P(n){\mathbf W}_{X,Z}(a(n)2^j)+ {\mathbf W}^*_{X,Z}(a(n)2^j)P^*(n)
\end{equation}

is well defined a.s.

Assumption $(A2)$: In (<ref>), for $H$ as in (<ref>),

\begin{equation}\label{e:assumptions_WZ=OP(1)}
\max \Big\{\|\E \mathbf{W}_Z(a(n)2^j)\|, \|\W_Z(a(n)2^j)\| , \|a(n)^{-H-(1/2)I} \mathbf{W}_{X,Z}(a(n) 2^j))\|\Big\} = O_{\bbP}(1).
\end{equation}

Assumption $(A3)$: The random matrix $\widehat{{\mathbf B}}_a(2^j)$ as in (<ref>) satisfies

\begin{equation}\label{e:sqrt(K)(B^-B)->N(0,sigma^2)}
\Big(\hspace{1mm}\sqrt{n_{a,j}}\hspace{0.5mm}(\vecoper_{{\mathcal S}}\widehat{{\mathbf B}}_a(2^j) - \vecoper_{{\mathcal S}}{\mathbf B}_a(2^j))\hspace{1mm}\Big)_{j=j_1,\hdots,j_m}\stackrel{d}\rightarrow {\mathcal N}(0,\Sigma_B(j_1,\hdots,j_m)), \quad n \rightarrow \infty.
\end{equation}

In addition,

\begin{equation}\label{e:|bfB_a(2^j)-B(2^j)|=O(shrinking)}
\|{{\mathbf B}}_a(2^j) - \mathbf B(2^j) \| = O(a(n)^{-\varpi}), \quad  n\to\infty, \quad j=j_1,\hdots,j_m,
\end{equation}

where $\varpi$ is the regularity parameter

\begin{equation}\label{e:varpi_parameter}
\varpi = \min \Big\{\hspace{0.5mm}h_1,\min_{1 \leq q_1 < q_2 \leq r}(h_{q_2}-h_{q_1}) \Big\}1_{\{\min_{1 \leq q_1 < q_2 \leq r}( h_{q_2}- h_{q_1}) >0\}}
+  h_1 1_{\{h_{q_1}=\ldots =  h_{q_r} \}},
\end{equation}

and also

\begin{equation}\label{e:B(2^j)_full_rank}
\mathbf B(2^j) \in \mathcal S_{>0}(r,\bbR) \quad \textnormal{and} \quad \det(\Sigma_B(j_1,\hdots,j_m)> 0.
\end{equation}

Assumption $(A4)$: The dimension $p(n)$ and the scaling factor $a(n)$ satisfy the relations

\begin{equation}\label{e:p(n),a(n)_conditions}
a(n) \leq \frac{n}{2^j}, \quad \frac{a(n)}{n}+ \frac{n}{a(n)^{2 \varpi +1}} \rightarrow 0, \quad \frac{p(n)}{n/a(n)} \rightarrow c \in [0,\infty), \quad n \rightarrow \infty,
\end{equation}

where $\varpi$ is given as in (<ref>).

Assumption $(A5)$: Let $P(n) \in M(p,r,\bbR)$ and $P_H \in GL(r,\bbR)$ be as in (<ref>) and (<ref>), respectively. Let $P(n)P_H=Q(n)R(n)$ be the $QR$ decomposition of $P(n)P_H$ and let $\varpi$ be as in (<ref>). Then, there exists a (deterministic) matrix $A \in \mathcal{S}_{>0}(r,\bbR)$ with Cholesky decomposition
$A = R^* R$
such that

\begin{equation}\label{e:<p1,p2>=c12_2}
\|R(n)-R\| = O(a(n)^{-\varpi}).
\end{equation}

Assumptions $(A1 - A3)$ are stated in the wavelet domain. Assumption $(A1)$ holds under very general conditions. In fact, since $X$ and $Z$ are assumed independent, it suffices that $\|{\mathbf W}_{X}(a(n)2^j)\| < \infty$ and $\|{\mathbf W}_{Z}(a(n)2^j)\| < \infty$ a.s., which holds as long as $X$ and $Z$ are well-defined discrete time stochastic processes. Assumption $(A2)$ ensures that the  influence of the random matrices ${\mathbf W}_{X,Z}(a(n)2^j)$ and ${\mathbf W}_{Z}(a(n)2^j)$ in the observed wavelet spectrum ${\mathbf W}(a(n)2^j)$ is not too large. Assumption $(A3)$ posits the asymptotic normality of the (wavelet domain) fractional component ${\mathbf W}_{X}(a(n)2^j)$ after compensating for scaling and non-canonical coordinates.

In turn, Assumption $(A4)$ controls the divergence rates among $n$, $a(n)$ and $p(n)$. In particular, it states that the scaling factor $a(n)$ must blow up slower than $n$, and that the three-component ratio $\frac{p(n) a(n)}{n}$ must converge to a constant (cf. the traditional ratio $\lim_{n \rightarrow \infty} p(n)/n$ for high-dimensional sample covariance matrices). Assumption $(A5)$ ensures that, asymptotically speaking, the angles between the column vectors of the matrix $P(n) P_H$ converge in such a way that the matrix of asymptotic angles $\lim_{n \rightarrow \infty}P^*_H P^*(n)P(n)P_H = A$ has full rank. This entails that $P(n)$ does not strongly perturb the scaling properties of the hidden random matrix ${\mathbf W}_{X}(a(n)2^j)$.

A discussion of broad contexts where assumptions $(A2)$ and $(A3)$ hold is deferred to Section <ref>. In the following example, to fix ideas we briefly illustrate assumption $(A3)$.

Recall that an ofBm is a multivariate $(i)$ Gaussian; $(ii)$ o.s.s.; $(iii)$ stationary-increment stochastic process. Suppose the $\bbR^r$-valued stochastic process $X$ is an ofBm, and rewrite $H = P_H J_H P^{-1}_H \in M(r,\bbR)$, where $P_H \in GL(r,\bbC)$ and $J_H $ is the Jordan form of $H$. It can be shown that most instances of ofBm admit the spectral density

\begin{equation}\label{e:ofBm_harmonizable_repres}
{\mathfrak{g}}_X(x) = x^{-(H+(1/2)I)}_+ \hspace{0.25mm} {\mathcal A}{\mathcal A}^* \hspace{0.25mm}x^{-(H^*+(1/2)I)}_+ + x^{-(H+(1/2)I)}_- \hspace{0.25mm}\overline{{\mathcal A}{\mathcal A}^*} \hspace{0.25mm}x^{-(H^*+(1/2)I)}_-, \quad x \in \bbR \backslash\{0\},
\end{equation}

for some ${\mathcal A} \in M(r,\bbC)$ (Didier and Pipiras <cit.>, Theorem 3.1).
So, fix $j \in \bbN$. Under mild assumptions on the ofBm, the associated auxiliary random matrix $\widehat{{\mathbf B}}_a(2^j)$ as in (<ref>) satisfies Assumption ($A3$) as a consequence of Theorem 3.1 and Lemma C.2 (extended to dimension $r$) in Abry and Didier <cit.>. In other words, relations (<ref>) and (<ref>) hold for a matrix sequence ${\mathbf B}_a(2^j) \in {\mathcal S}_{\geq0}(\bbR,r)$ and a matrix ${\mathbf B}(2^j) \in {\mathcal S}_{>0}(\bbR,r)$. Moreover,
if $H$ is diagonalizable with real eigenvalues, it can be easily shown that the matrix $\B(2^j)$ satisfies the entrywise (i.e., along canonical axes) scaling relations

\begin{equation}\label{e:B(2^j)_entrywise_scaling}
{\mathbf B}(2^j) = \Big( {\mathbf b}(2^j)_{ii'} \Big)_{i,i'=1,\hdots,r} = \Big( 2^{j(h_i + h_{i'}+1)}{\mathbf b}(1)_{ii'} \Big)_{i,i'=1,\hdots,r}.
\end{equation}

§ MAIN RESULTS

§.§ Asymptotic behavior of wavelet eigenvalues

We first establish that, after proper rescaling, the $r$ largest eigenvalues of a wavelet random matrix $\mathbf{W}(a(n)2^j)$ in high dimensions converge in probability to deterministic functions $\xi_q(2^j)$, $q = 1,\hdots,r$. Thus, such functions can be interpreted as asymptotic rescaled eigenvalues. Notably, they display a scaling property. Moreover, the remaining $p(n)-r$ eigenvalues of a wavelet random matrix are bounded in probability.

Assume $(W1 - W4)$ and $(A1 - A5)$ hold and fix $j \in \bbN$.
Then, for $p = p(n)$, the limits

\begin{equation}\label{e:lim_n_a*lambda/a^(2h+1)}
\plim_{n \rightarrow \infty}\frac{\lambda_{p-r+q}\big(\mathbf{W}(a(n)2^j)\big)}{a(n)^{2h_q+ 1}} =: \xi_q(2^j), \quad q=1,\ldots,r,
\end{equation}

exist, and the deterministic functions $\xi_q$ satisfy the scaling relation

\begin{equation}\label{e:xi-i0_scales}
\xi_{q}(2^j)  = 2^{j \hspace{0.5mm}(2 h_{q}+1)} \xi_{q}(1).
\end{equation}

In addition,

\begin{equation}\label{e:lim_lambda_p-r(W)}
0\leq \lambda_{1}\big(\mathbf{W}(a(n)2^j)\big) \leq \ldots \leq \lambda_{p-r}\big(\mathbf{W}(a(n)2^j)\big) = O_{\bbP}(1).
\end{equation}

In the following theorem, we establish the asymptotic normality of the $r$ largest wavelet log-eigenvalues in high dimensions. The statement covers the central cases of simple ($h_1 < \hdots < h_r$) and of identical ($h_1 = \hdots = h_r$) scaling eigenvalues. Note that, by comparison to Theorem <ref>, the asymptotic normality of wavelet log-eigenvalues requires some additional condition so as to ensure the finite-sample differentiability of the eigenvalues with respect to scaling parameters (see also Remark <ref>).

Suppose $(W1 - W4)$ and $(A1 - A5)$ hold, and fix $j_1,j_2 \in \bbN$, $j_1<j_2$. Further suppose one of the following conditions holds, namely,

(i) either

\begin{equation}\label{e:assumption:simple_Hurst}
-1/2 < h_1 < \ldots < h_r < \infty;
\end{equation}

(ii) or $h_1 = \ldots = h_r$ and the functions $\xi_q(2^j)$ in (<ref>) satisfy

\begin{equation}\label{e:xiq_distinct}
q_1\neq q_2 \Rightarrow \xi_{q_1}(1)\neq \xi_{q_2}(1).
\end{equation}

Then, for $p = p(n)$,

\begin{equation}\label{e:asympt_normality_lambda2}
\bbR^{m \times r} \ni \Big( \sqrt{n_{a,j}}\Big( \log \lambda_{p-r+q}\big({\mathbf W}(a(n)2^{j})\big)  - \log \lambda_{p-r+q}\big(\bbE {\mathbf W} (a(n)2^{j})\big) \Big)_{q=1,\hdots,r} \Big)_{j=j_1,\hdots,j_2} \stackrel{d}\rightarrow {\mathcal N}(0,\Sigma_{\lambda})
\end{equation}

as $n \rightarrow \infty$. If we write the asymptotic covariance matrix in block form $\Sigma_{\lambda} = \Big( \Sigma_{\lambda}(jj') \Big)_{j,j'= j_1,\hdots,j_2}$, then its main diagonal entries satisfy $\Sigma_{\lambda}(jj)_{ii}> 0$, $i = 1,\hdots,r$.

Under non-simple scaling matrix eigenvalues and without condition (<ref>), the asymptotic distribution of wavelet log-eigenvalues is generally non-Gaussian (in low dimensions, see Proposition A.1.1 in Li <cit.>; cf. Anderson <cit.>, Sections 13.5.1 and 13.5.2, on sample covariance matrices).

§.§ Consequences for statistical inference: a short discussion

Theorems <ref> and <ref> bear direct consequences for statistical inference. For example, fix $j$. In light of (<ref>), (<ref>) and (<ref>), the random vector

\Big\{\frac{\log \lambda_{p-r+q}\big({\mathbf W}(a(n)2^j)\big)}{2 \log a(n)}\Big\}_{q = 1,\hdots,r}

can be interpreted, in the language of Statistics, as consistent and asymptotically normal estimators of the vector of parameters $\{h_{q}\}_{q=1,\hdots,r}$ in high dimensions. Under the same conditions, the lowest $p-r$ wavelet log-eigenvalues stay bounded, whence

\begin{equation}\label{e:log_lambda_p-r/2log_a}
\frac{\log \lambda_{p-r}\big({\mathbf W}(a(n)2^j )\big)}{2 \log a(n)}
\end{equation}

converges to zero in probability (see Figure <ref>). In other words, Theorems <ref> and <ref> provide a framework for the high-dimensional estimation of the parameters $h_1,\hdots,h_r$ and $r$, and hence, of the scaling properties of the system (<ref>). Moreover, Theorem <ref> under condition $(i)$ can be used in testing the hypothesis of the equality of scaling eigenvalues.

§ EXAMPLES

Recall that assumptions $(A2)$ and $(A3)$ (i.e., (<ref>) and (<ref>)-(<ref>), respectively) are stated in the wavelet domain. In this section, we provide a Gaussian framework as well as Gaussian and non-Gaussian examples where such assumptions are satisfied.

§.§ A Gaussian framework

When $X$ and $Z$ are Gaussian, it can be shown that conditions (<ref>)–(<ref>) hold under very general assumptions. For the sake of generality, we consider an array of (colored) noise density functions. In other words, we assume that, for each $n$, a new set of such functions is randomly picked (see also Example <ref>). Also, to facilitate comparison with previous results in low dimensions (Abry and Didier <cit.>), we suppose the hidden process $X$ is an ofBm (see also Remark <ref>). Recall that, when $X$ is an ofBm, the auxiliary random matrix $\widehat{{\mathbf B}}_a(2^j)$ satisfies (<ref>) and (<ref>) (i.e., Assumption $(A3)$; see Example <ref>).

Suppose assumptions $(W1-W4)$ are in place. Further suppose the following conditions hold.

$(i)$ The $\bbR^p$-valued (noise) stochastic process $Z = \{Z(t)\}_{t \in \bbZ}$ as in (<ref>) has spectral density

\begin{equation}\label{e:Z(t)}
{\mathfrak g}_Z(x) \in {\mathcal S}_{\geq 0}(p,\bbR), \quad x \in [-\pi,\pi).
\end{equation}

In (<ref>), ${\mathfrak g}_Z = \textnormal{diag}({\mathfrak g}_{1,n}, \hdots, {\mathfrak g}_{p,n})$, where ${\mathfrak g}_{\ell,n}$, $\ell = 1,\hdots,p$, are drawn independently according to a distribution $F(d{\mathfrak g})$ on a class

\begin{equation}\label{e:G-class_in L2}
{\mathcal G} \subseteq L^1[-\pi,\pi)
\end{equation}

of univariate functions ${\mathfrak g}$ such that

\begin{equation}\label{e:sup_g_in_G|g|<infty}
\sup_{{\mathfrak g} \in {\mathcal G}} \hspace{0.5mm}\sup_{x \in [-\pi,\pi)} |{\mathfrak g}(x)| < \infty, \quad |\E {\mathfrak g}(x) - \E {\mathfrak g}(0) | = O(x^2), \quad x\to 0.
\end{equation}

In addition, suppose that, conditionally on $\mathfrak{g}$, the univariate covariance function $\gamma_{Z_1,{\mathfrak g}}$ satisfies

\begin{equation}\label{e:|gammaZ,g|=<exp}
|\gamma_{Z_1,{\mathfrak g}}(h)| \leq C_1 e^{- C_2 h}, \quad h \in \bbZ,
\end{equation}

where the constants $C_1,C_2 > 0$ do not depend on the (univariate) density $\mathfrak{g}$;

$(ii)$ the measure $F(d\mathfrak{g})$ is independent of $X$ in (<ref>);

$(iii)$ $0<h_1\leq\ldots\leq h_r<1$, and the stationary-increment, $\bbR^r$-valued stochastic process $X = \{X(t)\}_{t \in \bbZ}$ is an ofBm with Hurst matrix $H = P_H \textnormal{diag}(h_1,\hdots,h_n)P^{-1}_H$, $P_H \in GL(r,\bbR)$, whose spectral density (<ref>) is such that $\det\Re(\mathcal{A}\mathcal{A}^*) >0$.

Then, condition (<ref>) (i.e., Assumption $(A2)$) is satisfied.

Suppose that, for any fixed $n$, the $r$-variate ofBm $X$ and the $p(n)$-variate process $Z$ are independent. Further assume that, for fixed integers ${\mathfrak p}$ and ${\mathfrak q}$, $Z$ is made up of independent copies of an ARMA$({\mathfrak p},{\mathfrak q})$ processes (in particular, $F(d{\mathfrak g})$ is a Dirac delta on the associated spectral density). Then, the assumptions of Proposition <ref> are met.

The assumption in Proposition <ref> that the process $X$ is an ofBm or even that it has (first order) stationary increments is not crucial. It is well known that wavelet frameworks are suitable for stationary increment processes of any order. This is so because, for such processes, wavelet coefficients are, in general, stationary as long as the chosen number of vanishing moments of the underlying wavelet basis (see (<ref>)) is sufficiently large. Such topic has been broadly explored in the literature (e.g., Flandrin <cit.>, Wornell and Oppenheim <cit.>, Veitch and Abry <cit.>, Moulines et al. <cit.>, Roueff and Taqqu <cit.>; see Abry et al. <cit.> on the multivariate stochastic processes). In particular, condition (<ref>) on scaling eigenvalues allows for stationary fractional processes $X$ exhibiting long-range dependence (see, for instance, Embrechts and Maejima <cit.>, Pipiras and Taqqu <cit.>).

§.§ Non-Gaussian instances: a short discussion

Chracterizing the wavelet domain behavior of high-dimensional, non-Gaussian, second-order and fractional frameworks is a very broad topic that lies well outside the scope of this paper. In this section, we restrict ourselves to discussing certain non-Gaussian instances so as to help illustrate the fact that the wavelet domain properties established in Section <ref> do not fundamentally require the system (<ref>) to be Gaussian.

Let $X = \{X(t)\}_{t \in \bbZ}$ be a possibly non-Gaussian, $\bbR^r$-valued stochastic processes made up of independent linear fractional processes with finite second moments. Then, based on the framework constructed in Roueff and Taqqu <cit.>, we can show that the associated random matrix $\widehat{{\mathbf B}}_a(2^j) \in {\mathcal S}_{>0}(r,\bbR)$ satisfies conditions (<ref>) and (<ref>) (i.e., Assumption $(A3)$) under conditions $(W1-W4)$ and mild additional assumptions on $X$. This claim is made precise in Proposition <ref>.

Even though – in this particular case – the $r$ components of $X$ are assumed independent, note that statistically estimating the $r$ univariate scaling exponents based on the associated model (<ref>) is still, in general, a nontrivial problem. This is so due to the presence of the unknown coordinates matrix $P(n)P_H = P(n)$, as well as of the high-dimensional noise component $Z$. In fact, mathematically speaking, such model is closely related to the so-called blind source separation problems, which are of great interest in the field of signal processing (e.g., Naik and Wang <cit.>; in a fractional context, see Abry et al. <cit.>).

Recall that a distribution is called sub-Gaussian when its tails are no heavier than those of the Gaussian distribution (Vershynin <cit.>, Proposition 2.5.2). Sub-Gaussian distributions form a broad family that includes the Gaussian distribution itself, as well as compactly supported distributions, for example. Suppose the noise process $\{Z(t)\}_{t \in \bbZ}$ consists of (discrete-time) i.i.d. sub-Gaussian observations. Consider the Haar wavelet framework, where the wavelet coefficients are computed by means of Mallat's iterative procedure (<ref>). Then, for any $j \in \bbN$, the nonzero coefficients in the sequences $\{h_{j,2^{j}k - \ell}\}_{\ell \in \bbZ}$ and $\{h_{j,2^{j}k' - \ell}\}_{\ell \in \bbZ}$ do not overlap for $k \neq k'$ (cf. (<ref>)). For this reason, the wavelet coefficients are independent at any fixed scale. Although we do not provide a proof due to space constraints, it is then possible to use traditional concentration of measure techniques to show that the wavelet random matrix ${\mathbf W}_{Z}$ satisfies condition (<ref>) in Assumption $(A2)$. The study of the properties of wavelet random matrices ${\mathbf W}_{Z}$ under other wavelet bases or other classes of non-Gaussian observations is currently a broad topic of research.

In general, depending on the properties of the latent process $X$, the random matrix $\widehat{{\mathbf B}}_a(2^j)$ may not be asymptotically Gaussian, i.e., assumption $(A3)$ may not hold. In the univariate context ($r=1$), see, for instance, Bardet and Tudor <cit.>, Clausel et al. <cit.>. Nevertheless, for $r > 1$, the literature still lacks broad characterizations of conditions under which $\widehat{{\mathbf B}}_a(2^j)$ (or ${\mathbf W}_X(a(n)2^j)$) is asymptotically non-Gaussian.

§ CONCLUSION AND OPEN PROBLEMS

In this paper, we mathematically characterize the asymptotic and large scale behavior of the eigenvalues of wavelet random matrices in high dimensions. We assume that possibly non-Gaussian, finite-variance $p$-variate measurements are made of a low-dimensional $r$-variate ($r \ll p$) fractional stochastic process with unknown scaling coordinates and in the presence of additive high-dimensional noise. In the high-dimensional and large scale limit, we establish that the rescaled $r$ largest eigenvalues of the wavelet random matrices converge to scale invariant functions, whereas the remaining $p-r$ eigenvalues remain bounded. In addition, under slightly stronger assumptions, we show that the $r$ largest log-eigenvalues of wavelet random matrices exhibit asymptotically Gaussian distributions. The results bear direct consequences for high-dimensional modeling starting from measurements of the form of a signal-plus-noise system (<ref>), where $X$ is a latent process containing fractal information and $P$ (as well as $Z$) is unknown.

This research leads to many relevant open problems in high dimensions, some of which can be briefly described as follows. $(i)$ The results in Section <ref> build upon broad wavelet domain assumptions, hence providing a rich framework for future research pursuits involving wavelet random matrices. A natural direction of inquiry is the mathematical study of the properties of wavelet random matrices in second order non-Gaussian fractional instances, e.g., in regard to the mathematical control of the impact of heavier tails under the broad assumptions of Theorems <ref> and <ref>. $(ii)$ In turn, an interesting direction of extension is the characterization of sets of conditions under which the $r$ largest eigenvalues of wavelet random matrices exhibit non-Gaussian fluctuations (cf. Remarks <ref> and <ref>). $(iii)$ In modeling, starting from measurements of the general form (<ref>), the construction of an extended framework for the detection of scaling laws in high-dimensional systems calls for the investigation of the behavior of wavelet random matrices when $r \rightarrow \infty$. In particular, such undertaking requires a deeper study, in the wavelet domain, of the so-named eigenvalue repulsion effect (e.g., Tao <cit.>), which may severely skew the observed scaling laws (see Wendt et al. <cit.> on preliminary computational studies). $(iv)$ Superior finite-sample statistical performance can be attained by means of a wavelet eigenvalue regression procedure across scales (cf. Abry and Didier <cit.>). Namely, fix a range of scales
$j = j_1 , j_1 + 1, \hdots , j_2$ and let $m = j_2 - j_1 + 1$.

\begin{equation}\label{e:hl-hat}
\{ \widehat{\ell}_q \}_{q = 1,\hdots,r} :=  \Big\{\frac{1}{2}\Big(\sum_{j=j_1}^{j_2} w_j \log_2 \lambda_{p-r+q}\big({\mathbf W}(a(n)2^j)\big) -1\Big)\Big\}_{q = 1,\hdots,r}.
\end{equation}

In (<ref>), $w_j$, $j = j_1,\hdots,j_2$, are weights satisfying the relations $\sum^{j_2}_{j=j_1}w_j = 0$, $\sum^{j_2}_{j=j_1}j w_j = 1$,
where $w_j = 1$ if $m = 1$. As an immediate consequence of Theorem <ref>, $\{ \widehat{\ell}_q \}_{q = 1,\hdots,r}$ provides a consistent estimator of $\{h_q\}_{q = 1,\hdots,r}$. Moreover, as in low dimensions (Abry and Didier <cit.>), Theorem <ref> points to asymptotic normality. However, additional results would be required, which remains a topic for future work. $(v)$ In applications, the results in this paper naturally pave the way for the investigation of scaling behavior in high-dimensional (“Big") data from fields such as physics, neuroscience and signal processing.

§ PROOFS: SECTION <REF>

In this section, whenever convenient we simply write $p = p(n)$, $a = a(n)$ or $P = P(n)$. For the proofs, it is useful to first recall the variational characterization of the eigenvalues of a matrix $M \in {\mathcal S}(p,\bbC)$ provided by the Courant–Fischer principle. In other words, fix $q \in \{1,\hdots,p\}$ and consider ordered eigenvalues as in (<ref>). By the Courant–Fischer principle, we can express

\begin{equation}\label{e:Courant--Fischer}
\lambda_{q}(M) = \inf_{{\mathcal U}_q} \sup_{{\mathbf u} \in {\mathcal U}_q \cap S^{p-1}_{\bbC}}{\mathbf u}^* M {\mathbf u} = \sup_{{\mathcal U}_{p-q+1}}\inf_{{\mathbf u} \in {\mathcal U}_{p-q+1} \cap S^{p-1}_{\bbC}} {\mathbf u}^*M{\mathbf u},
\end{equation}

where ${\mathcal U}_q$ is a $q$-dimensional subspace of $\bbC^n$ (Horn and Johnson <cit.>, Chapter 4). A related useful idea is the following. For $M \in {\mathcal S}(p,\bbR)$, let ${\mathbf u}_1,\hdots, {\mathbf u}_p$ be unit eigenvectors of $M$ associated with the eigenvalues (<ref>) respectively. Then, for $q \in \{1,\hdots,p\}$, we can further express

\begin{equation}\label{e:lambdaq(M)_based_on_eigenvecs}
\lambda_{q}(M) = \inf_{{\mathbf u} \in \textnormal{span}\{{\mathbf u}_q,\hdots, {\mathbf u}_p\} \cap S^{p-1}}{\mathbf u}^* M {\mathbf u} = \sup_{{\mathbf u} \in \textnormal{span}\{{\mathbf u}_1,\hdots, {\mathbf u}_q\} \cap S^{p-1}}{\mathbf u}^* M {\mathbf u}.
\end{equation}

In proofs, we often use the following notation. For an arbitrary fixed $q\in \{1,\ldots,r\}$, $\mathcal{I}_-$, $\mathcal{I}_0$ and $\mathcal{I}_+$ are index sets given by the relations

\begin{equation} \label{e:def_indexsets}
\begin{array}{lll}
\mathcal{I}_- := \{\ell: h_\ell < h_q\}, \quad \mathcal{I}_0 := \{\ell: h_\ell =h_q\}, \quad \mathcal{I}_+ := \{\ell: h_\ell > h_q\}.
\end{array}
\end{equation}

Note that $\mathcal{I}_-$ and $\mathcal{I}_+$ are possibly empty. Also write

\begin{equation}\label{e:r1,r2,r3}
r_1:= \textnormal{card}(\mathcal{I}_-), \quad r_2:= \textnormal{card}(\mathcal{I}_0) \geq 1, \quad r_3:= \textnormal{card}(\mathcal{I}_+).
\end{equation}

Throughout this section, for notational simplicity, we write

\begin{equation}\label{e:P(n)PH_equiv_P}
P(n)P_H \equiv P(n) \equiv P,
\end{equation}

whose column vectors are denoted by ${\mathbf p}_{\ell}(n) \equiv {\mathbf p}_{\ell}$, $\ell = 1,\hdots,r$. Moreover, we use $\mathbf{h}$ as a shorthand for the matrix $\textnormal{diag}(h_1,\hdots,h_r)$.

Proof of Theorem <ref>:  Statement (<ref>) is a consequence of Lemma <ref>. The scaling relationship (<ref>) is an immediate consequence of the limit (<ref>), since if such a limit holds, writing $\widetilde a(n) = a(n)2^j$,

\xi_q(2^j)=\plim_{n \rightarrow \infty}\frac{\lambda_{p-r+q}\big(\mathbf{W}(a(n)2^j)\big)}{a(n)^{2h_q+ 1}} =\plim_{n \rightarrow \infty}\frac{\lambda_{p-r+q}\big(\mathbf{W}(\widetilde a(n))\big)}{(\widetilde a(n)2^{-j})^{2h_q+ 1}} = 2^{j(2h_q+1)}\xi_q(1).

So, we now show (<ref>).

Fix $q\in\{1,\ldots,r\}$. Let $\mathcal{I}_0$ and the (possibly empty) sets $\mathcal{I}_-,\mathcal{I}_+$  be given by the relations (<ref>), and consider $r_1$, $r_2$ and $r_3$ as defined in (<ref>). Let $P_H$ be as in (<ref>). Starting from conditions (<ref>) and (<ref>), consider the decomposition

\begin{equation}\label{e:rescaled_W}
\frac{{\mathbf W}(a2^j)}{a^{2h_q+1}} = \frac{Pa^{{\mathbf h}} \widehat{{\mathbf B}}_a(2^j)a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{O_{\bbP}(1)}{a^{2h_q+1}}+
\frac{Pa^{{\mathbf h}}O_{\bbP}(1)}{a^{2h_q+1/2}}+\frac{O_{\bbP}(1)a^{{\mathbf h}}P^*}{a^{2h_q+1/2}}.
\end{equation}

If $\mathcal{I}_+ = \emptyset$, then $\frac{\lambda_{p-r+q}({\mathbf W}(a2^j))}{a^{2h_q+1}}=O_{\bbP}(1)$. Thus, it suffices to show that every subsequence of $\frac{\lambda_{p-r+q}({\mathbf W}(a2^j))}{a^{2h_q+1}}$ that is convergent in probability converges to the same limit.  So, suppose two such limits exist along some subsequences $n_{\mathfrak a}$, $n_{\mathfrak b}$, and write

\begin{equation}\label{e:2lim_xi}
\plim_{n_i \rightarrow \infty} \frac{\lambda_{p-r+q}( {\mathbf W}(a(n_i)2^j))}{a(n_i)^{2h_q+1}} =: \xi_{q,i}(2^j), \quad i = {\mathfrak a}, \mathfrak b.
\end{equation}

The goal is to show that

\begin{equation}\label{e:xi_q,a(2^j)=<xi_a,b(2^j)}
\xi_{q,{\mathfrak a}}(2^j)\leq\xi_{q,{\mathfrak b}}(2^j).
\end{equation}

Once (<ref>) is established, we can exchange the roles of $n_{\mathfrak a},n_{\mathfrak b}$ in the proof and obtain the desired equality

\begin{equation}\label{e:xi_q,a(2^j)=xi_q,b(2^j)}
\xi_{q,{\mathfrak a}}(2^j)=  \xi_{q,{\mathfrak b}}(2^j).
\end{equation}

So, by passing to further subsequences if necessary, for $i=\mathfrak a, \mathfrak b$ and $\ell=1,\ldots,r$ we can assume that there are deterministic vectors ${\boldsymbol \gamma}_{\ell}^i$ such that

\begin{equation}\label{e:P^*(ni)u_p-r+ell(ni)->gamma^i_ell}
P^*(n_i) {\mathbf u}_{p-r+\ell}(n_i)\to {\boldsymbol \gamma}_{\ell}^i\in \bbR^r \quad \textnormal{a.s.}, \quad n_i \rightarrow \infty.
\end{equation}

\begin{equation}\label{e:Gamma-i_def}
\Gamma_i:=({\boldsymbol \gamma}_{1}^i,\ldots, {\boldsymbol \gamma}_{r}^i)\in GL(r,\bbR) \quad \textnormal{a.s.}, \quad i=\mathfrak a,\mathfrak b,
\end{equation}

by Lemma <ref>. Now, take any nonzero vector

\begin{equation}\label{e:w_in_the_middle}
\widetilde{ \boldsymbol{\gamma}} \in\textnormal{span}\{{\boldsymbol \gamma}_{q}^{\mathfrak a},\ldots,{\boldsymbol \gamma}_{r}^{\mathfrak a}\} \cap \textnormal{span}\{{\boldsymbol \gamma}_{1}^{\mathfrak b},\ldots,{\boldsymbol \gamma}_{q}^{\mathfrak b}\}%
\end{equation}

whose (almost sure) existence is guaranteed since $\dim \textnormal{span}\{{\boldsymbol \gamma}_{q}^{\mathfrak a},\ldots,{\boldsymbol \gamma}_{r}^{\mathfrak a}\}+ \dim\textnormal{span}\{{\boldsymbol \gamma}_{1}^{\mathfrak b},\ldots,{\boldsymbol \gamma}_{q}^{\mathfrak b}\}= r+1$.  In view of (<ref>), there are scalars $\alpha_{q}, \ldots, \alpha_{r}$ such that

\begin{equation}\label{e:gammatilde_a-basis}
\widetilde{ \boldsymbol{\gamma}}= \sum^{r}_{\ell = q}\alpha_{\ell}\hspace{0.5mm}{\boldsymbol \gamma}_{\ell}^{\mathfrak a} =: \Gamma_{{\mathfrak a}}\hspace{0.5mm}{\boldsymbol \alpha} \ \in \bbR^{r}.
\end{equation}

Along the subsequence $n_a$, define the unit (random) vector

{\mathbf v}_{\mathfrak a}(n_{\mathfrak a}) = \frac{\sum^{r}_{\ell=q}\alpha_{\ell}{\mathbf u}_{p-r+\ell}(n_{\mathfrak a})}{ \|\sum^{r}_{\ell=q}\alpha_{\ell}{\mathbf u}_{p-r+\ell}(n_{\mathfrak a})\|} = \frac{\sum^{r}_{\ell=q}\alpha_{\ell}{\mathbf u}_{p-r+\ell}(n_{\mathfrak a})}{ \|{\boldsymbol \alpha}\|}.

Then, by (<ref>) and for the deterministic vector $\widetilde{ \boldsymbol{\gamma}}$ as in (<ref>),

\begin{equation}\label{e:P*(n_a)v_a(n_a)->gamma-tilde/const}
P^*(n_{\mathfrak a}){\mathbf v}_{\mathfrak a}(n_{\mathfrak a}) \stackrel \bbP\to \frac{\widetilde{\boldsymbol \gamma}}{\|\boldsymbol \alpha\|}, \quad n_{\mathfrak a} \rightarrow \infty,
\end{equation}

where $\boldsymbol \alpha := (0,\ldots,0,\alpha_q,\ldots,\alpha_r)^* \in\bbR^r$. Moreover, by (<ref>),

\frac{\lambda_{p-r+q}({\mathbf W}(a(n_{\mathfrak a})2^j))}{a(n_{\mathfrak a})^{2 h_q + 1}} = \inf_{{\mathbf u} \in \textnormal{span}\{{\mathbf u}_{p-r+q}(n_{\mathfrak a}),\ldots,{\mathbf u}_{p}(n_{\mathfrak a})\}\cap S^{p-1}}  {\mathbf u}^*\frac{{\mathbf W}(a(n_{\mathfrak a})2^j)}{a(n_{\mathfrak a})^{2 h_q + 1}}{\mathbf u}

\leq  {\mathbf v}^*_{\mathfrak a}(n_{\mathfrak a}) \frac{{\mathbf W}(a(n_{\mathfrak a})2^j)}{a(n_{\mathfrak a})^{2 h_q + 1}} {\mathbf v}_{\mathfrak a}(n_{\mathfrak a}).

Taking limits, expressions (<ref>), (<ref>) and (<ref>) imply that

\begin{equation}\label{e:alpha^2q+...+alpha^2r*xi_q,a<=_bound}
\|\boldsymbol \alpha\|^2 \xi_{q,\mathfrak a}(2^j)\leq  \widetilde{ \boldsymbol{\gamma}}^*\Big[\textnormal{diag}(\underbrace{0,\ldots,0}_{r_1},\underbrace{1,\ldots,1}_{r_2}) \B(2^j)\textnormal{diag}(\underbrace{0,\ldots,0}_{r_1},\underbrace{1,\ldots,1}_{r_2})\Big]\widetilde{ \boldsymbol{\gamma}}.
\end{equation}

Now, again in view of (<ref>), there exist scalars $\beta_{1}, \hdots, \beta_{q}$ such that

\begin{equation}\label{e:gammatilde_b-basis}
\widetilde{ \boldsymbol{\gamma}}= \sum^{q}_{\ell=1}\beta_{\ell}\hspace{0.5mm}{\boldsymbol \gamma}_{\ell}^{\mathfrak b} =: \Gamma_{{\mathfrak b}}\hspace{0.5mm}{\boldsymbol \beta} \in \bbR^{r}.
\end{equation}

So, along the subsequence $n_{\mathfrak b}$, define the (random) vector

{\mathbf v}_{\mathfrak b}(n_{\mathfrak b}) =  \frac{\sum^{q}_{\ell=1}\beta_{\ell}\hspace{0.5mm}{\mathbf u}_{p-r+\ell}(n_{\mathfrak b})}{\|\sum^{q}_{\ell=1}\beta_{\ell}\hspace{0.5mm}{\mathbf u}_{p-r+\ell}(n_{\mathfrak b})\|}
= \frac{\sum^{q}_{\ell=1}\beta_{\ell}\hspace{0.5mm}{\mathbf u}_{p-r+\ell}(n_{\mathfrak b})}{\|{\boldsymbol \beta}\|}.

Again by (<ref>) and for $\widetilde{ \boldsymbol{\gamma}}$ as in (<ref>),

\begin{equation}\label{e:P*(n_b)v_b(n_b)->gamma-tilde/const}
P^*(n_{\mathfrak b}) {\mathbf v}_{\mathfrak b}(n_{\mathfrak b}) \stackrel \bbP \to \frac{\widetilde{ \boldsymbol{\gamma}}}{\|\boldsymbol \beta\|},\quad n_{\mathfrak b} \rightarrow \infty,
\end{equation}

where $\boldsymbol \beta = (\beta_1,\ldots,\beta_q,0,\ldots,0)^*\in \bbR^r$.

In addition, again by (<ref>),

\frac{\lambda_{p-r+q}({\mathbf W}(a(n_{\mathfrak b})2^j))}{a(n_{\mathfrak b})^{2 h_q + 1}} = \sup_{{\mathbf u} \in \textnormal{span}\{{\mathbf u}_{p-r+1}(n_{\mathfrak b}),\ldots,{\mathbf u}_{p-r+q}(n_{\mathfrak b})\}\cap S^{p-1}}  {\mathbf u}^*\frac{{\mathbf W}(a(n_{\mathfrak b})2^j)}{a(n_{\mathfrak b})^{2 h_q + 1}}{\mathbf u}

\geq {\mathbf v}_{\mathfrak b}^*(n_{\mathfrak b}) \frac{{\mathbf W}(a(n_{\mathfrak b})2^j)}{a(n_{\mathfrak b})^{2 h_q + 1}} {\mathbf v}_{\mathfrak b}(n_{\mathfrak b}).

Taking limits, (<ref>) and (<ref>) imply that

\begin{equation}\label{e:beta_1^2+...+beta_q^2*xi_{q,b}(2^j)>=_bound}
\|\boldsymbol \beta\|^2 \xi_{q,\mathfrak b}(2^j)\geq  \widetilde{ \boldsymbol{\gamma}}^*\Big[\textnormal{diag}(\underbrace{0,\ldots,0}_{r_1},\underbrace{1,\ldots,1}_{r_2}) \B(2^j)\textnormal{diag}(\underbrace{0,\ldots,0}_{r_1},\underbrace{1,\ldots,1}_{r_2})\Big] \widetilde{ \boldsymbol{\gamma}}.
\end{equation}

As a consequence of (<ref>) and (<ref>),

\begin{equation}\label{e:sum_alpha*xi=<sum_beta*xi}
\|\boldsymbol \alpha\|^2 \xi_{q,\mathfrak a}(2^j)\leq  \|\boldsymbol \beta\|^2 \xi_{q,\mathfrak b}(2^j).
\end{equation}

Now, recall that $\Gamma_i$ is given by (<ref>). Also note that, by (<ref>),

\begin{equation}\label{e:<p1,p2>=c12}
\|P^*_H P^*(n)P(n)P_H -A \| = O(a(n)^{-\varpi}).
\end{equation}

Then, by conditions (<ref>) and (<ref>), for some $A \in {\mathcal S}_{>0}(r,\bbR)$,

\begin{equation}\label{e:Gamma-i_Gamma*_i=A}
\Gamma_i\Gamma_i^*=\plim_{n_i\to\infty} P^*(n_i)P(n_i) = A, \quad i=\mathfrak a,\mathfrak b.
\end{equation}

So, recast

\begin{equation}\label{e:Gammai=Vi*Sigmai*W-transpi}
\Gamma_i = V_i \Sigma_i W^*_i, \quad i = \mathfrak a,\mathfrak b,
\end{equation}

where $\Sigma_i \in GL(r,\bbR)$ is a diagonal matrix containing the singular values of $\Gamma_i$ and $V_i,W_i\in U(r)$, $i=\mathfrak a,\mathfrak b$. By (<ref>),

\begin{equation}\label{e:Va_Sigma^2_V*a=A=Vb_Sigma^2_V*b}
V_{\mathfrak a} \Sigma^2_{\mathfrak a} V^*_{\mathfrak a} = \Gamma_{\mathfrak a}\Gamma^*_{\mathfrak a} = A= \Gamma_{\mathfrak b}\Gamma^*_{\mathfrak b} = V_{\mathfrak b} \Sigma^2_{\mathfrak b} V^*_{\mathfrak b}.
\end{equation}

Recall that the diagonal entries of $\Sigma_i$, $i = a,b$, are nonnegative. Thus, based on the two spectral decompositions of $A$ in (<ref>), without loss of generality we can assume that

\begin{equation}\label{e:Va=Vb_Sigma-a=Sigma-b}
V:= V_{\mathfrak a} = V_{\mathfrak b}, \quad \Sigma := \Sigma_{\mathfrak a} = \Sigma_{\mathfrak b}.
\end{equation}

However, in view of (<ref>), (<ref>) and (<ref>),

\Gamma_{\mathfrak a}\hspace{0.5mm} \boldsymbol{\alpha} = \widetilde{\boldsymbol \gamma} = \Gamma_{\mathfrak b} \hspace{0.5mm}\boldsymbol{\beta}.

So, relations (<ref>) and (<ref>) imply that

\begin{equation}\label{e:alphanorm=betanorm}
\|\boldsymbol{\alpha}\|= \|\Gamma_{\mathfrak a}^{-1}\widetilde{\boldsymbol \gamma} \|= \| W_{\mathfrak a} \Sigma^{-1} V^* \widetilde{\boldsymbol \gamma}\| = \|\Sigma^{-1} V^* \widetilde{\boldsymbol \gamma}\|
=\| W_{\mathfrak b} \Sigma^{-1} V^* \widetilde{\boldsymbol \gamma}\|= \|\Gamma_{\mathfrak b}^{-1}\widetilde{\boldsymbol \gamma} \|=\|\boldsymbol{\beta}\|.
\end{equation}

By (<ref>), expression (<ref>) holds, as claimed. As anticipated, we can now use the symmetric roles played by $n_{\mathfrak a},n_{\mathfrak b}$ in the proof of (<ref>) to further conclude that (<ref>) holds. This completes the proof for the case $\mathcal{I}_+ = \emptyset$.

Now, suppose that $\mathcal{I}_+ \neq \emptyset$ and fix $q \in \{1,\hdots,r\}$. By relation (<ref>) (itself a consequence of Lemma <ref>),
\frac{\lambda_{p-r+q}({\mathbf W}(a2^j))}{a^{2h_q+1}}=O_{\bbP}(1).
So, it suffices to show that every subsequence of $\frac{\lambda_{p-r+q}({\mathbf W}(a2^j))}{a^{2h_q+1}}$ that is convergent in probability converges to the same limit. Proceeding as for the case $\mathcal{I}_+ = \emptyset$, suppose two such limits exist along some subsequences $n_{\mathfrak a}$, $n_{\mathfrak b}$ as in (<ref>). As in (<ref>) and (<ref>), the goal is to show that $\xi_{q,{\mathfrak a}}(2^j)\leq\xi_{q,{\mathfrak b}}(2^j)$ and then use the symmetry in the argument to conclude that $\xi_{q,{\mathfrak a}}(2^j)=  \xi_{q,{\mathfrak b}}(2^j)$.

So, analogously to the argument for the case $\mathcal{I}_+ = \emptyset$, by passing to further subsequences if necessary, expressions (<ref>) and (<ref>) hold. For $i=\mathfrak a,\mathfrak b$, let $ \gamma^i_{\ell,k}$ be the $k$–th entry of the limit vector ${\boldsymbol \gamma}_{\ell}^i$ in (<ref>), $\ell = 1,\hdots,r$. As a consequence of expression (<ref>) in Lemma <ref>, for $i=\mathfrak a,\mathfrak b$ and $k\in\mathcal{I}_+$, the $k$–th entries of the limit vector ${\boldsymbol \gamma}_{\ell}^i$ in (<ref>) are zero whenever $\ell \in \mathcal I_- \cup \mathcal I_0$. In other words, we can write

\ell \in \mathcal{I}_-\cup \mathcal{I}_0 =\{1,\ldots,r-r_3\}\Rightarrow {\boldsymbol \gamma}^i_\ell =
\Big(\gamma^i_{\ell,1}, \hdots, \gamma^i_{\ell,r_1+r_2}, 0,\hdots, 0\Big)^*,\quad i=\mathfrak a,\mathfrak b.

Thus, by Lemma <ref>,

\begin{equation}\label{e:w_samedim}
\textnormal{span}\{{\boldsymbol \gamma}_{1}^{\mathfrak a},\ldots,{\boldsymbol \gamma}_{r-r_3}^{\mathfrak a}\} = \textnormal{span}\{e_{\ell}, \ell \in \mathcal{I}_+\}^\perp=\textnormal{span}\{{\boldsymbol \gamma}_{1}^{\mathfrak b},\ldots,{\boldsymbol \gamma}_{r-r_3}^{\mathfrak b}\} \quad \textnormal{a.s.},
\end{equation}

where $e_\ell$ is the $\ell$-th standard Euclidean basis vector in $\bbR^r$. So, take any nonzero vector

\begin{equation}\label{e:gamma_intersection}
\widetilde{\boldsymbol \gamma} \in\textnormal{span}\{{\boldsymbol \gamma}_{q}^{\mathfrak a},\ldots,{\boldsymbol \gamma}_{r-r_3}^{\mathfrak a}\} \cap \textnormal{span}\{{\boldsymbol \gamma}_{1}^{\mathfrak b},\ldots,{\boldsymbol \gamma}_{q}^{\mathfrak b}\}
\end{equation}

The (almost sure) existence of $\widetilde{\boldsymbol \gamma}$ is guaranteed since $\dim \textnormal{span}\{{\boldsymbol \gamma}_{q}^{\mathfrak a},\ldots,{\boldsymbol \gamma}_{r-r_3}^{\mathfrak a}\}+ \dim\textnormal{span}\{{\boldsymbol \gamma}_{1}^{\mathfrak b},\ldots,{\boldsymbol \gamma}_{q}^{\mathfrak b}\}= r-r_3+1$, and the spanning sets in (<ref>) generate the same $(r-r_3)$-dimensional space. Now, for

\begin{equation}\label{e:gamma-vec_and_x-vec}
{\boldsymbol \gamma}=(\gamma_1,\ldots,\gamma_r)\in \bbR^{r} \textnormal{ and } \mathbf{x}=(x_{r-r_3+1},\ldots,x_r)\in\bbR^{r_3},
\end{equation}

define the vector

\begin{equation}\label{e:def_y(w,x)}
{\mathbf y}({\boldsymbol \gamma};\mathbf{x}) = \Big(\underbrace{0,\ldots,0}_{r_1},\underbrace{{\gamma}_{r_1+1},\ldots,{ \gamma}_{r_1+r_2}}_{r_2},\underbrace{x_{r-r_3+1},\ldots,x_{r}}_{r_3}\Big)^* \in \bbR^{r}.
\end{equation}

Based on (<ref>), define the (deterministic) function

\begin{equation}\label{e:def_g}
g({\boldsymbol \gamma};\mathbf{x}) :=
{\mathbf y}^*({\boldsymbol \gamma};\mathbf{x}) \B(2^j) {\mathbf y}({\boldsymbol \gamma};\mathbf{x}) \in \bbR.
\end{equation}

For fixed ${\boldsymbol \gamma}\in \bbR^{r}$, as a function of ${\mathbf x}$ we can reexpress $g({\boldsymbol \gamma}; \cdot)$ as

g({\boldsymbol \gamma}; \mathbf{x}) = \mathbf{x}^*Q\mathbf{x} + \langle \mathbf{c},\mathbf{x}\rangle +d,

where $Q = (b_{ii'})_{i,i'\in\mathcal{I}_+}$, $\mathbf{c} \in \bbR^{r_3}$ and $d$ are constants. Note that $Q = (b_{ii'})_{i,i'\in\mathcal{I}_+}\in \mathcal{S}_{>0}(r,\bbR)$, since $ \B(2^j)\in \mathcal{S}_{>0}(r,\bbR)$, i.e., the mapping $\mathbf{x}\mapsto\mathbf{x}^*Q\mathbf{x}$ is strictly convex. Consequently, for a fixed ${\boldsymbol \gamma}$, $g({\boldsymbol \gamma}; \cdot )$ is also strictly convex. Since, in addition, the function $g({\boldsymbol \gamma}; \cdot )$ is quadratic in the argument ${\mathbf x}$,  it has a unique global minimum. Also, for fixed ${\boldsymbol \gamma} \in \bbR^r$ and $\mathbf{x} \in\bbR^{r_3}$ as in (<ref>), define the vector

\begin{equation}\label{e:def_yn}
{\mathbf y}_n(\boldsymbol \gamma ;\mathbf{x})=\Big(\underbrace{\gamma_1 a(n)^{h_1-h_q},\ldots,\gamma_{r_1} a(n)^{h_{r_1}-h_q}}_{r_1},\underbrace{\gamma_{r_1+1},\ldots,\gamma_{r_1+r_2}}_{r_2},\underbrace{x_{r-r_3+1},\ldots,x_{r}}_{r_3}\Big)^*.
\end{equation}

Based on (<ref>), analogously to (<ref>), define the random function

\begin{equation}\label{e:def_gn}
g_n(\boldsymbol \gamma;\mathbf{x}):={\mathbf y}^*_n(\boldsymbol \gamma;\mathbf{x}) \widehat \B_a(2^j) {\mathbf y}_n(\boldsymbol \gamma;\mathbf{x}).
\end{equation}

\begin{equation}\label{e:x*gamma}
\mathbf{x}_{*,\widetilde {\boldsymbol \gamma}}= (x_{*,\widetilde {\boldsymbol \gamma},\ell})_{\ell\in\mathcal{I}_+}\in \bbR^{r_3}
\end{equation}

be the unique (deterministic) minimizer of the function $g(\widetilde {\boldsymbol \gamma};\cdot)$ as defined by (<ref>). By (<ref>), there are scalars $\alpha_q,\hdots,\alpha_{r-r_3}$ based on which we may reexpress

\begin{equation}\label{e:expand_gamma_alpha}
\widetilde{\boldsymbol \gamma} = \sum^{r-r_3}_{\ell=q}\alpha_\ell {\boldsymbol \gamma}_\ell^{\mathfrak a}  \in \bbR^r.
\end{equation}

For fixed $x_{*,\widetilde {\boldsymbol \gamma},\bullet}$ as in (<ref>), Lemma <ref> implies that, along the subsequence $n_{\mathfrak a}$, there is a sequence of unit vectors

\begin{equation}\label{e:va(na)}
{\mathbf v}_{\mathfrak a}(n_{\mathfrak a})\in \text{span}\{{\mathbf u}_{p-r+q}(n_{\mathfrak a}),\ldots,{\mathbf u}_{p}(n_{\mathfrak a})\}
\end{equation}

such that, with probability tending to 1 as $n_{\mathfrak a} \rightarrow \infty$, their angles with respect to the (non-unit) column vectors of $P(n_{\mathfrak a})$ satisfy

\begin{equation}\label{e:<p,w>=infinitesimal}
\langle \mathbf{p}_\ell(n_{\mathfrak a}), {\mathbf v}_{\mathfrak a}(n_{\mathfrak a}) \rangle = \frac{x_{*,\widetilde {\boldsymbol \gamma},\ell}}{a(n_{\mathfrak a})^{h_\ell-h_q}}, \quad \ell\in \mathcal{I}_+.
\end{equation}

Moreover, still by Lemma <ref>,

\widetilde{\boldsymbol \gamma}_{\mathfrak a}(n_{\mathfrak a}):= P^*(n_{\mathfrak a}) {\mathbf v}_{\mathfrak a}(n_{\mathfrak a}) \stackrel \bbP \to \frac{\widetilde{{\boldsymbol \gamma}}}{ \|\boldsymbol \alpha\| }, \quad n_{\mathfrak a} \rightarrow \infty,

where, for $\alpha_q,\hdots,\alpha_{r-r_3}$ as in (<ref>), we write

\boldsymbol \alpha = (0,\ldots,0,\alpha_q,\ldots,\alpha_{r-r_3},0,\ldots,0)^*  \in \bbR^r.

So, for each $n_{\mathfrak a}$
consider the function $g_{n_{\mathfrak a}}(\widetilde{\boldsymbol \gamma}_{\mathfrak a}(n_{\mathfrak a});\mathbf{x}_{*,\widetilde {\boldsymbol \gamma}})$ as defined by (<ref>). Note that, based on (<ref>) and (<ref>), we can reexpress

g_{n_{\mathfrak a}}(\widetilde{\boldsymbol \gamma}_{\mathfrak a}(n_{\mathfrak a});\mathbf{x}_{*,\widetilde {\boldsymbol \gamma}})
= {\mathbf v}^*_{\mathfrak a}(n_{\mathfrak a})P(n_{\mathfrak a}) \frac{a(n_{\mathfrak a})^{{\mathbf h}} \widehat{\mathbf  B}_a(2^j) a(n_{\mathfrak a})^{{\mathbf h}}}{a(n_{\mathfrak a})^{2h_q}} P^*(n_{\mathfrak a}){\mathbf v}_{\mathfrak a}(n_{\mathfrak a}).

Then, by (<ref>) and again with probability tending to 1,

\frac{\lambda_{p-r+q}(\mathbf{W}(a(n_{\mathfrak a})2^j))}{a(n_{\mathfrak a})^{2h_q + 1}}  = \inf_{{\mathbf u}\in \text{span}\{{\mathbf u}_{p-r+q}(n_{\mathfrak a}),\ldots,{\mathbf u}_{p}(n_{\mathfrak a})\}}{\mathbf u}^*\frac{{\mathbf W}(a(n_{\mathfrak a})2^j)}{a(n_{\mathfrak a})^{2h_q + 1}} {\mathbf u}

\begin{equation}\label{e:lambda_p-r+q<v*Wv}
\leq {\mathbf v}^*_{\mathfrak a}(n_{\mathfrak a}) \frac{ \mathbf{W}(a(n_{\mathfrak a})2^j))}{a(n_{\mathfrak a})^{2h_q + 1}}{\mathbf v}_{\mathfrak a}(n_{\mathfrak a})
\end{equation}

$$= g_{n_{\mathfrak a}}(\widetilde{\boldsymbol \gamma}_{\mathfrak a}(n_{\mathfrak a});\mathbf{x}_{*,\widetilde {\boldsymbol \gamma}}) +
{\mathbf v}^*_{\mathfrak a}(n_{\mathfrak a})
\Big( \frac{\W(a(n_{\mathfrak a})2^j)}{a(n_{\mathfrak a})^{2h_q+1}}-\frac{P(n_{\mathfrak a})a(n_{\mathfrak a})^{{\mathbf h}} \widehat{\mathbf  B}_a(2^j) a(n_{\mathfrak a})^{{\mathbf h}}P^*(n_{\mathfrak a})}{a(n_{\mathfrak a})^{2h_q}}\Big){\mathbf v}_{\mathfrak a}(n_{\mathfrak a}) .

However, by (<ref>),

a(n_{\mathfrak a})^{\mathbf h-h_q I}P^*(n_{\mathfrak a})\mathbf v_a(n_{\mathfrak a}) = {\mathbf y}_{n_{\mathfrak a}}(\widetilde{\boldsymbol \gamma}_{\mathfrak a}(n_{\mathfrak a});\mathbf{x}_{*,\widetilde {\boldsymbol \gamma}})\stackrel \bbP \to \|\boldsymbol \alpha\|^{-1} \mathbf y(\widetilde{\boldsymbol \gamma}, \mathbf{x}_{*,\widetilde {\boldsymbol \gamma}}), \quad n_{\mathfrak a} \rightarrow \infty.

In particular, $a(n_{\mathfrak a})^{\mathbf h-h_q I}P^*(n_{\mathfrak a})\mathbf v_a(n_{\mathfrak a}) =O_\bbP(1)$ along the subsequence $n_{\mathfrak a}$.  Thus, in view of expression (<ref>) and still along $n_{\mathfrak a}$,

{\mathbf v}_{\mathfrak a}^*(n_{\mathfrak a}) \Big( \W(a2^j)-\frac{Pa^{{\mathbf h}} \widehat{\mathbf  B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}}\Big){\mathbf v}_{\mathfrak a}(n_{\mathfrak a})

= o_\bbP(1) + \mathbf v_a^*(n_{\mathfrak a})\Big(\frac{Pa^{{\mathbf h}-h_q I}O_{\bbP}(1)}{a^{h_q+1/2}}+\frac{O_{\bbP}(1)a^{{\mathbf h}-h_q I}P^*}{a^{h_q+1/2}}\Big){\mathbf v}_{\mathfrak a}(n_{\mathfrak a})
= o_\bbP(1).

\begin{equation}\label{e:g_na(va,x*)_conv}
g_{n_{\mathfrak a}}(\widetilde{\boldsymbol \gamma}_{\mathfrak a}(n_{\mathfrak a});\mathbf{x}_{*,\widetilde {\boldsymbol \gamma}}) +
{\mathbf v}^*_{\mathfrak a}(n_{\mathfrak a})
\Big( \frac{\W(a2^j)}{a^{2h_q+1}}-\frac{Pa^{{\mathbf h}} \widehat{\mathbf  B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}}\Big){\mathbf v}_{\mathfrak a}(n_{\mathfrak a})
\stackrel \bbP \to \frac{g(\widetilde{{\boldsymbol \gamma}};\mathbf{x}_{*,\widetilde {\boldsymbol \gamma}}) }{ \|\boldsymbol \alpha\|^2 }, \quad n_{\mathfrak a} \rightarrow \infty.
\end{equation}

So, taking $\lim_{n_{\mathfrak a} \rightarrow \infty}$ in (<ref>), expressions (<ref>) and (<ref>) imply that

\begin{equation}\label{e:xi_a_leq_g}
\xi_{q,{\mathfrak a}}(2^j) \leq \frac{ g(\widetilde{{\boldsymbol \gamma}};\mathbf{x}_{*,\widetilde {\boldsymbol \gamma}}) }{ \|\boldsymbol \alpha\|^2 }.
\end{equation}

On the other hand, proceeding in the same fashion as with expression (<ref>), there are scalars $\beta_{1}, \hdots, \beta_{q}$ based on which we can reexpress

\begin{equation}\label{e:gamma=sum_beta*gamma-ell}
\widetilde{{\boldsymbol \gamma}}= \sum^{q}_{\ell=1}\beta_{\ell}{\boldsymbol \gamma}_{\ell}^{\mathfrak b}.
\end{equation}

Moreover, set

\begin{equation}\label{e:vb(nb)}
{\mathbf v}_{\mathfrak b}(n_{\mathfrak b}) =  \frac{\sum^{q}_{\ell=1}\beta_{\ell}{\mathbf u}_{p-r+\ell}(n_{\mathfrak b})}{\|\sum^{q}_{\ell=1}\beta_{\ell}{\mathbf u}_{p-r+\ell}(n_{\mathfrak b})\|}.
\end{equation}

Then, by (<ref>),

\begin{equation}\label{e:P*(nb)vb}
P^*(n_{\mathfrak b}){\mathbf v}_{\mathfrak b}(n_{\mathfrak b})\stackrel{\bbP}\to \frac{\widetilde{\boldsymbol \gamma}}{ \|\boldsymbol \beta\| }, \quad n_{\mathfrak b} \rightarrow \infty,
\end{equation}

where, for $\beta_1,\ldots,\beta_q$ as in (<ref>), we write

\boldsymbol \beta := (\beta_1,\ldots,\beta_q,0,\ldots,0)^* \in \bbR^r.

Moreover, as a consequence of relation (<ref>) in Lemma <ref>, by passing to a subsequence if necessary we can assume that there is ${\mathbf x}_0 \in \bbR^{r_3}$ such that

\begin{equation}\label{e:power_law*angle->x0/norm_beta}
\Big(a(n_{\mathfrak b})^{h_\ell - h_q} \langle {\mathbf p}_\ell(n_{\mathfrak b}), {\mathbf v}_{\mathfrak b}(n_{\mathfrak b}) \rangle\Big)_{\ell\in{\mathcal{I}_+}} \stackrel{\bbP}\rightarrow  \frac{{\mathbf x}_0}{ \|\boldsymbol \beta\| }, \quad n_b \rightarrow \infty.
\end{equation}

So, write $\widetilde{\boldsymbol \gamma}_{\mathfrak b}(n_{\mathfrak b})= P^*(n_{\mathfrak b}) {\mathbf v}_{\mathfrak b}(n_{\mathfrak b})$.  Then, by (<ref>),

\frac{\lambda_{p-r+q}({\mathbf W}(a(n_{\mathfrak b})2^j))}{a(n_{\mathfrak b})^{2 h_q + 1}} =  \sup_{{\mathbf u} \in \text{span}\{{\mathbf u}_{p-r+1}(n_{\mathfrak b}),\ldots,{\mathbf u}_{p-r+q}(n_{\mathfrak b})\}}{\mathbf u}^*\frac{{\mathbf W}(a(n_{\mathfrak b})2^j)}{a(n_{\mathfrak b})^{2 h_q + 1}}{\mathbf u}

\begin{equation}\label{e:lambda2/a^(2h1+1)_bounded_below}
\geq {\mathbf v}_{\mathfrak b}^*(n_{\mathfrak b}) \frac{{\mathbf W}(a(n_{\mathfrak b})2^j)}{a(n_{\mathfrak b})^{2 h_q + 1}} {\mathbf v}_{\mathfrak b}(n_{\mathfrak b})
\end{equation}

= g_{n_{\mathfrak b}}\Big(\widetilde{\boldsymbol \gamma}_{\mathfrak b}(n_{\mathfrak b});
\big(a(n_{\mathfrak b})^{h_\ell - h_q} \langle {\mathbf p}_\ell(n_{\mathfrak b}), {\mathbf v}_{\mathfrak b}(n_{\mathfrak b}) \rangle\big)_{\ell\in{\mathcal{I}_+}}\Big)
+ {\mathbf v}_{\mathfrak b}^*(n_{\mathfrak b}) \Big( \W(a2^j)-\frac{Pa^{{\mathbf h}} \widehat{\mathbf  B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}}\Big){\mathbf v}_{\mathfrak b}(n_{\mathfrak b}).

However, by expressions (<ref>) and (<ref>) – and still along the subsequence $n_{\mathfrak b}$ –, the right-hand side of (<ref>) can be rewritten as

{\mathbf y}^*_n\Big(\widetilde{\boldsymbol \gamma}_{\mathfrak b}(n_{\mathfrak b});
\big(a(n_{\mathfrak b})^{h_\ell - h_q} \langle {\mathbf p}_\ell(n_{\mathfrak b}), {\mathbf v}_{\mathfrak b}(n_{\mathfrak b}) \rangle\big)_{\ell\in{\mathcal{I}_+}}\Big) \widehat \B_a(2^j) {\mathbf y}_n\Big(\widetilde{\boldsymbol \gamma}_{\mathfrak b}(n_{\mathfrak b});
\big(a(n_{\mathfrak b})^{h_\ell - h_q} \langle {\mathbf p}_\ell(n_{\mathfrak b}), {\mathbf v}_{\mathfrak b}(n_{\mathfrak b}) \rangle\big)_{\ell\in{\mathcal{I}_+}}\Big)

\begin{equation}\label{e:g_nb(vb,x*)_conv}
+ o_{\bbP}(1) + {\mathbf v}^*_{\mathfrak b}(n_{\mathfrak b})
\Big(\frac{P a^{{\mathbf h}-h_qI}O_{\bbP}(1)}{a^{h_q+1/2}}+\frac{O_{\bbP}(1)a^{{\mathbf h}-h_qI}P^*}{a^{h_q+1/2}}\Big){\mathbf v}_{\mathfrak b}(n_{\mathfrak b})
\stackrel{\bbP} \rightarrow \frac{ g( \widetilde{{\boldsymbol \gamma}};\mathbf{x}_{0} )}{ \|\boldsymbol \beta\|^2 } \geq  \frac{ g ( \widetilde{{\boldsymbol \gamma}};\mathbf{x}_{*,\widetilde{\boldsymbol \gamma}})}{ \|\boldsymbol \beta\|^2 },
\end{equation}

as $n_b \rightarrow \infty$, where $\mathbf{x}_{*,\widetilde{\boldsymbol \gamma}}$ is defined as in (<ref>). In (<ref>), the limit is a consequence of relations (<ref>) and (<ref>). By taking $\lim_{n_{\mathfrak b}\to\infty}$ in (<ref>), expressions (<ref>) and (<ref>) imply that
\xi_{q,{\mathfrak b}}(2^j) \geq \frac{g(\widetilde {\boldsymbol \gamma};\mathbf{x}_{*,\widetilde {\boldsymbol \gamma}})}{ \|\boldsymbol \beta\|^2  }.
By repeating the argument leading to (<ref>), we find that $\|\boldsymbol \alpha\| = \|\boldsymbol \beta\|$, implying
\xi_{q,{\mathfrak b}}(2^j)\geq \frac{g(\widetilde {\boldsymbol \gamma};\mathbf{x}_{*,\widetilde {\boldsymbol \gamma}})}{\|\boldsymbol \alpha\|^2} \geq \xi_{q,{\mathfrak a}}(2^j).
Analogously, by exchanging the roles of the subsequences $n_{\mathfrak a}$ and $n_{\mathfrak b}$ in the argument, we find that

\begin{equation}\label{e:xi_q=g(gamma,x)}
\xi_{q,{\mathfrak a}}(2^j)= \frac{g(\widetilde {\boldsymbol \gamma},\mathbf{x}_{*,\widetilde {\boldsymbol \gamma}})}{\|\boldsymbol \alpha\|^2}= \xi_{q,{\mathfrak b}}(2^j). \hspace{2mm}
\end{equation}

This finishes the proof for the case where ${\mathcal I}_+ \neq \emptyset$. Thus, (<ref>) is established. $\Box$

Proof of Theorem <ref>:
Recall that, in this section, for notational simplicity we work under (<ref>), namely, $P(n)P_H \equiv P(n) \equiv P$. In view of (<ref>), for notational simplicity we write

\begin{equation}\label{e:WZ=OP(1),a^(-h-1/2I)WXZ}
{\mathbf W}_Z(a(n)2^{j})=O_{\bbP}(1), \quad a(n)^{-{\mathbf h}-(1/2)I}{\mathbf W}_{X,Z}(a(n)2^{j})=O_{\bbP}(1),
\end{equation}

where the appropriate matrix dimensions are implicit. To avoid needless duplication in proving the result under conditions (<ref>) and (<ref>) separately, we also continue to make use of the notation (<ref>). So, rewrite

\log \lambda_{p-r+q}\Big( \frac{{\mathbf W}(a(n)2^{j})}{a^{2h_q+1}} \Big)  - \log \lambda_{p-r+q}\Big( \frac{\bbE {\mathbf W} (a(n)2^{j})}{a^{2h_q+1}} \Big)

= \Big\{\log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} \widehat{{\mathbf B}}_a(2^j)a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{O_{\bbP}(1)}{a^{2h_q+1}}+
\frac{Pa^{{\mathbf h}}O_{\bbP}(1)}{a^{2h_q+1/2}}+\frac{O_{\bbP}(1)a^{{\mathbf h}}P^*}{a^{2h_q+1/2}}\Big)

- \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{O_{\bbP}(1)}{a^{2h_q+1}}+
\frac{Pa^{{\mathbf h}}O_{\bbP}(1)}{a^{2h_q+1/2}}+\frac{O_{\bbP}(1)a^{{\mathbf h}}P^*}{a^{2h_q+1/2}}\Big) \Big\}

+ \Big\{\log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{O_{\bbP}(1)}{a^{2h_q+1}}+
\frac{Pa^{{\mathbf h}}O_{\bbP}(1)}{a^{2h_q+1/2}}+\frac{O_{\bbP}(1)a^{{\mathbf h}}P^*}{a^{2h_q+1/2}}\Big)

- \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{O_{\bbP}(1)}{a^{2h_q+1}}\Big) \Big\}

\begin{equation}\label{e:loglambda-loglambdaE_3_terms}
+ \Big\{\log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{O_{\bbP}(1)}{a^{2h_q+1}}\Big)
- \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} \Big)\Big\}.
\end{equation}

We consider each term on the right-hand side of (<ref>) separately. Note that there is functional dependence between the matrices $\widehat{{\mathbf B}}_a(2^j) \in {\mathcal S}_{\geq 0}(r,\bbR)$, $O_{\bbP}(1) \in {\mathcal S}_{\geq 0}(p,\bbR)$ and $O_{\bbP}(1) \in M(r,p,\bbR)$. However, we will construct partial mean value theorem-type expansions of each term in (<ref>) where each of these matrix arguments is first treated as functionally independent variables, and then plug back in the actual value of the matrix argument.

For ${\mathbf B} \in {\mathcal S}_{> 0}(r,\bbR)$, define the functions

\begin{equation}\label{e:fv}
f_{n,q,1}({\mathbf B}) := \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{O_{\bbP}(1)}{a^{2h_q+1}}+
\frac{Pa^{{\mathbf h}}O_{\bbP}(1)}{a^{2h_q+1/2}}+\frac{O_{\bbP}(1)a^{{\mathbf h}}P^*}{a^{2h_q+1/2}}\Big),
\end{equation}

\begin{equation}\label{e:fv_2}
f_{n,q,2}(R)  := \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{O_{\bbP}(1)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} R + R^*\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)
\end{equation}

\begin{equation}\label{e:fv_3}
f_{n,q,3}(R):= \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} +R\Big).
\end{equation}

The layout of the proof is as follows.  We will establish the  convergence

\begin{equation}\label{e:f_1_weakconv}
\Big(\sqrt{n_{a,j}}\hspace{0.5mm}\Big(f_{n,q,1}(\widehat{{\mathbf B}}_a(2^j)) - f_{n,q,1}({\mathbf B}_a(2^j))\Big)_{q=1,\hdots,r}\Big)_{j=j_1,\hdots,j_2}\stackrel{d}\rightarrow {\mathcal N}(0,\Sigma_{\lambda}),
\end{equation}

and also that

\begin{equation}\label{e:f_2_to0}
\sqrt{\frac{n}{a}}\Big(f_{n,q,2}\Big( \frac{a^{- {\mathbf h}- (1/2)I}{\mathbf W}_{X,Z}(a2^j)}{a^{h_q + 1/2}} \Big) - f_{n,q,2}(
\mathbf{0})\Big) = o_{\bbP}(1),
\end{equation}

\begin{equation}\label{e:f_3_to0}
\sqrt{\frac{n}{a}}\Big(f_{n,q,3}\Big( \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}}\Big) -f_{n,q,3}({\mathbf 0})\Big) = o_{\bbP}(1).
\end{equation}

Then, as a consequence of (<ref>), (<ref>) and (<ref>),

\Big( \sqrt{n_{a,j}}\Big( \log \lambda_{p-r+q}({\mathbf W}(a(n)2^{j}))  - \log \lambda_{p-r+q}(\bbE {\mathbf W} (a(n)2^{j})) \Big)_{q=1,\hdots,r} \Big)_{j=j_1,\hdots,j_2}

=  \Big(\sqrt{n_{a,j}}\hspace{0.5mm}\Big(f_{n,q,1}(\widehat{{\mathbf B}}_a(2^j)) - f_{n,q,1}({\mathbf B}_a(2^j))\Big)_{q=1,\hdots,r}\Big)_{j=j_1,\hdots,j_2}+ o_{\bbP}(1) \stackrel{d}\rightarrow {\mathcal N}(0,\Sigma_{\lambda}),

as $n \rightarrow \infty$, which proves (<ref>).

We proceed first to establish (<ref>).  Recall that, for any $M \in {\mathcal S}(p,\bbR)$, the differential of a simple eigenvalue $\lambda_{q}(M)$ exists in a connected vicinity of $M$ and is given by

\begin{equation}\label{e:dlambdal}
d\lambda_{q}(M) = {\mathbf u}_{q}^* \hspace{0.5mm}\{d M \} \hspace{0.5mm}{\mathbf u}_{q},
\end{equation}

where ${\mathbf u}_{q}$ is a unit eigenvector of $M$ associated with $\lambda_{q}(M)$ (Magnus <cit.>, p. 182, Theorem 1). Now, fix an arbitrary $q\in \{1,\ldots,r\}$, and let $\mathcal{I}_0$ and (the possibly empty) sets $\mathcal{I}_-,\mathcal{I}_+$ be as in (<ref>). The eigenvalue

\begin{equation}\label{e:EW_has_pairwise_distinct_eigens}
\lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}}  + \frac{O_{\bbP}(1)}{a^{2h_q+1}}+
\frac{Pa^{{\mathbf h}}O_{\bbP}(1)}{a^{2h_q+1/2}}+\frac{O_{\bbP}(1)a^{{\mathbf h}}P^*}{a^{2h_q+1/2}} \Big)
\end{equation}

must be simple for large enough $n$ with probability tending to 1. In fact, for $M_n = a^{2h_q+1}\big(\frac{O_{\bbP}(1)}{a^{2h_q+1}}+
\frac{Pa^{{\mathbf h}}O_{\bbP}(1)}{a^{2h_q+1/2}}+\frac{O^*_{\bbP}(1)a^{{\mathbf h}}P^*}{a^{2h_q+1/2}}\big) $ (see expression (<ref>)), by Corollary <ref>,

\lambda_{p-r+\ell}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}}  + \frac{O_{\bbP}(1)}{a^{2h_q+1}}+
\frac{Pa^{{\mathbf h}}O_{\bbP}(1)}{a^{2h_q+1/2}}+\frac{O^*_{\bbP}(1)a^{{\mathbf h}}P^*}{a^{2h_q+1/2}} \Big)
\stackrel{\bbP}\to
\begin{cases}
0, & \ell\in \mathcal I_-;\\
\xi_\ell(2^j), & \ell\in \mathcal I_0;\\
\infty, & \ell\in \mathcal I_+,
\end{cases}

as $n \rightarrow \infty$, where $\xi_\ell(2^j)$, $\ell \in \mathcal{I}_0$, are given by (<ref>) and are pairwise distinct by condition (<ref>). So, for any fixed (and small enough) $\varepsilon > 0$, let
{\mathcal O}_{n} = \{ {\mathbf B} \in {\mathcal S}_{> 0}(r,\bbR): \|{\mathbf B}-{\mathbf B}_a(2^j)\| < \varepsilon\}.
Then, also for large $n$, by (<ref>) with probability going to 1 the derivative of the function $f_{n,q,1}$ in (<ref>) exists in the open and connected set ${\mathcal O}_n  \subseteq {\mathcal S}_{> 0}(r,\bbR)$. For any ${\mathbf B} \in {\mathcal O}_n$, an application of Proposition 3 in Abry and Didier <cit.> yields

\begin{equation}\label{e:Taylor_determ}
f_{n,q,1}({\mathbf B}) - f_{n,q,1}({\mathbf B}_a(2^j)) = \sum^{r}_{\ell = 1}\sum^{r}_{\ell' = 1}\frac{\partial}{\partial b_{\ell,\ell'}}f_{n,q}(\breve{{\mathbf B}}) \hspace{1mm} \pi_{\ell,\ell'}({\mathbf B} - {\mathbf B}_a(2^j))
\end{equation}

for some matrix $\breve{{\mathbf B}} \in {\mathcal S}_{>0}(r,\bbR)$ lying in a segment connecting ${\mathbf B}$ and ${\mathbf B}_a(2^j)$ across ${\mathcal  S}_{>0}(r,\bbR)$. Define the event $A_n = \big\{ \omega:  \widehat{{\mathbf B}}_{a}(2^j) \in {\mathcal O}_n \big\}$. By (<ref>),

\begin{equation}\label{e:W_has_pairwise_distinct_eigens}
\bbP (A_n) \rightarrow 1, \quad n \rightarrow \infty.
\end{equation}

By (<ref>), for large enough $n$ and in the set $A$, the expansion

\begin{equation}\label{e:Taylor}
f_{n,q,1}(\widehat{{\mathbf B}}_{a} (2^j)) - f_{n,q,1}({\mathbf B}_a(2^j)) = \sum^{r}_{\ell = 1}\sum^{r}_{\ell' = 1}\frac{\partial}{\partial b_{\ell,\ell'}}f_{n,q,1}(\breve{{\mathbf B}}_{a} (2^j)) \hspace{1mm} \pi_{\ell,\ell'}(\widehat{{\mathbf B}}_{a}(2^j) - {\mathbf B}_a(2^j))
\end{equation}

holds for some matrix $\breve{{\mathbf B}}_{a}(2^j)$ lying in a segment connecting $\widehat{{\mathbf B}}_{a}(2^j)$ and ${\mathbf B}_a(2^j)$ across ${\mathcal  S}_{>0}(r,\bbR)$.

Consider the matrix of derivatives

\Big\{\frac{\partial}{\partial b_{\ell,\ell'}}f_{n,q,1}(\breve{{\mathbf B}})\Big\}_{\ell,\ell'= 1,\hdots,r}
= \Big\{ \lambda^{-1}_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} \breve{{\mathbf B}} a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{O_{\bbP}(1)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \frac{O_{\bbP}(1)}{a^{h_q+1/2}} + \frac{O^*_{\bbP}(1)}{a^{h_q+1/2}}\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)

\begin{equation}\label{e:df1/dbi1,i2}
\times \frac{\partial}{\partial b_{\ell,\ell'}}\lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} \breve{{\mathbf B}} a^{{\mathbf h}}P^*}{a^{2h_q}} +\frac{O_{\bbP}(1)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \frac{O_{\bbP}(1)}{a^{h_q+1/2}} + \frac{O^*_{\bbP}(1)}{a^{h_q+1/2}}\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)\Big\}_{\ell,\ell'= 1,\hdots,r}.
\end{equation}

In (<ref>), the differential of the eigenvalue
\lambda_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} \breve{{\mathbf B}} a^{{\mathbf h}}P^*}{a^{2h_q}} +\frac{O_{\bbP}(1)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \frac{O_{\bbP}(1)}{a^{h_q+1/2}} + \frac{O^*_{\bbP}(1)}{a^{h_q+1/2}}\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)
is given by expression (<ref>) with

\begin{equation}\label{e:MVT_arg_fnq1}
\frac{\breve{{\mathbf W}}}{a^{2h_q+1}}:= \frac{Pa^{{\mathbf h}} \breve{{\mathbf B}} a^{{\mathbf h}}P^*}{a^{2h_q}} +\frac{O_{\bbP}(1)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \frac{O_{\bbP}(1)}{a^{h_q+1/2}} + \frac{O^*_{\bbP}(1)}{a^{h_q+1/2}} \frac{a^{{\mathbf h}}P^*}{a^{h_q}}
\end{equation}

in place of $M$ and ${\mathbf u}_{p-r+q}(n)$ denoting a unit eigenvector of (<ref>) associated with its $(p-r+q)$–th eigenvalue. Moreover,

\frac{\partial}{\partial b_{\ell,\ell'}}\lambda_{p-r+q}\Big(\frac{\breve{{\mathbf W}}}{a^{2h_q+1}}\Big)

= {\mathbf u}^*_{p-r+q}(n)\frac{\partial}{\partial b_{\ell,\ell'}}\Big[\frac{Pa^{{\mathbf h}} \breve{{\mathbf B}} a^{{\mathbf h}}P^*}{a^{2h_q}} +\frac{O_{\bbP}(1)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \frac{O_{\bbP}(1)}{a^{h_q+1/2}} + \frac{O^*_{\bbP}(1)}{a^{h_q+1/2}}\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big]{\mathbf u}_{p-r+q}(n)

= {\mathbf u}^*_{p-r+q}(n)\Big( \frac{Pa^{{\mathbf h}} {\mathbf 1}_{\ell,\ell'} a^{{\mathbf h}}P^*}{a^{2h_q}}\Big){\mathbf u}_{p-r+q}(n), \quad \quad  \ell,\ell'= 1,\hdots,r,

where ${\boldsymbol 1}_{\ell,\ell'}$ is a matrix with 1 on entry $(\ell,\ell')$ and zeroes elsewhere. Therefore, by Proposition <ref>, for $1 \leq \ell \leq \ell' \leq r$, and by using $\breve{{\mathbf B}}_{a}(2^j)$ in place of $\breve{{\mathbf B}}$ and $\breve{{\mathbf W}}(a2^j)$ in place of $\breve{{\mathbf W}}$ in (<ref>), we can pick the sequence ${\mathbf u}_{p-r+q}(n)$ so as to obtain

a^{-(2h_{q}+1)} {\mathbf u}^*_{p-r+q}(n)\Big\{\frac{\partial}{\partial b_{\ell,\ell'}}\breve{{\mathbf W}}(a2^j)\Big\}{\mathbf u}_{p-r+q}(n)

= {\mathbf u}^*_{p-r+q}(n)P \textnormal{diag}(a^{h_1 - h_{q}},\hdots,1,\hdots,a^{h_n - h_{q}}){\boldsymbol 1}_{\ell,\ell'}
\textnormal{diag}(a^{h_1 - h_{q}},\hdots,1,\hdots,a^{h_n - h_{q}})P^*{\mathbf u}_{p-r+q}(n)

\begin{equation}\label{e:limit_derivative_eigenvalue}
=\langle {\mathbf p}_{ \ell}(n),{\mathbf u}_{p-r+q}(n)\rangle a^{h_{\ell} - h_{q}} \langle {\mathbf p}_{ \ell'}(n),{\mathbf u}_{p-r+q}(n)\rangle a^{h_{\ell'} - h_{q}}
\stackrel{\bbP}\rightarrow \left\{\begin{array}{cc}
0, & \ell \in\mathcal{I}_-;\\
\gamma_{{\ell}q }\gamma_{\ell'q}, & \ell,\ell'\in\mathcal{I}_0,\\
\gamma_{\ell q} x_{\ell',*}, & \ell \in \mathcal{I}_0, \ell'\in \mathcal{I}_+;\\
x_{\ell,*}x_{\ell',*}, &  \ell,\ell'\in \mathcal{I}_+.\\
\end{array}\right.
\end{equation}

In (<ref>), the entries $x_{\ell,*}$ (depending on $q$), $\ell \in \mathcal{I}_+$, of the vector $\textbf{x}_{q,*}(2^j)$ are given by expression (<ref>) in Lemma <ref>, and each $\gamma_{\ell,q}$ is as in Proposition <ref>. Then, by (<ref>) and (<ref>) in Theorem <ref>, expression (<ref>) converges in probability to the matrix

\begin{equation}\label{e:limit_derivative}
\frac{1}{\xi_q(2^j)}\times
\left(\begin{array}{cccccc}
0 &0& 0 & \\
0& (\gamma_{{\ell}q }\gamma_{\ell'q})_{\ell,\ell'\in\mathcal{I}_0} & (\gamma_{\ell q} x_{\ell',*})_{\ell \in \mathcal{I}_0, \ell'\in \mathcal{I}_+}\\
0 &  (x_{\ell,*}\gamma_{\ell' q} )_{ \ell \in \mathcal{I}_+,\ell' \in \mathcal{I}_0}& (x_{\ell,*}x_{\ell',*})_{\ell,\ell'\in \mathcal{I}_+}
\end{array}\right),
\end{equation}

where $\xi_q(2^j) > 0$. Turning back to (<ref>), expression (<ref>) and condition (<ref>) imply that

\bbR \ni \sqrt{n_{a,j}}\hspace{0.5mm}(f_{n,q,1}(\widehat{{\mathbf B}}_a(2^j)) - f_{n,q,1}({\mathbf B}_a(2^j)))

\begin{equation}\label{e:term1_conv_in_law}
= \sum^{r}_{\ell = 1}\sum^{r}_{\ell' = 1}\frac{\partial}{\partial b_{\ell,\ell'}}f_{n,q,1}(\breve{{\mathbf B}}) \hspace{1mm} \sqrt{n_{a,j}}\hspace{0.5mm}\pi_{\ell,\ell'}(\widehat{{\mathbf B}}_a(2^j) - {\mathbf B}_a(2^j)) \stackrel{d}\rightarrow  {\mathcal N}(0,\sigma^2_{q}(j)),
\end{equation}

as $n \rightarrow \infty$, where $\sigma^2_{q}(j) > 0$ as a consequence of the fact that $\Sigma_B(j)$ in (<ref>) has full rank. Since, for $j = j_1,\hdots,j_2$ and $q = 1,\hdots,r$, the asymptotic normality of each individual log-eigenvalue results from the factor (<ref>), then the joint asymptotic normality (<ref>) holds, as claimed.

We now turn to (<ref>). The eigenvalue $\log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} \Big)$
must be simple for large enough $n$ with probability tending to 1. In fact, by Corollary <ref> with $M_n \equiv {\mathbf 0}$,

\lambda_{p-r+\ell}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}}\Big) \to
\begin{cases}
0, & \ell \in {\mathcal I}_-;\\
\xi_\ell(2^j), & \ell \in {\mathcal I}_0;\\
\infty, & \ell \in {\mathcal I}_+,
\end{cases}

where the functions $\xi_\ell(2^j)$ are given by (<ref>) and are pairwise distinct by condition (<ref>) when $h_1 = \hdots = h_r$. In other words, for any small $\delta>0$, there is

\begin{equation}\label{e:n0}
n_0 \in \bbN
\end{equation}

such that, for $n \geq n_0$, the mapping

M(p,\bbR) \ni R\mapsto\lambda_{p-r+q}\left(\frac{Pa^{{\mathbf h}}\mathbf{B}_a(2^j)a^{{\mathbf h}}P^*}{a^{2h_q}} + R\right)

is differentiable for all $\|R\|<\delta$. Therefore, by Proposition 3 in Abry and Didier <cit.>, we can write

f_{n,q,3}(R) - f_{n,q,3}({\mathbf 0}):= \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + R\Big)
- \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} \Big)

\begin{equation}\label{e:fnq3(R)-fnq3(0)}
= \sum^{p}_{i_1 = 1}\sum^{p}_{i_2 = 1}\frac{\partial}{\partial r_{i_1,i_2}}f_{n,q,1}(\breve{R}) \hspace{1mm} \pi_{i_1,i_2}(R)
\end{equation}

for some matrix $\breve{R} \in {\mathcal S}_{>0}(p,\bbR)$ lying in a segment connecting $R$ and ${\mathbf 0}$ across ${\mathcal S}_{>0}(p,\bbR)$. Moreover,

\Big\{\frac{\partial}{\partial r_{i_1,i_2}}f_{n,q,3}(\breve{R})\Big\}_{i_1,i_2= 1,\hdots,p}
= \Big\{ \lambda^{-1}_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \breve{R} \Big)

\begin{equation}\label{e:df/dbi1,i2}
\times \frac{\partial}{\partial r_{i_1,i_2}}\lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \breve{R} \Big)\Big\}_{i_1,i_2= 1,\hdots,p},
\end{equation}

where the differential of the eigenvalue $\lambda_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + R \Big)$
is given by expression (<ref>) with

\begin{equation}\label{e:MVT_arg_fnq3}
\frac{Pa^{{\mathbf h}}{\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + R
\end{equation}

in place of $M$ and ${\mathbf u}_{p-r+q}(n)$ denoting a unit eigenvector of (<ref>) associated with its $(p-r+q)$–th eigenvalue. In addition,

\frac{\partial}{\partial r_{i_1,i_2}}\lambda_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + R\Big)
= {\mathbf u}^*_{p-r+q}(n)\frac{\partial}{\partial r_{i_1,i_2}}\Big[\frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + R\Big]{\mathbf u}_{p-r+q}(n)

= {\mathbf u}^*_{p-r+q}(n) {\mathbf 1}_{i_1,i_2}{\mathbf u}_{p-r+q}(n), \quad \quad  i_1, i_2= 1,\hdots,p.

\begin{equation}\label{e:varphi1>0}
\varphi_1 > 0.
\end{equation}

Under condition (<ref>), for any small $\delta>0$, there is $n_1 \in \bbN$ such that, for $n \geq n_1$, $\|{\mathbf W}_Z(a2^j)/a^{2h_q + 1}\| \leq \delta$ with probability $1 - \varphi_1$. So, let $n_0$ be as in (<ref>). Then, for $n \geq \max\{n_0,n_1\}$, expression (<ref>) implies that

\sqrt{\frac{n}{a}}\Big(f_{n,q,3}\Big( \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}}\Big) -f_{n,q,3}({\mathbf 0})\Big)
= \sum^{p}_{i_1 = 1}\sum^{p}_{i_2 = 1}\frac{\partial}{\partial r_{i_1,i_2}}f_{n,q,3}(\breve{R}_n) \hspace{1mm} \sqrt{\frac{n}{a}}\pi_{i_1,i_2}\Big( \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} \Big)

\begin{equation}\label{e:sqrt(n/a)(fnq3(WZ/a^2hq+1)-fnq(0))}
= \lambda^{-1}_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \breve{R}_n\Big)\sum^{p}_{i_1 = 1}\sum^{p}_{i_2 = 1} {\mathbf u}^*_{p-r+q}(n) {\mathbf 1}_{i_1,i_2}{\mathbf u}_{p-r+q}(n) \hspace{1mm} \sqrt{\frac{n}{a}}\pi_{i_1,i_2}\Big( \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} \Big)
\end{equation}

with probability $1 - \varphi_1$. In (<ref>), the matrix $\breve{R}_n$ lies between $O_{\bbP}(1)/a^{2h_q+1}$ and $\mathbf{0}$. Since $\|\breve{R}_n\|=o_{\bbP}(1)$, by Corollary <ref>, $\lambda_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \breve{R}_n\Big) \stackrel{\bbP}\rightarrow \xi_{q}(2^j)$, $n \rightarrow \infty$. Now recall that, for $M = (m_{i_1,i_2}) \in M(p_1,p_2,\bbR)$,

\begin{equation}\label{e:max|sil(p)|=<C'}
\max_{i_1=1,\hdots,p_1; \hspace{.5mm}i_2 =1,\hdots,p_2}|m_{i_1, i_2}(n)| \leq \| M\|.
\end{equation}

Thus, by condition (<ref>), $\max_{i_1,i_2=1,\hdots,p}|\pi_{i_1,i_2}( {\mathbf W}_Z(a2^j) )| = O_{\bbP}(1)$. Also recall that, for a vector $\textbf{x} \in \bbR^p$,
\begin{equation}\label{e:|x|1=<sqrt(p)|x|2}
\|\textbf{x}\|_{1}\leq \sqrt{p} \hspace{0.5mm}\|\textbf{x}\|_{2}.
\end{equation}

Then, expression (<ref>) is bounded in absolute value by

\frac{C}{a^{2h_q+1}}\sqrt{\frac{n}{a}} \sum^{p}_{i_1 = 1}\sum^{p}_{i_2 = 1} |{\mathbf u}_{p-r+q}(n)_{i_1}{\mathbf u}_{p-r+q}(n)_{i_2}| =
\frac{C}{a^{2h_q+1}}\sqrt{\frac{n}{a}} \hspace{0.5mm}\|{\mathbf u}_{p-r+q}(n)\|^2_{1}

\begin{equation}\label{e:term3_bound}
\leq \frac{C}{a^{2h_q+1}}\sqrt{\frac{n}{a}} \hspace{0.5mm}p = \frac{C}{a^{2h_q+1}}\Big(\frac{n}{a}\Big)^{3/2} \frac{p}{n/a} = C \frac{a^{3h_1}}{a^{2h_q+1}}\Big(\frac{n}{a^{2h_1 +1}}\Big)^{3/2} \frac{p}{n/a} \rightarrow 0, \quad n \rightarrow \infty,
\end{equation}

where the limit follows from condition (<ref>). Since $\varphi_1 > 0$ is arbitrary,
\sqrt{\frac{n}{a}}\Big(f_{n,q,3}\Big( \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}}\Big) -f_{n,q,3}({\mathbf 0})\Big) = o_{\bbP}(1),
which is (<ref>).

We now turn to (<ref>). Note that, again by Corollary <ref> with $M_n = {\mathbf W}_Z(a2^j) = O_{\bbP}(1)$ (see (<ref>)),

\begin{equation}\label{e:lambda(WX+WZ)_simple}
\lambda_{p-r+\ell}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}}\Big) \stackrel{\bbP}\to
\begin{cases}
0, & \ell \in {\mathcal I}_-;\\
\xi_\ell(2^j), & \ell \in {\mathcal I}_0;\\
\infty, & \ell \in {\mathcal I}_+.
\end{cases}
\end{equation}

In (<ref>), the functions $\xi_\ell(2^j)$, $\ell \in \mathcal{I}_0$, are given by (<ref>) and are pairwise distinct by condition (<ref>) when $h_1 = \hdots = h_r$. For $R \in {\mathcal M}(r,p,\bbR)$, recast

f_{n,q,2}(R)  = \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} R + R^*\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)

= \log \Big(\hspace{0.5mm} {\mathbf u}^*_{p-r+q}(n,R)\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} R + R^*\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big){\mathbf u}_{p-r+q}(n,R)\hspace{0.5mm}\Big)

\begin{equation}\label{e:log_lambda(W(a2j,R))/a^(2hq+1)}
=: \log \lambda_{p-r+q}(\boldsymbol{\mathcal W}(a(n)2^j,R)/a^{2h_q+1}).
\end{equation}

We claim that for a sufficiently small $\varepsilon>0$, the eigenvalues

\begin{equation}\label{e:lambda(WX+W_Z+R-term)_simple_claim}
\lambda_{p-r+\ell}(\boldsymbol{\mathcal W}(a(n)2^j,R)), \quad \ell \in1,\ldots r,
\end{equation}

must be simple for large enough $n$ with probability tending to 1 uniformly over $\|R\|\leq \varepsilon$ for all $\ell\in\mathcal{I}_0$. Indeed, continuity of eigenvalues implies that for every $\varepsilon>0$, $R\mapsto \lambda_{p-r+\ell}(\boldsymbol{\mathcal W}(a(n)2^j,R)) $ attains its maximum over $\{R:\|R\|\leq \varepsilon\}$ at some $R_n(\ell)$ for each $n$.  Clearly, $R_n(\ell)=O_\bbP(1),$  so Corollary <ref> with $M_n= a^{2h_q+1}\big(\frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}}+\frac{Pa^{{\mathbf h}}}{a^{h_q}} R_n(\ell) + R^*_n(\ell) \frac{a^{{\mathbf h}}P^*}{a^{h_q}}\big)$ implies

\frac{\lambda_{p-r+\ell}(\boldsymbol{\mathcal W}(a(n)2^j,R_n(\ell)))}{a^{2h_\ell+1}} \to \xi_\ell(2^j),\quad \ell=1\ldots,r.

In particular, for $\ell_*\notin \mathcal I_0$,

\sup_{\|R\|\leq \varepsilon}\lambda_{p-r+\ell_*}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} R + R^*\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)

=:\lambda_{p-r+\ell_*}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} R_n(\ell_*) + R_n(\ell_*)^*\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big) \stackrel
\bbP \to\begin{cases}
\infty, & \ell_*\in \mathcal{I}_+;\\
0, & \ell_*\in \mathcal{I}_- \cup \{1,\ldots,p-r\},
\end{cases}

where the statement for $\ell_*\in\{1,\ldots,p-r\}$ follows from Lemma <ref> since $\lambda_{p-r} \big(\boldsymbol{\mathcal W}(a(n)2^j,R_n(\ell))\leq\sup_{\|R\|\leq \varepsilon}\lambda_{p-r} \big(\boldsymbol{\mathcal W}(a(n)2^j,R)\big) =\frac{O_{\bbP}(1)}{a^{2h_q+1}}$. Therefore, claim (<ref>) is established.

Thus, for $R$ with small norm $\|R\|$, an adaptation of Proposition 3 in Abry and Didier <cit.> to rectangular matrices implies

f_{n,q,2}(R) - f_{n,q,2}(0) = \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} R + R^* \frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)

- \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} +\frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}}  \Big)
= \sum^{r}_{i_1 = 1}\sum^{p}_{i_2 = 1}\frac{\partial}{\partial r_{i_1,i_2}}f_{n,q,2}(\breve{R}) \hspace{1mm} \pi_{i_1,i_2}(R)

for some matrix $\breve{R} \in {\mathcal M}(r,p,\bbR)$ lying in a segment connecting $R$ and ${\mathbf 0}$ across ${\mathcal M}(r,p,\bbR)$. Consider the matrix of derivatives

\Big\{\frac{\partial}{\partial r_{i_1,i_2}}f_{n,q,2}(\breve{R})\Big\}_{i_1= 1,\hdots,r; \hspace{1mm}i_2 = 1,\hdots,p}

= \Big\{ \lambda^{-1}_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \breve{R} + \breve{R}^* \frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)

\begin{equation}\label{e:df/dri1,i2}
\times \frac{\partial}{\partial r_{i_1,i_2}}\lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} +\frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \breve{R} + \breve{R}^*  \frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)\Big\}_{i_1= 1,\hdots,r; \hspace{1mm}i_2 = 1,\hdots,p}.
\end{equation}

In (<ref>), the differential of the eigenvalue
\lambda_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} +\frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \breve{R} + \breve{R}^* \frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)
is given by expression (<ref>) with

\begin{equation}\label{e:MVT_arg_fnq2}
\frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} +\frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \breve{R} + \breve{R}^* \frac{a^{{\mathbf h}}P^*}{a^{h_q}}
\end{equation}

in place of $M$ and ${\mathbf u}_{p-r+q}(n,\breve{R})$ denoting a unit eigenvector of (<ref>) associated with its $(p-r+q)$–th eigenvalue. Moreover,

\frac{\partial}{\partial r_{i_1,i_2}}\lambda_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} \mathbf{B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} +\frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \breve{R} + \breve{R}^* \frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)

= {\mathbf u}^*_{p-r+q}(n, \breve{R})\frac{\partial}{\partial r_{i_1,i_2}}\Big[\frac{Pa^{{\mathbf h}} \mathbf{B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} +\frac{{\mathbf W}_Z(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \breve{R} + \breve{R}^* \frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big]{\mathbf u}_{p-r+q}(n,\breve{R})

= {\mathbf u}^*_{p-r+q}(n,\breve{R})\Big( \frac{Pa^{{\mathbf h}}}{a^{h_q}} {\mathbf 1}_{i_1,i_2}  + {\mathbf 1}_{i_2,i_1} \frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big){\mathbf u}_{p-r+q}(n,\breve{R}), \quad \quad  i_1= 1,\hdots,r, \quad i_2= 1,\hdots,p.

\begin{equation}\label{e:sqrt(n)(fnq2(R)-fnq2(0))}
\sqrt{\frac{n}{a}}\Big(f_{n,q,2}(R) - f_{n,q,2}(\mathbf{0})\Big)
\end{equation}

= \lambda^{-1}_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} \mathbf{B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_{Z}(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \breve{R} + \breve{R}^* \frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)

\sum^{r}_{i_1 = 1}\sum^{p}_{i_2 = 1}\Big\{ {\mathbf u}^*_{p-r+q}(n,\breve R)\Big( \frac{Pa^{{\mathbf h}}}{a^{h_q}} {\mathbf 1}_{i_1,i_2} + {\mathbf 1}_{i_2,i_1}\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big){\mathbf u}_{p-r+q}(n,\breve R) \Big\}_{i_1 = 1,\hdots,r; \hspace{1mm}i_2 = 1,\hdots,p}
\hspace{1mm}\sqrt{\frac{n}{a}} \pi_{i_1,i_2} (R).

Note that, by condition (<ref>),

\begin{equation}\label{e:|WXZ|=<delta/2}
\frac{1}{a^{h_q+1/2}}\hspace{0.5mm}\|a^{- {\mathbf h}- (1/2)I}{\mathbf W}_{X,Z}(a2^j)\| = o_{\bbP}(1),
\end{equation}

namely, $\|\frac{a^{- {\mathbf h}- (1/2)I}{\mathbf W}_{X,Z}(a2^j)}{a^{h_q+1/2}} \|$ is arbitrarily small, in probability, for large enough $n$. Thus, by (<ref>) and (<ref>),

\begin{equation}\label{e:sqrt(n/a)(fnq2(o(1))-fnq0(0))}
\sqrt{\frac{n}{a}}\Big(f_{n,q,2}\Big( \frac{a^{- {\mathbf h}- (1/2)I}{\mathbf W}_{X,Z}(a2^j)}{a^{h_q + 1/2}} \Big) - f_{n,q,2}(\mathbf{0})\Big)
\end{equation}

= \lambda^{-1}_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} \mathbf{B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_{Z}(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \breve{R}_n + \breve{R}^*_n \frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big)

\times \sum^{r}_{i_1 = 1}\sum^{p}_{i_2 = 1}{\mathbf u}^*_{p-r+q}(n)\Big( \frac{Pa^{{\mathbf h}}}{a^{h_q}} {\mathbf 1}_{i_1,i_2} + {\mathbf 1}_{i_2,i_1}\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big){\mathbf u}_{p-r+q}(n)
\times \sqrt{\frac{n}{a}} \pi_{i_1,i_2} \Big( \frac{a^{- {\mathbf h}- (1/2)I}{\mathbf W}_{X,Z}(a2^j)}{a^{h_q + 1/2}} \Big).

Note that $\breve R_n$ lies in a segment connecting $\mathbf 0$ and  the matrix $ \frac{a^{- {\mathbf h}- (1/2)I}{\mathbf W}_{X,Z}(a2^j)}{a^{h_q + 1/2}}$, whose norm tends to $0$ in probability due to condition (<ref>).  Hence, $\breve R_n= o_{\bbP}(1)$. On the other hand,
Corollary <ref> implies

\begin{equation}\label{e:max_(up-r+q,p_ell(n))a^(h_ll-hq)}
\lambda_{p-r+q}\Big(\frac{Pa^{{\mathbf h}} \mathbf{B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_{Z}(a2^j)}{a^{2h_q+1}} + \frac{Pa^{{\mathbf h}}}{a^{h_q}} \breve{R}_n + \breve{R}^*_n \frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big) \rightarrow \xi_{q}(2^j), \quad n \rightarrow \infty.
\end{equation}

As a consequence of condition (<ref>) and of (<ref>),
\max_{i_1 = 1,\hdots,r; \hspace{0.5mm}i_2=1,\hdots,p}|\pi_{i_1,i_2}( a^{-\mathbf{h}-(1/2)I}{\mathbf W}_{X,Z}(a2^j) )| = O_{\bbP}(1).
Therefore, by expressions (<ref>), (<ref>) and by condition (<ref>), the right-hand side of (<ref>) is bounded, in absolute value, by

O_{\bbP}(1) \Big|\sum^{r}_{i_1 = 1}\sum^{p}_{i_2 = 1} {\mathbf u}^*_{p-r+q}(n)\Big( \frac{Pa^{{\mathbf h}}}{a^{h_q}} {\mathbf 1}_{i_1,i_2} + {\mathbf 1}_{i_2,i_1}\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big){\mathbf u}_{p-r+q}(n) \hspace{1mm} \sqrt{\frac{n}{a}} \pi_{i_1,i_2} \Big( \frac{a^{- {\mathbf h}- (1/2)I}{\mathbf W}_{X,Z}(a2^j)}{a^{h_q + 1/2}} \Big)\Big|

\leq O_{\bbP}(1) \sqrt{\frac{n}{a}} \Big( \frac{1}{a^{h_q + 1/2}} \Big)\sum^{r}_{i_1 = 1}\sum^{p}_{i_2 = 1} \Big|{\mathbf u}^*_{p-r+q}(n)\Big( \frac{Pa^{{\mathbf h}}}{a^{h_q}} {\mathbf 1}_{i_1,i_2} + {\mathbf 1}_{i_2,i_1}\frac{a^{{\mathbf h}}P^*}{a^{h_q}}\Big){\mathbf u}_{p-r+q}(n)\Big|

= O_{\bbP}(1) \sqrt{\frac{n}{a}} \Big( \frac{1}{a^{h_q + 1/2}} \Big)\sum^{r}_{i_1 = 1}\Big| \langle{\mathbf p}_{i_1},{\mathbf u}_{p-r+q}(n)\rangle a^{h_{i_1}-h_q} \Big| \hspace{1mm}\sum^{p}_{i_2 = 1}  | {\mathbf u}_{p-r+q}(n)_{i_2} |

\begin{equation}\label{e:sqrt(n/a)(fnq2(rescaled_WXZ)-fnq2(0))_bound}
\leq O_{\bbP}(1) \sqrt{\frac{n}{a}} \Big( \frac{1}{a^{h_q + 1/2}} \Big) \hspace{1mm} \sqrt{p}
\leq O_{\bbP}(1) \sqrt{\frac{p}{n/a}} \Big(\frac{n}{a}\Big) \frac{1}{a^{h_q + 1/2}}
\leq O_{\bbP}(1) \frac{n}{a^{2\varpi + 3/2}} \stackrel{\bbP}\rightarrow 0.
\end{equation}

In (<ref>), the last inequality results from the bound $(h_q - h_1) + h_1 \geq 2 \varpi$. In other words, we can write

\sqrt{\frac{n}{a}}\Big\{\log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B} (2^j)a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_{Z}(a2^j)}{a^{2h_q+1}}+
\frac{Pa^{{\mathbf h}}{\mathbf W}_{X,Z}(a2^j)}{a^{2h_q+1/2}}+\frac{{\mathbf W}^*_{X,Z}(a2^j)a^{{\mathbf h}}P^*}{a^{2h_q+1/2}}\Big)

- \log \lambda_{p-r+q}\Big( \frac{Pa^{{\mathbf h}} {\mathbf B}_a(2^j) a^{{\mathbf h}}P^*}{a^{2h_q}} + \frac{{\mathbf W}_{Z}(a2^j)}{a^{2h_q+1}}\Big)\Big\}

= \sqrt{\frac{n}{a}}\Big\{f_{n,q,2}\Big( \frac{a^{- {\mathbf h}- (1/2)I}{\mathbf W}_{X,Z}(a2^j)}{a^{h_q + 1/2}} \Big) - f_{n,q,2}(0)\Big)\Big\} = o_{\bbP}(1).

Hence, (<ref>) holds. $\Box$

§ AUXILIARY STATEMENTS

The proofs of Section <ref> depend on a series of auxiliary statements that are established in this section. As in the previous section, whenever convenient we simply write $p = p(n)$, $a = a(n)$ or $P = P(n)$. Also, for notational simplicity, as in (<ref>) we write $P(n)P_H \equiv P(n) \equiv P$. Several lemmas in this section are stated in terms of matrices $\widetilde{\W}(a2^j)$, which are slightly generalized versions of $\W(a2^j)$ (such matrices appear throughout the proof of Theorem <ref>).

The following lemma is used in Theorem <ref>. It is also applied in Lemma <ref>, which in turn is needed to establish Theorem <ref>.

Fix $j\in \bbN$, and $q\in\{1,\ldots,r\}$. Let $\mathbf B_n$ denote either $\widehat{\mathbf B}_a(2^j)$ or ${\mathbf B}_a(2^j)$, and let

$$\widetilde{\W}(a2^j) =Pa^{{\mathbf h + \frac{1}{2}I}}\hspace{0.25mm}\mathbf B_n\hspace{0.25mm} a^{{\mathbf h + \frac{1}{2}I}} P^* + M_n,$$

\begin{equation}\label{e:def_Mn}
\frac{Pa^{{\mathbf h}}O_{\bbP}(1)}{a^{2h_q+1/2}}+\frac{O_{\bbP}(1)^*a^{{\mathbf h}}P^*}{a^{2h_q+1/2}}\big).
\end{equation}

In (<ref>), the terms $O_\bbP(1)$ denote any matrix bounded in probability of appropriate dimension.
Then, under $(W1 - W4)$ and $(A1 - A5)$, as $n \rightarrow \infty$,

\begin{equation}\label{e:lambdap-r=O_P(1)}
\lambda_{p-r}(\widetilde{\W}(a2^j)) = O_\bbP(1)
\end{equation}

\begin{equation}\label{e:log_lambdaq_consistency}
\frac{\lambda_{p-r+1}(\widetilde{\W}(a2^j))}{a(n)^{2h_1 + 1}} = \Omega_{\bbP}(1).
\end{equation}

In particular, statements (<ref>), (<ref>) hold for the matrix $\mathbf{W}(a2^j)$, and $\lambda_{p-r}(\mathbf{W}(a2^j)) = O_\bbP(1)$, $\lambda_{p-r+1}(\mathbf{W}(a 2^j))\stackrel \bbP \to \infty$.

To establish (<ref>), just note that, by (<ref>),

0 \leq \lambda_{p-r}(\widetilde{\W}(a2^j)) \leq \sup_{{\mathbf u} \in \{\mathbf{p}_{1},\hdots,\mathbf{p}_{r}\}^{\perp} \cap S^{p-1}}{\mathbf u}^* \widetilde{\W}(a2^j){\mathbf u}
= \sup_{{\mathbf u} \in \{\mathbf{p}_{1},\hdots,\mathbf{p}_{r}\}^{\perp}\cap S^{p-1}}{\mathbf u}^* O_\bbP(1){\mathbf u} = O_\bbP(1).

To establish (<ref>), first rewrite

\begin{equation}\label{e:lambda_p-q+1(W(a2^j))}
\lambda_{p-r+1}(\widetilde{\W}(a2^j))=\sup_{{\mathcal U}_{r}}\inf_{{\mathbf u} \in {\mathcal U}_{r} \cap S^{p-1}_{\bbC}} {\mathbf u}^*\widetilde{\W}(a2^j){\mathbf u}
\geq \inf_{{\mathbf u} \in \textnormal{span}\{ \mathbf p_{1},\ldots,\mathbf p_r \}}{\mathbf u}^*\widetilde{\W}(a2^j){\mathbf u}=: \widetilde{{\mathbf u}}^*(n)\widetilde{\W}(a2^j)\widetilde{{\mathbf u}}(n).
\end{equation}

In (<ref>),

\begin{equation}\label{e:w(n)_in_span(p1,...,pr)}
\widetilde{{\mathbf u}} \equiv \widetilde{{\mathbf u}}(n) \in \textnormal{span}\{ \mathbf p_{1},\ldots,\mathbf p_r \}
\end{equation}

is an appropriately chosen unit vector, and, like the vectors ${\mathbf p}_{\ell}={\mathbf p}_{\ell}(n)$, $\ell = 1,\hdots,r$, it is a function of $n$. Now,

\begin{equation}\label{e:W=PWXP*+WZ+PWXZ+WXZ*P*_2}
a^{-1}\widetilde{\W}(a2^j)= Pa^{\mathbf h} \mathbf B_n a^{\mathbf h}P^*+ O_\bbP(1) + P a^{\mathbf h} O_\bbP(1) + O^*_\bbP(1)a^{\mathbf h}P^*.
\end{equation}

In view of (<ref>), there is some constant $C>0$ such that $\|P^*\widetilde{{\mathbf u}}\|>C$ with probability tending to 1. In turn, this implies that, for some $C'>0$, $\|a^{\mathbf h}P^*\widetilde{{\mathbf u}}\|\geq C' a^{h_1}$, again with probability tending to $1$. Therefore, by (<ref>), as well as by conditions (<ref>) and (<ref>), there is some constant $C'' >0$ such that

a^{-1}\widetilde{{\mathbf u}}^*\widetilde{\W}(a2^j)\widetilde{{\mathbf u}}

=\|a^{\mathbf h}P^*\widetilde{{\mathbf u}}\|^2 \Big\{ \frac{\widetilde{{\mathbf u}}^*Pa^{\mathbf h}}{\|\widetilde{{\mathbf u}}^*Pa^{\mathbf h}\|} \mathbf B_n \frac{a^{\mathbf h}P^*\widetilde{{\mathbf u}}}{\|a^{\mathbf h}P^*\widetilde{{\mathbf u}}\|}
+ \frac{\widetilde{{\mathbf u}}^* O_\bbP(1)\widetilde{{\mathbf u}}}{\|a^{\mathbf h}P^*\widetilde{{\mathbf u}}\|^2} + \frac{2}{\|\widetilde{{\mathbf u}}^*P a^{\mathbf h}\|} \widetilde{{\mathbf u}}^*O_\bbP(1)\frac{a^{\mathbf h}P^*\widetilde{{\mathbf u}}}{\|a^{\mathbf h}P^*\widetilde{{\mathbf u}}\|} \Big\}

\geq  \|a^{\mathbf h}P^*\widetilde{{\mathbf u}}\|^2 \big\{\lambda_1(\mathbf B_n)+ o_\bbP(1) \big\} \geq C'' a^{2h_1}

with probability tending to 1, since $\plim_{n\to\infty} \mathbf B_n=\mathbf B(2^j) \in\mathcal{S}_{>0}(\bbR)$ almost surely due to assumpmtion $(A3)$ (see  (<ref>)). Hence, (<ref>) is established.

The following lemma is used in the proof of Theorem <ref>, Proposition <ref> and Lemma <ref>.

Fix $q \in\{1,\hdots,r\}$ and consider the index sets $\mathcal{I}_-,\mathcal{I}_0,$ and $\mathcal{I}_+$ as in (<ref>). Suppose $\mathcal{I}_+ \neq \emptyset$, and let $\widetilde{\mathbf W}(a2^j)$ be given as in Lemma <ref>. Also, let ${\boldsymbol {\mathfrak u}}_{p-r+q}(n)$ be a unit eigenvector associated with the $(p-r+q)$–th eigenvalue of $\widetilde{\mathbf W}(a2^j)$. Then, for every $\varepsilon>0$,

\begin{equation}\label{e:supR'max(angles*powerlaws)=O_P(1)}
\max_{\ell\in \mathcal{I}_+}\{ |\langle {\mathbf p}_{  \ell}(n),{\boldsymbol {\mathfrak u}}_{p-r+q}(n)\rangle| \hspace{0.5mm}a(n)^{h_{\ell}-h_q}\} = O_{\bbP}(1), \quad n \rightarrow \infty.
\end{equation}

In particular, for unit eigenvectors $\mathbf u_{p-r+q}(n)$ of $\W(a2^j)$,

\begin{equation}\label{e:max|<up-r+q(n),pi(n)>|a(n)^(hi-hq)=OP(1)}
\max_{\ell\in\mathcal{I_+}} \big\{ |\langle{\mathbf p}_{  \ell}(n), \mathbf u_{p-r+q}(n) \rangle| \hspace{0.5mm}a(n)^{h_{\ell}-h_q} \big\} = O_{\bbP}(1), \quad n \rightarrow \infty.
\end{equation}

First note that

\begin{equation}\label{e:lambdaq/a^(2h1+1)=O(1)}
\frac{\lambda_{p-r+q}(\widetilde{\mathbf W}(a2^j))}{a(n)^{2h_q + 1}} = \inf_{{\mathcal U}_{p-r+q}}\sup_{{\mathbf u} \in {\mathcal U}_{p-r+q} \cap S^{p-1}} {\mathbf u}^*\frac{\widetilde{\mathbf W}(a2^j)}{a(n)^{2h_q + 1}}{\mathbf u} \leq \sup_{{\mathbf u} \in {\textnormal{span}\{\mathbf p_{q+1},\ldots,\mathbf{p}_r\}}^{\perp} \cap S^{p-1}} {\mathbf u}^*\frac{\widetilde{\mathbf W}(a2^j)}{a(n)^{2h_q + 1}}{\mathbf u}
\leq \sup_{{\mathbf u} \in {\textnormal{span}\{\mathbf p_\ell,\ell\in \mathcal{I}_+\}}^{\perp} \cap S^{p-1}} {\mathbf u}^*\frac{\widetilde{\mathbf W}(a2^j)}{a(n)^{2h_q + 1}}{\mathbf u} = O_{\bbP}(1).
\end{equation}

By way of contradiction,
suppose that

\mathcal{I}_*:=\Big\{\ell\in\{1,\ldots,r\}: \limsup_{n\to\infty} |\langle {\mathbf p}_\ell,{\boldsymbol {\mathfrak u}}_{p-r+q}(n)\rangle | \hspace{1mm} a(n)^{h_\ell-h_q}= \infty \Big\}\subseteq \mathcal{I}_+

is nonempty with non-vanishing probability. Equivalently, we can pick some $\ell_1\in \mathcal{I}_*$, and a subsequence $n_1$ such that $ |\langle {\mathbf p}_{\ell_1},{\boldsymbol {\mathfrak u}}_{p-r+q}(n_1)\rangle | \hspace{1mm} a(n_1)^{h_{\ell_1}-h_q}\to\infty$.  If there exists some $\ell_2\in\mathcal{I}_*$ such that
\limsup_{n_1\to\infty}\frac{|\langle{\mathbf p}_{\ell_2},{\boldsymbol {\mathfrak u}}_{p-r+q}(n_1)\rangle| a(n_1)^{h_{\ell_2}}}{|\langle{\mathbf p}_{\ell_1},{\boldsymbol {\mathfrak u}}_{p-r+q}(n_1)\rangle| a(n_1)^{h_{\ell_1}}} = \infty,
then we can take a subsequence $n_2$ of $n_1$ such that
\lim_{n_2\to\infty}\frac{|\langle{\mathbf p}_{\ell_2},{\boldsymbol {\mathfrak u}}_{p-r+q}(n_2)\rangle| a(n_2)^{h_{\ell_2}}}{|\langle{\mathbf p}_{\ell_1},{\boldsymbol {\mathfrak u}}_{p-r+q}(n_2)\rangle| a(n_2)^{h_{\ell_1}}} = \infty
(n.b.: $\limsup$ was replaced with $\lim$). Repeating this process, if there exists some $\ell_{i+1}\in\mathcal{I}_*$ such that

\begin{equation}\label{e:ith_subseq}
\limsup_{n_i\to\infty}\frac{|\langle{\mathbf p}_{\ell_{i+1}},{\boldsymbol {\mathfrak u}}_{p-r+q}(n_i)\rangle| a(n_i)^{h_{\ell_{i+1}}}}{|\langle{\mathbf p}_{\ell_i},{\boldsymbol {\mathfrak u}}_{p-r+q}(n_i)\rangle| a(n_i)^{h_{\ell_i}}} = \infty,
\end{equation}

take a refinement $n_{i+1}$ of the previous subsequence $n_{i}$ so that (<ref>) holds as an (infinite) limit, as opposed to a $\limsup$. Let $r_*\leq |\mathcal{I}_*|$ be the maximal index from this construction such that (<ref>) holds as a limit. By definition, such limit holds along the sequence $n_{r_*}$ (denoted in the sequel by $n'$, for notational simplicity).   Then, clearly, if $r_*>1$,
\max_{i=1,\ldots,r_*-1} \lim_{n'\to\infty}\frac{|\langle{\mathbf p}_{\ell_{i}},{\boldsymbol {\mathfrak u}}_{p-r+q}(n')\rangle| a(n')^{h_{\ell_{i}}}}{|\langle{\mathbf p}_{\ell_{r_*}},{\boldsymbol {\mathfrak u}}_{p-r+q}(n')\rangle| a(n')^{h_{\ell_{r_*}}}} =0.
In particular, for  $r_*\geq 1$, $\frac{|\langle{\mathbf p}_{\ell},{\boldsymbol {\mathfrak u}}_{p-r+q}(n')\rangle| a(n')^{h_{\ell}}}{|\langle{\mathbf p}_{\ell_{r_*}},{\boldsymbol {\mathfrak u}}_{p-r+q}(n')\rangle| a(n')^{h_{\ell_{r_*}}}} =O_{\bbP}(1)$, $\ell\in \mathcal{I}_+$. Hence, refining the subsequence $n'$ further if necessary we may assume that, for each $\ell = 1,\hdots,r$,

\begin{equation}\label{e:lim_ratio_angle*power_law=y_ell}
\lim_{n'\to\infty}\frac{\langle{\mathbf p}_{\ell},{\boldsymbol {\mathfrak u}}_{p-r+q}(n')\rangle a(n')^{h_{\ell}}}{\langle{\mathbf p}_{\ell_{r_*}},{\boldsymbol {\mathfrak u}}_{p-r+q}(n')\rangle a(n')^{h_{\ell_{r_*}}}}=:y_\ell \in \bbR.
\end{equation}

Note that, in (<ref>), $y_\ell \neq 0 \Rightarrow \ell \in \mathcal I_*\subseteq \mathcal{I}_+$. Moreover, there must be at least one such $\ell$. In other words, we can write $\mathbf{y}^* = (y_1,\ldots,y_r)\neq\mathbf{0}$, where, with non-vanishing probability,

\begin{equation}\label{e:a^h_P*_u/fastest->y}
\frac{a^{{\mathbf h}}P^*{\boldsymbol {\mathfrak u}}_{p-r+q}(n')}{|\langle{\mathbf p}_{\ell_{r_*}},{\boldsymbol {\mathfrak u}}_{p-r+q}(n')\rangle| a(n')^{h_{\ell_{r_*}}}} \to \mathbf{y},\qquad n'\to\infty.%
\end{equation}

So, consider the (divergent) sequence $\vartheta(n'):=|\langle{\mathbf p}_{\ell_{r_*}},{\boldsymbol {\mathfrak u}}_{p-r+q}(n')\rangle| a(n')^{h_{\ell_{r_*}}-h_q}$. Then, by relation (<ref>),

\frac{1}{\vartheta(n')^2}\frac{\lambda_{p-r+q}(\widetilde{\mathbf W}(a2^j))}{a(n')^{2h_q + 1}} = {\mathbf u}^*_{p-r+q}(n')\frac{\widetilde{\mathbf W}(a2^j)}{\vartheta(n')^2a(n')^{2h_q + 1}}{\boldsymbol {\mathfrak u}}_{p-r+q}(n')

\begin{equation}\label{e:lambda/power_law/fastest}
={\mathbf u}^*_{p-r+q}(n',R_{n'})\frac{Pa^{{\mathbf h}}{\mathbf B}_{n'} a^{{\mathbf h}}P^*}{\vartheta(n')^2a^{2h_q}}{\boldsymbol {\mathfrak u}}_{p-r+q}(n',R_{n'}) \hspace{1mm}\{1 + o_{\bbP}(1)\}
\to \mathbf{y}^* \mathbf{B}(2^j) \mathbf{y}>0,
\end{equation}

due to assumption $(A3)$. The limit (<ref>) implies that $\frac{\lambda_{p-r+q}(\widetilde{\mathbf W}(a2^j))}{a(n)^{2h_q + 1}} $ is unbounded with positive probability, which contradicts (<ref>). Hence, (<ref>) is established. $\Box$

Under the convention $P(n)\equiv P(n)P_H$ (see (<ref>)), recall that

\W(a(n)2^j)=\frac{P(n)a(n)^{{\mathbf h}} \widehat{{\mathbf  B}}_a(2^j) a(n)^{{\mathbf h}}P^*(n)}{a(n)^{2h_q}} + \frac{{\mathbf W}_Z(a(n)2^j)}{a(n)^{2h_q+1}}

+ \frac{P(n)a(n)^{{\mathbf h}}}{a(n)^{h_q}} \big[ a(n)^{-{\mathbf h}} P^{-1}_H  \W_{X,Z}(a(n)2^j) \big] + \big[ a(n)^{-{\mathbf h}} P^{-1}_H  \W_{X,Z}(a(n)2^j) \big]^*\frac{a(n)^{{\mathbf h}}P^*(n)}{a(n)^{h_q}}.

Thus, as a consequence of expression (<ref>) in Lemma <ref>,

\begin{equation}\label{e:lambda_p-r+q/scaling=OP(1)}
\frac{\lambda_{p-r+q}({\mathbf W}(a(n)2^j))}{a(n)^{2h_q+1}}=O_{\bbP}(1), \quad q=1,\ldots,r.
\end{equation}

Relation (<ref>) is invoked in the proof of Theorem <ref>.

The following lemma is used in the proofs of Lemma <ref> and Proposition <ref>.

For each $n$, let $\p_\ell(n)$, $\ell=r+1,\ldots,p$, be an orthonormal basis for the nullspace of $P(n)$. Let $\varpi$ be as in (<ref>) and fix any $q \in\{ 1,\hdots,r\}$, and let ${\boldsymbol {\mathfrak u}}_{p-r+q}(n)$ denote a unit eigenvector associated with the $(p-r+q)$–th eigenvalue of $\widetilde{ \mathbf W}(a2^j)$ as given in Lemma <ref>. Then,

\begin{equation}\label{e:sum<pi(n),up-r+q(n)>2-o(a-varpi)_first}
\sum^{p}_{i=r+1}\langle {\mathbf p}_{  i}(n), {\boldsymbol {\mathfrak u}}_{p-r+q}(n)\rangle^2 = o_\bbP\big(a(n)^{-2\varpi}\big).
\end{equation}

In particular, (<ref>) holds for unit eigenvectors $ \mathbf u_{p-r+q}(n)$ associated with the matrix ${\mathbf W}(a2^j)$.

By way of contradiction, suppose that, for some $\varepsilon_0 > 0$ and some subsequence (still denoted $n \in \bbN$, for notational simplicity),

\begin{equation}\label{e:sum_r+1...p_<pi,u>^2>0}
a(n)^{2\varpi}\sum^{p}_{i=r+1} \langle {\mathbf p}_{  i}(n), {\boldsymbol {\mathfrak u}}_{p-r+q}(n)\rangle^2  \geq \varepsilon_0 > 0
\end{equation}

with non-vanishing probability. Consider the projection $\mathbf v(n)$ of ${\boldsymbol{\mathfrak u}}_{p-r+q}(n)$ onto the nullspace of $P(n)$, i.e.,

\begin{equation}\label{e:v*(n)}
\mathbf v (n) =  \sum^{p}_{i=r+1} \langle \p_{i}(n),{\boldsymbol {\mathfrak u}}_{p-r+q}(n)\rangle \hspace{0.5mm} \p_{i}(n).
\end{equation}

Clearly, $ \mathbf v^*(n)\{P(n)O_\bbP(1)+ O_\bbP(1)^* P^*(n)\}\mathbf v(n) = 0$. In view of condition (<ref>), with non-vanishing probability,

\begin{equation}\label{e:v*(n)W(a2j)v(n)=OP(1)}
\mathbf v^*(n)\widetilde {\mathbf W}(a2^j)\mathbf v(n) = O_{\bbP}(1).
\end{equation}

However, let $O(n) = ( \mathbf {\boldsymbol {\mathfrak u}}_{1}(n)\hdots \mathbf {\boldsymbol {\mathfrak u}}_{p}(n) )$ be a matrix of eigenvectors of $\widetilde {\mathbf W}(a2^j)$. Then, for any $q \in\{ 1,\hdots,r\}$,

\mathbf v^*(n)\widetilde{\mathbf W}(a2^j)\mathbf v(n)= \mathbf v^*(n)O(n) \textnormal{diag}\Big(\lambda_1(\widetilde {\mathbf W}(a2^j)),\hdots,\lambda_{p}(\widetilde {\mathbf W}(a2^j))\Big)O^*(n)\mathbf v(n)

= \sum^{p}_{i=1} \langle{\boldsymbol {\mathfrak u}}_{i}(n), \mathbf v(n)\rangle^2 \lambda_i(\widetilde {\mathbf W}(a2^j)) \geq \langle {\boldsymbol {\mathfrak u}}_{p-r+q}(n),\mathbf  v(n)\rangle^2 \lambda_{p-r+q}(\widetilde {\mathbf W}(a2^j))

\begin{equation}\label{e:v*(n)W(a2j)v(n)_lower_bound}
= \Big( \sum^{p}_{i=r+1} \langle \p_{i}(n), {\boldsymbol {\mathfrak u}}_{p-r+q}(n)\rangle^2 \Big)^2 \lambda_{p-r+q}(\widetilde{\mathbf W}(a2^j)) \geq \varepsilon^2_0 \hspace{0.5mm} a(n)^{-2\varpi}\lambda_{p-r+q}(\widetilde {\mathbf W}(a2^j)),
\end{equation}

where the last inequality follows from (<ref>). As a consequence of relation (<ref>) in Lemma <ref>, and of the fact that $2\varpi\leq 2h_1 <2h_q+1$, the right-hand side of (<ref>) blows up in probability, which contradicts (<ref>). This shows (<ref>).

The following lemma is used in the proofs of Theorem <ref>, Proposition <ref>, and Lemma <ref>. It contains statements about the properties of limits of convergent subsequences of matrices of inner products between coordinates and random eigenvectors (on sufficient conditions for the convergence of the main sequence, see Proposition <ref>).

For $\ell = 1,\hdots,p$, let ${\boldsymbol {\mathfrak u}}_{\ell}(n)$ be an eigenvector associated with the $\ell$–th eigenvalue of the random matrix $\widetilde {\mathbf W}(a(n)2^j)$ as defined in  Lemma <ref>. Define

\begin{equation}\label{e:U_r(n)}
\mathbf{U}_r(n) =[\boldsymbol{\mathfrak u}_{p-r+1}(n),\ldots,\boldsymbol{\mathfrak u}_{p}(n)].
\end{equation}

\begin{equation}\label{e:Gamma_n=P*(n)U_r(n)}
\Gamma_n = P^*(n) \mathbf U_r(n)
\end{equation}

be the matrix containing the inner products between the coordinate vectors in $P(n)$ and the eigenvectors in (<ref>). Consider any subsequence $\{\Gamma_{n'}\}_{n'}$ of $\{\Gamma_n\}_{n \in \bbN}$ along which the limit

\begin{equation}\label{e:def_Gamma}
\plim_{n'\to\infty}\Gamma_{n'} =:\Gamma \in M(r,\bbR)
\end{equation}

exists, where $\Gamma$ is a constant matrix a.s.  Then,

\begin{equation}\label{e:Gamma_in_GL(r,R)}
\Gamma \in GL(r,\bbR) \quad \textnormal{a.s.}
\end{equation}

Also, for each $q\in\{1,\ldots,r\}$, if $\mathcal{I}_+=\{\ell:h_\ell>h_q\} \neq \emptyset$ (cf. (<ref>)), then the submatrix

\begin{equation}\label{e:def_Gamma_plusblock}
\Gamma_+:=(\Gamma_{\ell,\ell'})_{\ell,\ell'\in\mathcal{I}_+}
\end{equation}

of $\Gamma$ is nonsingular.

Let $q\in\{1,\ldots,r\}$. For notational simplicity, we write $n = n'$ to index the subsequence. So, for each $n$, let $\mathbf{p}_{r+1}(n),\ldots,\mathbf{p}_{p}(n)$ be an orthonormal basis for the nullspace of ${P(n)}$.  As a consequence of Lemma <ref>,

\begin{equation}\label{e:sum<pi(n),up-r+q(n)>2->0}
\sum^{p}_{i=r+1} \langle {\mathbf p}_{  i}(n), {\boldsymbol {\mathfrak u}}_{p-r+q}(n)\rangle^2 \stackrel{\bbP}\rightarrow 0, \quad n \rightarrow \infty.
\end{equation}

Also, for all $n \in \bbN$ and some $Q(n)\in M(p,r,\bbR)$, $R(n)\in GL(r,\bbR)$, let

\begin{equation}\label{e:P(n)=Q(n)R(n)}
\end{equation}

be the $QR$ decomposition of $P(n)$ (e.g., Horn and Johnson <cit.>, Theorem 2.1.14, (a)). Write $\widetilde \p_1(n),\ldots,\widetilde \p_r(n)$ for the (orthonormal) columns of $Q(n)$. Recall that, for a matrix $M=(m_{i,i'})\in M(r,\bbR)$,

\begin{equation}\label{e:tr(M*M)=sum_singval^2}
\sum^{r}_{i,i'=1} m_{i,i'}^2= \textnormal{tr}(M^* M) = \sum_{i=1}^r\sigma_i^2(M),
\end{equation}

where $\sigma_i(M)$, $i = 1,\hdots,r$, are the singular values of $M$. Before showing (<ref>), we establish that, for $\mathbf{U}_r(n)$ as in (<ref>), there exists some constant $C>0$ such that

\begin{equation}\label{e:sigma2_1(Q*(n)Ur(n))>=C>0}
\sigma_{1}^2\big(Q^*(n)\mathbf{U}_r(n) \big)\geq C>0
\end{equation}

for all large $n$ with probability going to 1. In fact, (<ref>) implies that, as $n \rightarrow \infty$,
$\sum_{i=1}^r \langle \widetilde \p_i(n), {\mathbf u}_{p-r+q}(n)\rangle^2  \stackrel \bbP \rightarrow 1$, $q=1,\ldots,r$. Then, by relation (<ref>) and again as $n \rightarrow \infty$,

\begin{equation}\label{e:sum_of_singular_values_tends_to_r}
\sum^{r}_{s=1}\sigma_{s}^2\Big(Q^*(n)\mathbf{U}_r(n)\Big) = \sum_{i,\ell=1}^r \langle \widetilde{\p}_i(n),{\boldsymbol {\mathfrak u}}_{p-r+\ell}(n)\rangle^2 \stackrel \bbP \rightarrow r.
\end{equation}

However, the largest singular value $\sigma_{r}^2\big(Q^*(n)\mathbf{U}_r(n)\big)$ satisfies $\sigma_{r}^2\big(Q^*(n)\mathbf{U}_r(n)\big)\leq \|Q^*(n)\|^2\|\mathbf{U}_r(n)\|^2 = 1$. Hence, by (<ref>), $\sigma_{1}^2\big(Q^*(n)\mathbf{U}_r(n)\big)\stackrel{\bbP}\to 1$. Thus, (<ref>) holds, as claimed.

We are now in a position to show (<ref>). Suppose $\Gamma\mathbf{x}=0$ for some $\mathbf{x}\in\bbR^r$. Note that

\begin{equation*}
\|R^*(n)Q^*(n) \mathbf{U}_r(n) \mathbf{x}\| \geq \sigma_1(R^*(n)) \hspace{0.5mm} \sigma_1\big(Q^*(n)\mathbf{U}_r(n) \big) \hspace{0.5mm} \|\mathbf{x}\| \geq C^{1/2} \hspace{0.5mm} \sigma_1(R^*(n)) \hspace{0.5mm} \|\mathbf{x}\| \quad \textnormal{a.s.},
\end{equation*}

where the second inequality stems from (<ref>). Thus, since $P(n)$ and $R(n)$ in (<ref>) share the same singular values, condition (<ref>) implies that

\begin{equation}\label{e:||R*(n)Q*(n)Ur(n)x||}
\liminf_{n \rightarrow \infty}\|R^*(n)Q^*(n) \mathbf{U}_r(n) \mathbf{x}\| \geq C^{1/2}\lim_{n\rightarrow \infty}\sigma_1(R^*(n)) \hspace{0.5mm} \|\mathbf{x}\| \quad \textnormal{a.s.}
\end{equation}

In (<ref>), $\lim_{n\rightarrow \infty}\sigma_1(R^*(n)) > 0$ exists as a consequence of (<ref>). However, by condition (<ref>),

\begin{equation}\label{e:|Gamma*x|}
\|\Gamma\mathbf{x}\|= \plim_{n\to\infty} \|R^*(n)Q^*(n) \mathbf{U}_r(n) \mathbf{x}\|.
\end{equation}

Therefore, by expressions (<ref>), (<ref>) and (<ref>), as $n \rightarrow \infty$, for some $C' > 0$, $0 = \|\Gamma\mathbf{x}\| \geq C'\|\mathbf{x}\|$ with probability going to 1. Since $\Gamma$ is a constant matrix a.s., this implies that $\mathbf{x}=0$. Hence, (<ref>) is established.

To show (<ref>), once again let $\Gamma$ as in (<ref>) and fix $q \in \{1,\hdots,r\}$. If $\mathcal{I}_+=\{\ell:h_\ell>h_q\} \neq \emptyset$, Lemma <ref> implies that $\Gamma_{\ell,\ell'}=0$ whenever $\ell\in \mathcal{I}_-\cup\mathcal{I}_0$ and $\ell' \in \mathcal{I}_+$, since, in this case, $\ell<\ell'$ and $|\langle \p_{\ell'}(n), \mathbf u_{p-r+\ell}(n)\rangle a(n)^{h_{\ell'}-h_\ell}|$ is bounded, but $a(n)^{h_{\ell'}-h_\ell}\to\infty$. In other words, we can recast

\begin{equation}\label{e:Gamma_uppertri}
\Gamma = \begin{pmatrix}
(\Gamma_{\ell,\ell'})_{\ell,\ell'\in\mathcal{I}_-\cup\mathcal{I}_0} & (\Gamma_{\ell,\ell'})_{\ell\in\mathcal{I}_-\cup\mathcal{I}_0,\ell'\in\mathcal{I}_+} \\
\mathbf{0}& (\Gamma_{\ell,\ell'})_{\ell,\ell'\in\mathcal{I}_+}
\end{pmatrix}.
\end{equation}

For $\Gamma_+$ as in (<ref>), suppose $\Gamma^*_+\mathbf{y} = \mathbf{0}$ for some $\mathbf{y}\in\bbR^{r_3}$. Now let $\mathbf{x}=(\mathbf{0},\mathbf{y})^*\in \bbR^{r}$. Then, by (<ref>), $\Gamma^*\mathbf{x}=\mathbf{0}$. Thus, $\mathbf{x}=\mathbf{0}$. In other words, (<ref>) holds. $\Box$

The following lemma is used in the proofs of Theorem <ref> and Proposition <ref>.

Fix $q\in\{1,\ldots,r\}$, and for $\ell = 1,\hdots,p$, let ${\boldsymbol {\mathfrak u}}_{\ell}(n)$ be a unit eigenvector associated with the $\ell$–th eigenvalue of the random matrix $\widetilde {\mathbf W}(a(n)2^j)$ as defined in  Lemma <ref>. For each $n$, let

\begin{equation}\label{e:def_gamma_n}
\boldsymbol \gamma_\ell(n):=P^*(n){\boldsymbol {\mathfrak u}}_{p-r+\ell}(n),\quad \ell=1,\ldots,r,
\end{equation}

be the column vectors of the matrix $\Gamma_n$ as in (<ref>). Suppose $n_{\mathfrak a}$ is a subsequence along which the limits

\begin{equation}\label{e:lim_na_P*(na)u(na)}
\plim_{n_{\mathfrak a} \to\infty} {\boldsymbol \gamma}_\ell (n_{\mathfrak a}) = {\boldsymbol \gamma}^{\mathfrak a}_\ell, \quad \ell=1,\ldots,r,
\end{equation}

exist, where ${\boldsymbol \gamma}^{\mathfrak a}_\ell$ is a constant vector a.s.  For any fixed scalars

\begin{equation}\label{e:alpha_q,...,alpha_r-r3}
\alpha_q,\ldots,\alpha_{r-r_3},
\end{equation}

\begin{equation}\label{e:gamma-tilde}
\widetilde{\boldsymbol \gamma} = \sum^{r-r_3}_{\ell=q}\alpha_\ell {\boldsymbol \gamma}_\ell^{\mathfrak a}.
\end{equation}

Fix a vector

\begin{equation}\label{e:x*}
\mathbf{x}_{*} = (x_{*,{r-r_3+1}},\hdots,x_{*,r})^* \in \bbR^{r_3}.
\end{equation}

Then, there exists a sequence of unit vectors ${\mathbf v}_{\mathfrak a}(n_{\mathfrak a})\in \textnormal{span}\{{\mathbf u}_{p-r+q}(n_{\mathfrak a}),\ldots,{\mathbf u}_{p}(n_{\mathfrak a})\}$ satisfying

\begin{equation}\label{e:<p,w>=infinitesimal_2}