High dimensional asymptotics of likelihood ratio tests in the Gaussian sequence model under convex constraints

Consider the Gaussian sequence model

where $\mu\in\mathbb{R}^{n}$ is unknown and $\xi=(\xi_{1},\ldots,\xi_{n})$ is an $n$-dimensional standard Gaussian vector. In a variety of applications, prior knowledge on the mean vector $\mu$ can be naturally translated into the constraint $\mu\in K$, where $K$ is a closed convex set in $\mathbb{R}^{n}$. Two such important examples that will be considered in this paper are: (i) Lasso in the constrained form [Tib96], where $K$ is an $\ell_{1}$-norm ball, and (ii) isotonic regression [CGS15], where $K$ is the cone consisting of monotone sequences. We also refer the readers to [Bar02, JN02, BHL05, Cha14, GS18] and many references therein for a diverse list of further concrete examples of $K$. In this paper, we will be interested in the following ‘goodness-of-fit’ testing problem:

where $\mu_{0}\in K\subset\mathbb{R}^{n}$, and $K$ is an arbitrary closed and convex subset of $\mathbb{R}^{n}$. Throughout the manuscript, the asymptotics will take place as $n\to\infty$, and the explicit dependence of $\mu,\mu_{0},K$ and related quantities on the dimension $n$ will be suppressed for ease of notation.

Given observation $Y$ generated from model (1.1), arguably the most natural and generic test for (1.2) is the likelihood ratio test (LRT). Under the Gaussian model (1.1), the log-likelihood ratio statistic (LRS) for (1.2) takes the form

Here $\widehat{\mu}_{K}\equiv\Pi_{K}(Y)\equiv\arg\min_{\nu\in K}\lVert Y-\nu\rVert^{2}$ is the metric projection of $Y$ onto the constraint set $K$ with respect to the canonical Euclidean $\ell_{2}$ norm $\|\cdot\|$ on $\mathbb{R}^{n}$. As $K$ is both closed and convex, $\Pi_{K}$ is well-defined, and the resulting $\widehat{\mu}_{K}$ is both the least squares estimator (LSE) and the maximum likelihood estimator for the mean vector $\mu$ under the Gaussian model (1.1). The risk behavior of $\widehat{\mu}_{K}$ is completely characterized in the recent work [Cha14].

The LRT for (1.2) and its generalizations thereof have gained extensive attention in the literature, see e.g. [Che54, Bar59a, Bar59b, Bar61a, Bar61b, Kud63, BBBB72, KC75, RW78, WR84, Sha85, RLN86, RWD88, Sha88, MS91, MRS92a, MRS92b, DT01, Mey03, SM17, WWG19] for an incomplete list. In our setting, an immediate way to use the LRS $T(Y)$ in (1.1) to form a test is to simulate the critical values of $T(Y)$ under $H_{0}$. More precisely, for any confidence level $\alpha\in(0,1)$, we may determine through simulations an acceptance region $\mathcal{I}_{\alpha}\subset\mathbb{R}$ such that the LRS satisfies $\mathbb{P}\big{(}T(Y)\in\mathcal{I}_{\alpha}\big{)}=1-\alpha$ under $H_{0}$, and then formulate the LRT accordingly. In some special cases, including the classical setting where $K$ is a subspace, the distribution of $T(Y)$ under the null is even known in closed form, so the simulation step can be skipped.

Clearly, the almost effortless LRT as described above already gives an exact type I error control at the prescribed level for a generic pair $(\mu_{0},K)$. The equally important question of its power behavior, however, is more complicated and requires a much deeper level of investigation. In the classical setting of parametric models and certain semiparametric models, the power behavior of the LRT can be precisely determined, at least asymptotically, for contiguous alternatives in the corresponding parameter spaces, cf. [vdV98, vdV02]. An important and basic ingredient for the success of power analysis in these settings is the existence of a limiting distribution of the LRT under $H_{0}$ that can be ‘perturbed’ in a large number of directions of alternatives.

Unfortunately, the distribution of the LRS $T(Y)$ in (1.1) under the null, for both finite-sample and asymptotic regimes, is only understood in very few cases. One such case is, as mentioned above, the classical setting where $K$ is a subspace of dimension $\dim(K)$. Then the null distribution of $T(Y)$ is a chi-squared distribution with $\dim(K)$ degrees of freedom. Another case is $\mu_{0}=0$ and $K$ is a closed convex cone. In this case, the null distribution of $T(Y)$ is the chi-bar squared distribution, cf. [Bar61a, Kud63, BBBB72, KC75, Sha85, RWD88], which can be expressed as a finite mixture of chi-squared distributions. Apart from these special cases, little next to nothing is known about the distribution of the LRS $T(Y)$ for a general pair of $(\mu_{0},K)$ under the null $H_{0}$, owing in large part to the fact that the null distribution of $T(Y)$ highly depends on the exact location of $\mu_{0}$ with respect to $K$ and is thus intractable in general. Consequently, the lack of such a general description of the limiting distribution of $T(Y)$ causes a fundamental difficulty in obtaining precise characterizations of the power behavior of the LRT for a general pair $(\mu_{0},K)$. On the other hand, such generality is of great interest as it allows us to consider several significant examples, for instance testing general signals in isotonic regression and with constrained Lasso. See Section 4 for more details.

The unifying theme of this paper starts from the simple observation that in the classical setting where $K$ is a subspace, as long as $\dim(K)$ diverges as $n\rightarrow\infty$, the distribution of $T(Y)$ has a progressively Gaussian shape, under proper normalization. Such a normal approximation in ‘high dimensions’ also holds for the more complicated chi-bar squared distribution; see [Dyk91, GNP17] for a different set of conditions. One may therefore hope that normal approximation of $T(Y)$ under the null would hold in a far more general context than just these cases. More importantly, such a distributional approximation would ideally form the basis for power analysis of the LRT.

A particularly important special setting for (1.2) is the case of testing $H_{0}:\mu=0$ versus $H_{1}:\mu\in K$, where $K$ is assumed to be a closed convex cone in $\mathbb{R}^{n}$. We perform a detailed case study of the following slightly more general testing problem:

where $K_{0}\subset K\subset\mathbb{R}^{n}$ is a subspace, and $K$ is a closed convex cone. The primary motivation to study (1.10) arises from the problem of testing a global polynomial structure versus its shape-constrained generalization; concrete examples include constancy versus monotonicity, linearity versus convexity, etc.; see Section 4.5 for details. The LRS for (1.10) takes the slightly modified form:

The dependence in the notation of the LRS $T(Y)$ on $K_{0}$ will usually be suppressed when no confusion could arise.

Specializing our first main result to this testing problem, we show in Theorem 3.8 that normal approximation of $T(Y)$ under $H_{0}$ holds essentially under the minimal growth condition that $\delta_{K}-\dim(K_{0})\to\infty$, where $\delta_{K}$ is the statistical dimension of $K$ (formally defined in Definition 2.2). Similar to (1.8), the normal approximation makes possible the following precise characterization of the power behavior of the LRT under the prescribed growth condition (see Theorem 3.9):

As $\sigma_{0}^{2}=\operatorname{Var}(T(\xi))\asymp\delta_{K}-\dim(K_{0})$ (cf. Lemma 2.4) for the modified LRS $T(Y)$ in (1.3), the LRT is power consistent under $\mu$ if and only if

Formula (1.13) shows that power consistency of the LRT is determined completely by (a complicated expression involving) the ‘distance’ of the alternative $\mu\in K$ to its projection onto $K_{0}$ in the problem (1.10). Compared to the uniform $\|\cdot\|$-separation rate derived in the recent work [WWG19] (cf. (3.19) below), (1.3)-(1.13) provide asymptotically precise power characterizations of the LRT for a sequence of point alternatives. This difference is indeed crucial as (1.13), similar to (1.9), cannot be equivalently inverted into a lower bound on $\lVert\mu-\Pi_{K_{0}}(\mu)\rVert$ alone. This illustrates that the non-uniform power behavior of the LRT is not an aberration in certain artificial testing problems, but is rather a fundamental property of the LRT in the high dimensional regime that already appears in the special yet important setting of testing subspace versus a cone.

As an illustration of the scope of our theoretical results, we validate the normal approximation of the LRT and exemplify its power behavior in two classes of problems:

1. (1)

   Testing in orthant/circular cone, isotonic regression and Lasso;

2. (2)

   Testing parametric assumptions versus shape-constrained alternatives, e.g., constancy versus monotonicity, linearity versus convexity, and generalizations thereof.

The results in this paper are related to the vast literature on nonparametric testing in the Gaussian sequence model, or more general Gaussian models, under a minimax framework. We refer the readers to the monographs [IS03, GN16] for a comprehensive treatment on this topic, and [Bar02, JN02, DJ04, BHL05, ITV10, ACCP11, Ver12, CD13, Car15, CCT17, CCC⁺19, CV19, CCC⁺20, MS20] and references therein for some recent papers on a variety of testing problems. Many results in these references establish minimax separation rates under a pre-specified metric, with the Euclidean metric $\lVert\cdot\rVert$ being a popular choice in the Gaussian sequence model. In particular, for the testing problem (1.10), this minimax approach with $\lVert\cdot\rVert$ metric is adopted in the recent work [WWG19], which derived minimax lower bounds for the separation rates and the uniform separation rate of the LRT. These results show that the LRT is minimax rate-optimal in a number of examples, while being sub-optimal in some other examples.

Our results are of a rather different flavor and give a precise distributional description of the LRT. Such a description is made possible by the central limit theorems of the LRS under the null proved in Theorems 3.1 and 3.8. It also allows us to make two significant further steps beyond the work [WWG19]:

1. (1)

   For the testing problem (1.10), we provide asymptotically exact power formula of the LRT in (1.3) for each and every alternative, as opposed to lower bounds for the uniform separation rates of the LRT in $\lVert\cdot\rVert$ as in [WWG19]. As a result, the main results for the separation rates of the LRT in [WWG19] follow as a corollary of our main results (see Corollary 3.11 for a formal statement).

2. (2)

   Our theory applies to the general testing problem (1.2) which allows for a general pair $(\mu_{0},K)$ without a cone structure. This level of generality goes much beyond the scope of [WWG19] and covers several significant examples, including testing in isotonic regression and Lasso problems.

The precise power characterization we derive in (1.3) has interesting implications when compared to the minimax results derived in [WWG19]. In particular, as discussed in Section 1.4.1, (i) it is possible that the LRT beats substantially the minimax separation rates in $\lVert\cdot\rVert$ metric for individual alternatives, and (ii) the sub-optimality of the LRT in [WWG19] is actually only witnessed at alternatives along some ‘bad’ directions. In this sense, our results not only give precise understanding of the power behavior of the canonical LRT in this testing problem, but also highlights some intrinsic limitations of the popular minimax framework under the Euclidean metric in the Gaussian sequence model, both in terms of its optimality and sub-optimality criteria.

From a technical perspective, our proof technique differs significantly from the one adopted in [WWG19]. Indeed, the proofs of the central limit theorems and the precise power formulae in this paper are inspired by the second-order Poincaré inequality due to [Cha09] and related normal approximation results in [GNP17]. These technical developments are of independent interest and have broader applicability; see for instance [HJS21] for further developments in the context of testing high dimensional covariance matrices.

The rest of the paper is organized as follows. Section 2 reviews some basic facts on metric projection and conic geometry. Section 3 studies normal approximation for the LRS $T(Y)$ and the power characterizations of the LRT both in the general setting (1.2) and the more structured setting (1.10). Applications of the abstract theory to the examples mentioned above are detailed in Section 4. Proofs are collected in Sections 5 and 6 and the appendix.

For any positive integer $n$, let $[1:n]$ denote the set $\{1,\ldots,n\}$. For $a,b\in\mathbb{R}$, $a\vee b\equiv\max\{a,b\}$ and $a\wedge b\equiv\min\{a,b\}$. For $a\in\mathbb{R}$, let $a_{\pm}\equiv(\pm a)\vee 0$. For $x\in\mathbb{R}^{n}$, let $\lVert x\rVert_{p}$ denote its $p$-norm $(0\leq p\leq\infty)$, and $B_{p}(r;x)\equiv\{z\in\mathbb{R}^{n}:\lVert z-x\rVert_{p}\leq r\}$. We simply write $\lVert x\rVert\equiv\lVert x\rVert_{2}$, $B(r;x)\equiv B_{2}(r;x)$, and $B(r)\equiv B(r;0)$ for notational convenience. By $\bm{1}_{n}$ we denote the vector of all ones in $\mathbb{R}^{n}$. For a matrix $M\in\mathbb{R}^{n\times n}$, let $\lVert M\rVert$ and $\lVert M\rVert_{F}$ denote the spectral and Frobenius norms of $M$ respectively.

For a multi-index $\bm{k}=(k_{1},\ldots,k_{n})\in\mathbb{Z}_{\geq 0}^{n}$, let $\lvert\bm{k}\rvert\equiv\sum_{i=1}^{n}k_{i}$. For $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$, and $\bm{k}=(k_{1},\ldots,k_{n})\in\mathbb{Z}_{\geq 0}^{n}$, let $\partial_{\bm{k}}f(z)\equiv\frac{\partial^{\lvert\bm{k}\rvert}f(z)}{\partial_{k_{1}}z_{1}\cdots\partial_{k_{n}}z_{n}}$ for $z\in\mathbb{R}^{n}$ whenever definable. A vector-valued map $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ is said to have sub-exponential growth at $\infty$ if $\lim_{\lVert x\rVert\to\infty}\lVert f(x)e^{-\lVert x\rVert}\rVert=0$. For $f=(f_{1},\ldots,f_{n}):\mathbb{R}^{n}\to\mathbb{R}^{n}$, let $J_{f}(z)\equiv(\partial f_{i}(z)/\partial z_{j})_{i,j=1}^{n}$ denote the Jacobian of $f$ and

for $z\in\mathbb{R}^{n}$ whenever definable.

We use $C_{x}$ to denote a generic constant that depends only on $x$, whose numeric value may change from line to line unless otherwise specified. $a\lesssim_{x}b$ and $a\gtrsim_{x}b$ mean $a\leq C_{x}b$ and $a\geq C_{x}b$ respectively, and $a\asymp_{x}b$ means $a\lesssim_{x}b$ and $a\gtrsim_{x}b$ ($a\lesssim b$ means $a\leq Cb$ for some absolute constant $C$). For two nonnegative sequences $\{a_{n}\}$ and $\{b_{n}\}$, we write $a_{n}\ll b_{n}$ (respectively $a_{n}\gg b_{n}$) if $\lim_{n\rightarrow\infty}(a_{n}/b_{n})=0$ (respectively $\lim_{n\rightarrow\infty}(a_{n}/b_{n})=\infty$). We follow the convention that $0/0=0$. $\mathcal{O}_{\mathbf{P}}$ and $\mathfrak{o}_{\mathbf{P}}$ denote the usual big and small O notation in probability.

We reserve the notation $\xi=(\xi_{1},\ldots,\xi_{n})$ for an $n$-dimensional standard normal random vector, and $\varphi,\Phi$ for the density and the cumulative distribution function of a standard normal random variable. For any $\alpha\in(0,1)$, let $z_{\alpha}$ be the normal quantile defined by $\mathbb{P}(\mathcal{N}(0,1)\geq z_{\alpha})=\alpha$. For two random variables $X,Y$ on $\mathbb{R}$, we use $d_{\mathrm{TV}}(X,Y)$ and $d_{\mathrm{Kol}}(X,Y)$ to denote their total variation distance and Kolmogorov distance defined respectively as

Here $\mathcal{B}(\mathbb{R})$ denotes all Borel measurable sets in $\mathbb{R}$.

We start by presenting the normal approximation result for $T(Y)$ in (1.1) under the null hypothesis (1.2); see Section 5.1 for a proof. This will serve as the basis for the size and power analysis of the LRT (1.6) in the testing problem (1.2).

The bound (3.1) is obtained by a generalization of [GNP17, Theorem 2.1] using the second-order Poincaré inequality [Cha09], together with a lower bound for $\sigma_{\mu_{0}}^{2}$ using Fourier analysis in the Gaussian space [NP12, Section 1.5]. The Fourier expansion can be performed up to the second order thanks to the absolute continuity of the first-order partial derivatives of $T(Y)$ (cf. Lemma 2.1).

We now comment on the structure of (3.1). The first term $\lVert\mathbb{E}_{\mu_{0}}\widehat{\mu}_{K}-\mu_{0}\rVert^{2}$ in the denominator is the squared bias of the projection estimator $\widehat{\mu}_{K}$, while the second term $\lVert\mathbb{E}_{\mu_{0}}J_{\widehat{\mu}_{K}}\rVert_{F}^{2}$, which depends on the magnitudes of the first-order partial derivatives of $\widehat{\mu}_{K}$, can be roughly understood as the ‘variance’ of $\widehat{\mu}_{K}$. Consequently, one may expect that the denominator is of the order $\mathbb{E}_{\mu_{0}}\lVert\widehat{\mu}_{K}-\mu_{0}\rVert^{2}$, so the overall bound scales as $\mathcal{O}\big{(}1/\sqrt{\mathbb{E}_{\mu_{0}}\lVert\widehat{\mu}_{K}-\mu_{0}\rVert^{2}}\big{)}$. As will be clear in Section 4, this is indeed the case in all the examples worked out, and the major step in applying (3.1) to concrete problems typically depends on obtaining sharp lower bounds for the ‘variance’ term $\lVert\mathbb{E}_{\mu_{0}}J_{\widehat{\mu}_{K}}\rVert_{F}^{2}$, which may require non-trivial problem-specific techniques.

Using Theorem 3.1, we may characterize the size and power behavior of the LRT. For a possibly unbounded interval $I\subset\mathbb{R}$, let $\Delta_{I}:\overline{\mathbb{R}}\to[0,1]$ be defined as follows: For $w\in\mathbb{R}$,

and $\Delta_{I}(\pm\infty)\equiv\lim_{w\to\pm\infty}\Delta_{I}(w)$, which is clearly well-defined. $\Delta_{I}$ is either monotonic or unimodal, so $\Delta_{I}^{-1}(\beta)$ contains at most two elements for any $\beta\in[0,1]$. Two primary examples of $I$ are $\mathcal{A}^{\textrm{os}}_{\alpha}\equiv(-\infty,z_{\alpha}]$ and $\mathcal{A}^{\textrm{ts}}_{\alpha}\equiv[-z_{\alpha/2},z_{\alpha/2}]$ — the acceptance regions for the one- and two-sided LRTs respectively, where we have

It is clear that $\Delta_{\mathcal{A}^{\textrm{os}}_{\alpha}}(0)=\Delta_{\mathcal{A}^{\textrm{ts}}_{\alpha}}(0)=\alpha$, $\Delta^{-1}_{\mathcal{A}^{\textrm{os}}_{\alpha}}(1)=\{+\infty\}$, and $\Delta^{-1}_{\mathcal{A}^{\textrm{ts}}_{\alpha}}(1)=\{\pm\infty\}$. Recall the definitions of $m_{\mu}$ and $\sigma^{2}_{\mu}$ for general $\mu\in K$ in (1.5). The following result (see Section 5.2 for a proof) characterizes the power behavior of the LRT.

As a simple toy example of Theorem 3.2, we consider the testing problem (1.2) in the linear regression case, where $K\equiv K_{X}\equiv\{X\theta:\theta\in\mathbb{R}^{p}\}$ for some fixed design matrix $X\in\mathbb{R}^{n\times p}$, with $p\leq n$. We will be interested in the high dimensional regime $\operatorname{rank}(X)\to\infty$ where the normal approximation for the LRT holds under the null.