Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
42,458 characters · 9 sections · 41 citation commands
A Uniform Bound on the Operator Norm of Sub-Gaussian Random Matrices and Its Applications
Keywords Random Matrix Theory, Operator Norm, Uniform Bound, Operator Norm Minimizing Estimator, Functional Factor Models.
\onehalfspacing
Since its introduction in nuclear physics Wigner1955 and mathematical statistics Wishart1928, random matrix theory has been developed to understand the properties of the spectra of large dimensional random matrices generated by various distributions. These include the asymptotic theory of the empirical distribution of the eigenvalues of large dimensional random matrices and bounds on the extreme eigenvalues. For detailed results on these topics, readers can refer to recent surveys like Bai2008, EdelmanRao2005, BaiSilverstein2010, and Tao2012, among others.
In random matrix theory the study of the asymptotics of the largest eigenvalue of large dimensional random matrices goes back to Geman1980. Suppose that $X$ is an $N \times T$ matrix consisting of random variables $x_{it}$. Many researchers have derived the limit of the largest eigenvalue of the sample covariance matrix, $\lambda_1(X'X)$\footnote{$X'$ denotes the transpose of $X$.}, under various distributional assumptions on the random matrix $X$. For example, when $X_{it}$ are iid $N(0,1)$ and $\kappa := \lim \frac{N}{T}$, Geman1980 showed that $\frac{1}{N} \lambda_1(X'X) \rightarrow_{a.s.} (1 + \kappa^{1/2})^2$. Johnstone2001 obtained a stronger result that the properly normalized largest eigenvalue, $\frac{\lambda_1(X'X) - \mu_{NT}}{\sigma_{NT}}$ with $\mu_{NT} = (\sqrt{N-1} + \sqrt{T} )^2$ and $\sigma_{NT} = (\sqrt{N-1} + \sqrt{T} )(1/\sqrt{N-1} + 1/\sqrt{T})^{1/3}$, converges to the Tracy--Widom law; this has been later shown to hold under more general distributional assumptions by khorunzhiy2012high and tao2011random, among many others.
The aforementioned results imply that $\lambda_1(X'X)$ is stochastically bounded\footnote{A sequence of random variables $\xi_n$ is said to be stochastically bounded or order $a_n$, $\xi_n=O_p(a_n)$, if for any $\varepsilon>0$ there exists $M>0$ such that $\mathbb{P}(|\xi_n/a_n|\geq M) \leq \varepsilon$ for all large enough values of $n$.} of order $\max(N,T)$, or equivalently, the operator norm $\| X \| := \sqrt{\lambda_1(X'X)}$ is stochastically bounded of order $\sqrt{\max(N,T)}$. In fact, such bound does not require that the underlying distribution is Gaussian and can be derived under much weaker conditions. For example, Latala2005 showed that the bound holds if $x_{it}$ are independent across $(i,t)$ with mean zero and uniformly bounded fourth moments. MoonWeidner2017 extended this result for the cases where $x_{it}$ are weakly correlated across $i$ or $t$. Other papers that have established similar bounds on $E\| X \| $ include Bandeira2016, Guedonetal2017 and latala2018dimension.
In the case where $X$ consists of independent sub-Gaussian entries, the $\sqrt{\max(N,T)}$ order for the operator norm may be obtained using a powerful way of bounding sub-Gaussian stochastic processes called generic chaining, which was developed in fernique1975regularite and advanced later by M. Talagrand in a series of papers. Indeed, note that $||X||=\max_{u\in\mathbf{U}}\max_{v\in\mathbf{V}} u'Xv$, where maxima are taken over the unit spheres $\mathbf{U} \subset \mathbb{R}^N$ and $\mathbf{V}\subset \mathbb{R}^T$, respectively. The process $Z(u,v)=u'Xv$ defined on $\mathbf{U}\times \mathbf{V}$ can be shown to be sub-Gaussian and so we can invoke generic chaining to get the bound for its expected maximum in terms of a certain measure of metric complexity of $\mathbf{U}\times \mathbf{V}$ called Talagrand's functional $\gamma_2(\mathbf{U}\times\mathbf{V})$ (see definition in the next section). It turns out that $\gamma_2(\mathbf{U}\times\mathbf{V})$ has exact order $\sqrt{\max(N,T)}$.
In this paper we extend existing nonasymptotic bounds on the operator norm of a high-dimensional random matrix to the case of elements that are allowed to be weakly dependent and to be functions of a possibly infinite-dimensional parameter. Specifically, let $x_{it}(\beta)$ be weakly dependent over $t$, sub-Gaussian stochastic processes indexed by parameter $\beta$ belonging to a (pseudo-)metric space $(\mathbf{B},d_\mathbf{B})$. Let $X(\beta)$ be the $N \times T$ matrix consisting of $x_{it}(\beta)$ and let $\gamma_2(\mathbf{B},d_\mathbf{B})$ be Talagrand's functional of $\mathbf{B}$ w.r.t. $d_{\mathbf{B}}$ . Our main contribution is to show that $\mathbb{E} \sup_{\beta \in \mathbf{B}} \| X(\beta) \| $ is of order $\sqrt{\max(N,T)}+ \gamma_2(\mathbf{B},d_\mathbf{B})$.
We illustrate usefulness of this uniform bound with two examples. In one, we propose and show consistency of a new estimator that minimizes the operator norm of a matrix that consists of moment functions. In the other, we consider the generalization of the standard factor model to the case of functional data and suggest a new estimator of the maximal number of factors.
The paper is organized as follows. Section (ref) introduces our uniform bound along with the techniques necessary for its derivation. Section (ref) contains two applications of our theoretical result. Finally, Section (ref) concludes the paper. The appendix contains two technical proofs of the results in the main text.
Throughout the paper, $C$ will denote a universal positive constant that may not be the same at each occurrence, but may never depend on sample sizes, dimensions or any other features of the modeling framework.
Our main result is based on the general bound on suprema of sub-Gaussian processes called the generic chaining bound. We discuss this classic technique in this section and provide a proof in the appendix for completeness.
First, we need the following definitions. The $\psi$-Orlicz norm of a random variable $Y$ is defined as
where $\psi:\mathbb{R}_{+}\to \mathbb{R}_{+}$ is a convex function satisfying $\lim_{x\to \infty} \psi(x)/x = \infty$ and $\lim_{x\to 0} \psi(x)/x = 0$, and the convention that the infimum of an empty set is $+\infty$. In this paper, we let $\psi=\psi_2$, where $\psi_2(x)=\exp(x^2)-1$, and call $||\cdot||_{\psi_2}$ just “the Orlicz norm”. A random variable with finite ($\psi_2$-)Orlicz norm is called sub-Gaussian.
Intuitively, the Orlicz norm quantifies the decay speed for the tails of the distribution of $Y$. In fact, $||Y||_{\psi_2} \leq K$ is equivalent to\footnote{See e.g. vershynin2018high, Proposition 2.5.2.} \[ \mathbb{P}(|Y|\leq t) \geq 1-2e^{-t^2/K^2} \text{ for all } t\geq 0. \] Hence, for example, Gaussian distributions and distributions with bounded support are all sub-Gaussian.
Note also that the last inequality implies
Now let $T$ be a set and $d$ be a (pseudo-)metric on this set such that $(T,d)$ is a (pseudo-)metric space\footnote{Throughout the paper, “metric” can be replaced by a less restrictive notion of “pseudometric”, a distinction we omit from now on.}. Consider a zero mean stochastic process $(Z_t)$ indexed by the elements of $T$. The process $(Z_t)$ is said to have sub-Gaussian increments if there exists a constant $K>0$ such that
It has long been understood that behavior of sub-Gaussian processes is intimately connected to the metric complexity of its index set. In particular, the conventional bound on the expected supremum of $(Z_t)$ (see e.g. van1996weak Corollary 2.2.8.) is
where $N(T,d,\varepsilon)$ is the covering number of $(T,d)$ (i.e. the minimal number of $\varepsilon$-balls that is sufficient to cover $T$ in metric $d$) and $C$ is an absolute constant. The integral on the right hand side is sometimes called Dudley's entropy of $(T,d)$ and quantifies complexity of $(T,d)$ across multiple scales.
It turns out, however, that Dudley's entropy bound is not optimal, even for Gaussian processes. In fact, the entropy may be infinite when the expected supremum is not, rendering the bound uninformative\footnote{For an illustrative example, see Exercise 8.1.12 in vershynin2018high.}.
This led to the development of more precise ways to control suprema of sub-Gaussian processes in fernique1975regularite and talagrand2006generic. The generic chaining bound is stronger than (ref) and is sharp for Gaussian processes\footnote{See Section 8.6 in vershynin2018high.}. To introduce it, we need another definition.
For a metric space $(T,d)$, a sequence of finite subsets $T_0\subset T_1\subset \cdots \subset T$ is admissible if their cardinalities satisfy
Let the distance from the point $t\in T$ to the set $T_k \subset T$ be \[ d(t,T_k) = \inf_{t'\in T_k} d(t,t'). \]
Talagrand's functional $\gamma_2$ is then defined by the formula
where the infimum is taken over all admissible sequences $(T_k)$. Note that we can restrict our attention to only those admissible sequences that eventually come arbitrarily close to any point $t\in T$, which is possible provided $(T,d)$ is separable\footnote{A metric space $(T,d)$ is separable if it has a countable subset that is dense in $T$.}.
To understand the relation between Talagrand's functional and Dudley's entropy, let us provide the discussion from talagrand2006generic pp.12--13 here.
Denote $N_0=1$, $N_k=2^{2^k}$ for $k\geq 1$, and \[ e_k(T)=\inf_{S \subset T: \,\, |S|\leq N_k} \sup_{t\in T} d(t,S). \] Note that
where the second equality holds because minimizing the sum $ \sum_{k=0}^\infty 2^{k/2} \sup_{t\in T} d(t,T_k)$ w.r.t. all admissible sequences $(T_k)$ can be performed by separately minimizing each term $\sup_{t\in T} d(t,T_k)$ w.r.t. subsets $T_k\subset T$ satisfying $|T_k|\leq N_k$.
The definition of $e_k(T)$ involves choosing at most $N_k$ points $S$ in $T$ such that the balls with radius $e_k(T)$ and centers in $S$ cover $T$; moreover, $e_k(T)$ is the minimal such radius, i.e. \[ e_k(T) = \inf\left\{\varepsilon>0:\,\, N(T,d,\varepsilon)\leq N_k\right\}. \] It follows that if $e_k(T)<\varepsilon$, then $N(T,d,\varepsilon)>N_k$ or $N(T,d,\varepsilon)\geq N_k+1$. Hence we can write
Since $\log(N_k+1)\geq 2^k \log 2$ for $k\geq 0$, summation over $k\geq 0$ yields
where, of course, $e_0(T)=\text{diam}(T)=\sup_{t,s\in T} d(t,s)$.
The term on the left hand side of this inequality satisfies
Combining this with (ref) and (ref) yields the key relation
Hence, when used as an upper bound, Talagrand's functional is sharper than Dudley's entropy.
We are now ready to state the generic chaining bound for sub-Gaussian processes, see e.g. Theorem 8.5.3 in vershynin2018high.
We impose the following assumptions.
where $\psi_{i\tau}(\beta)$ are nonrandom coefficients such that, for all $i=1,\dots,N$ and $\beta\in\mathbf{B}$,
(ref) is very weak and only imposes separability of the metric space $\mathbf{B}$ which holds for most parameter spaces encountered in practice such as Euclidean spaces and spaces of integrable functions. (ref) is similar to case (ii) in Lemma S.2.1 of MoonWeidner2017 and allows $x_{it}(\beta)$ to be weakly dependent over time. (ref) and (ref) impose uniform sub-Gaussianity on the innovations $\varepsilon_{it}(\beta)$ and their increments $\varepsilon_{it}(\beta_1)-\varepsilon_{it}(\beta_2)$, respectively. Note that (ref) is equivalent to the tail bound \[ \mathbb{P}\left( |\varepsilon_{it}(\beta_1)-\varepsilon_{it}(\beta_2)| \leq t \cdot d_\mathbf{B}(\beta_1,\beta_2) \right) \geq 1-2e^{-\frac{t^2}{K_2^2}} \text{ for all } t\geq 0. \]
Denote $\psi_\tau(\beta) = (\psi_{1\tau}(\beta),\dots,\psi_{N\tau}(\beta))'$ and let $\Xi_{-\tau}(\beta)$ the $N\times T$ matrix consisting of $\varepsilon_{it}(\beta)$, $i=1,\dots,N$, $t=1-\tau,\dots,T-\tau.$ Equation (ref) can be rewritten in the matrix form as
Suppose for a moment that we have a bound on $\Xi_{-\tau}(\beta)$ of the form \[ \mathbb{E} \sup_\beta ||\Xi_{-\tau}(\beta)|| \leq \varphi(N,T,\mathbf{B}), \] where $\varphi$ does not depend on $\tau$. Then
This shows that the bound on $\mathbb{E} \sup_\beta ||X(\beta)||$ is, up to the absolute constant $D$, the same as the bound on $\mathbb{E} \sup_\beta ||\Xi_{-\tau}(\beta)||$. Hence we can focus on obtaining the latter bound from now on. It will be clear from the proof that the bound will not depend on $\tau$, so we consider the case $\tau=0$ and denote $\Xi=\Xi_0$ for brevity.
The operator norm of $\Xi(\beta)$ can be expressed as \[ ||\Xi(\beta)||=\sup_{u\in \mathbf{U},v \in \mathbf{V}} Z(u,v,\beta), \] where $\mathbf{U}$ and $\mathbf{V}$ are unit spheres in $\mathbb{R}^N$ and $\mathbb{R}^T$, respectively, and the process \[ Z(u,v,\beta):=u'\Xi(\beta)v = \sum_{i=1}^N \sum_{t=1}^T u_i v_t \varepsilon_{it}(\beta). \] Define the $L_1$ product metric on $\mathbf{U}\times\mathbf{V}\times\mathbf{B}$ by \[ d((\tilde{u},\tilde{v}, \tilde{\beta}), (u,v,\beta)) = d_{\mathbb{R}^N}(\tilde{u},u) + d_{\mathbb{R}^T}(\tilde{v},v) + d_\mathbf{B}(\tilde{\beta},\beta). \] where $d_{\mathbb{R}^d}$ denotes the standard Euclidean metric on $\mathbb{R}^d$.
To obtain a uniform bound on $||\Xi(\beta)||$, we would like to apply (ref) to the process $Z(\cdot)$ defined on the metric space $(\mathbf{U}\times \mathbf{V}\times \mathbf{B}, d)$. Our first lemma asserts that $Z$ has sub-Gaussian increments.
Our second lemma establishes the bound on Talagrand's functional of a product space in terms of Talagrand's functionals of component spaces.
Finally, by (ref), we can apply the generic chaining bound of (ref) to $Z(u,v,\beta)$ defined on the separable metric space $T=\mathbf{U}\times \mathbf{V} \times \mathbf{B}$ with the $L_1$ metric $d$. (ref) then yields
For the unit sphere $S^{d-1}$ in $\mathbb{R}^d$, its Dudley's entropy satisfies \[ \int_0^{\text{diam}(S^{d-1})} \sqrt{\log N(S^{d-1},||\cdot||,\varepsilon)} \,d\varepsilon \leq C \sqrt{d}. \] Besides, Talagrand's functional is bounded from above by Dudley's entropy (e.g. Exercise 8.5.7 in vershynin2018high), up to absolute constant factors.
Applying these bounds to unit spheres $\mathbf{U}\subset \mathbb{R}^N$ and $\mathbf{V}\subset \mathbb{R}^T$ gives \[ \mathbb{E} \sup_{\beta \in \mathbf{B}} ||\Xi(\beta)|| \leq CK\left( \sqrt{\max(N,T)} + \gamma_2(\mathbf{B},d_\mathbf{B}) \right). \]
Finally, taking into account the inequality (ref), we obtain the main theoretical result of this paper.
{\bf Remarks}
In this section, we investigate a new estimator that minimizes the operator norm of the moment function matrix. Suppose that $\varepsilon_{it}(\beta) \in \mathbb{R}^L $ are $L$ moment functions of $\beta \in \mathbf{B} \subset \mathbb{R}^K$ such that $\mathbb{E}(\varepsilon_{it}(\beta_0)) = 0$. For simplicity, assume that $L = K = 1$. Let $\varepsilon(\beta) = [\varepsilon_{it}(\beta)]$, the $N \times T$ matrix of moment functions.
The conventional method of moment estimator solves \[ \tilde{\beta} = \arg\min_{\beta \in \mathbf{B}} \left| \frac{1}{NT} \sum_{i,t} \varepsilon_{it}(\beta) \right| = \arg\min_{\beta \in \mathbf{B}} \left| \frac{\mathbf{1}_N^{\prime}}{\sqrt{N}} \left(\frac{ \varepsilon(\beta)}{\sqrt{NT}} \right) \frac{\mathbf{1}_T}{\sqrt{T}} \right|, \] where $\mathbf{1}_N$ is the $N$-vector of ones.
The new estimator we propose minimizes the operator norm of the moment function matrix $\varepsilon(\beta)$,
In this section we establish consistency of $\widehat{\beta}$ using our main result of the previous section.
Conditions (i)-(ii) of Assumption (ref) ensure that $\varepsilon_{it}(\beta) - \mathbb{E}( \varepsilon_{it}(\beta))$ satisfies Assumptions (ref)-(ref). The last condition (iii) corresponds to the identification condition of the extremum estimator.
For consistency of $\widehat{\beta}$, it suffices to show that for any $\epsilon>0$, there exists $\delta > 0$ such that
with probability approaching one.
First, note that, since $\mathbb{E}( \varepsilon(\beta_0))=0$, the triangle inequality yields
On the other hand,
Combine (ref) and (ref) to obtain
Finally, choose $\delta$ as in Assumption (ref)(iii) to guarantee \[ \inf_{ | \beta - \beta_0 | \geq \epsilon} \frac{\| \mathbb{E}(\varepsilon(\beta)) \|}{\sqrt{NT}} \geq 2\delta \] and note that (ref) gives \[ \sup_{\beta \in \mathbf{B}} \frac{\| \varepsilon(\beta) - \mathbb{E}(\varepsilon(\beta) \|}{\sqrt{NT}} = o_p(1). \] Then (ref) implies
which finishes the proof of consistency of $\hat{\beta}$.
{\bf Remarks}
Consider a generic factor model for functional data
where $\beta$ belongs to a separable metric space $(\mathbf{B},d_\mathbf{B})$, $Y(\beta)\in \mathbb{R}^{N\times T}$ is the observation matrix of functional outcomes $\beta \mapsto y_{it}(\beta)$, and $\lambda(\beta)\in \mathbb{R}^{N\times R(\beta)}$, $f(\beta)\in \mathbb{R}^{ T \times R(\beta)}$ such that for all $\beta\in \mathbf{B}$ the probability limits of $\lambda(\beta)'\lambda(\beta)/N$ and $f(\beta)'f(\beta)/T$ exist and are positive definite deterministic matrices such that
The object of interest is the maximal rank $R=\max_{\beta\in \mathbf{B}}R(\beta)$.
To illustrate applicability of this model, suppose that the outcome variable is intraday pollution levels $y_{it}(\beta)$, where $\beta$ is the time within a day, across counties $i$ and time $t$, as in aue2015prediction. It is plausible to assume that counties with higher population density and dependence on automobiles will have higher average levels of pollution. At the same time, pollution patterns on weekdays and on weekends may differ in a systematic way. Hence it is reasonable to model the intraday pollution curve $y_{it}(\cdot)$ as the interaction of the county fixed effect $\lambda_i(\cdot)$ and the time effect $f_t(\cdot)$, plus independent noise, arriving at model (ref). A related approach to modeling functional time series can be found in kargin2008curve, whose empirical objective is to predict the contract rate curves of daily Eurodollar futures.
Of course, arguments similar to those outlined above may be applied to modeling of numerous other functional quantities, from mortality as a function of age to crop yields as a function of spatial location. For more examples and an overview of functional data analysis, see e.g. wang2016functional and kowal2019functional.
Let us now show heuristically how to derive a consistent estimator of the maximal rank $R$.
Note that the model assumptions imply
If $U(\beta)$ satisfies the conditions of (ref), we have $\sup_\beta ||U(\beta)|| =O_p\left( \sqrt{\max(N,T)} + \gamma_2(\mathbf{B},d_\mathbf{B})\right)$ and so \[ \sup_\beta ||U(\beta)|| = O_p\left(\sqrt{\max(N,T)}\right). \] Denote $s_i(A)$ the $i$-th largest singular value of matrix $A$. The Ky Fan inequality for singular values asserts that for $A,B \in \mathbb{R}^{N \times T}$ \[ |s_i(A+B)-s_i(A)| \leq s_1(B)=||B|| \text{ for all } i=1,\dots,\min(N,T). \]
Using this inequality, for a fixed $\beta$ we obtain
Therefore, there exists a positive constant $C>0$ such that
where the last inequality holds by ((ref)).
On the other hand,
This establishes consistency of the following natural estimator of $R$,
where $\psi_{NT}$ is a sequence of real numbers satisfying $\psi_{NT}\to 0$ and $\psi_{NT} \sqrt{\min(N,T)} \to \infty$.
Empirical practice calls for an automatic procedure for choosing the tuning parameter $\psi_{NT}$. One may consider one of the following three options, using the penalty term from bai2002determining:
where $\hat{\sigma}^2=\sup_\beta \hat{\sigma}^2(\beta)$ is a consistent estimator of \[ \sigma^2=\sup_\beta \sigma^2(\beta) = \sup_\beta \frac{1}{NT}\sum_{i=1}^N \sum_{t=1}^T \mathbb{E}\left[ u_{it}(\beta)^2 \right]. \] In applications, $\hat{\sigma}^2(\beta)$ can be replaced by the residual variance of $Y(\beta)$ after partialling out $k_{\max}$ factors using principle component analysis, where $k_{\max}$ is a pre-specified upper bound on the true maximal number of factors $R$.
Here we illustrate the performance of the maximal rank estimator in the functional factor model described in the previous section with a simple simulation design.
The data generating process is the functional factor model (ref), where, for simplicity, we let the loadings $\lambda(\beta)$ and the factors $f(\beta)$ to be independent of $\beta$. In scalar form, the model is
where $\lambda_{ir}, f_{tr} \sim \text{iid} \, N(0,1)$ and
The chosen specification for $u_{it}(\cdot)$ comes from a generic representation of any Gaussian stochastic process as an infinite trigonometric series, in which we only retain one term. Clearly, the error variance $\mathbf{V}(u_{it}(\beta))=\sigma$ for all $\beta$ and there is nontrivial dependence of $u_{it}(\beta)$ across values of $\beta$. We set $\sigma=1$. The results do not change substantially when larger values of $\sigma$ are used.
We choose the range of parameter $\beta$ to be $\mathbf{B} = \{0,0.1,\dots,0.9,1\}$ and the corresponding ranks \[ (R(0),R(0.1),\dots,R(0.9), R(1)) = (4,4,1,4,3,1,2,3,4,4,1), \] so that the true value of interest is $R=\max_\beta R(\beta)=4$.
The simulated bias and root MSE for the maximal rank estimator (ref) are shown in (ref). Clearly, the choice $\psi_{NT}=\psi_{NT,3}$ for the tuning parameter leads to poor small sample performance, which is similar to the results of bai2002determining. However, under the other two choices $\psi_{NT,1}, \psi_{NT,2}$, bias and RMSE are modest even in small samples and become essentially zero when $N,T \ge 50$.
Given these simulation results, we are convinced that our generalization (ref) of the estimator of bai2002determining will be useful for practitioners who are interested in estimating factor models with functional data.
In this paper, we derive a novel uniform stochastic bound on the operator norm of sub-Gaussian random matrices. We use it to establish consistency of a new estimator that minimizes the operator norm of the matrix of moment functions as well as to introduce an estimator of the maximal number of factors in a functional interactive fixed effects model.