The Likelihood Ratio Test in High-Dimensional Logistic Regression Is Asymptotically a Rescaled Chi-Square

TL;DR

This paper demonstrates that in high-dimensional logistic regression, the likelihood ratio test statistic converges to a rescaled chi-square distribution, challenging classical assumptions and providing new correction factors for p-value calculations.

Contribution

It proves that Wilks' theorem fails in high-dimensional settings and derives a new asymptotic distribution for the likelihood ratio test, including a method to compute the scaling factor.

Findings

Likelihood ratio statistic is asymptotically a rescaled chi-square distribution.

The scaling factor exceeds one when the ratio of variables to samples is positive.

The new distribution accurately predicts finite-sample behavior.

Abstract

Logistic regression is used thousands of times a day to fit data, predict future outcomes, and assess the statistical significance of explanatory variables. When used for the purpose of statistical inference, logistic models produce p-values for the regression coefficients by using an approximation to the distribution of the likelihood-ratio test. Indeed, Wilks' theorem asserts that whenever we have a fixed number $p$ of variables, twice the log-likelihood ratio (LLR) $2\Lambda$ is distributed as a $\chi^2_k$ variable in the limit of large sample sizes $n$; here, $k$ is the number of variables being tested. In this paper, we prove that when $p$ is not negligible compared to $n$, Wilks' theorem does not hold and that the chi-square approximation is grossly incorrect; in fact, this approximation produces p-values that are far too small (under the null hypothesis). Assume that $n$ and $p$ grow large in such a way that $p/n\rightarrow\kappa$ for some constant $\kappa < 1/2$. We prove that for a class of logistic models, the LLR converges to a rescaled chi-square, namely, $2\Lambda~\stackrel{\mathrm{d}}{\rightarrow}~\alpha(\kappa)\chi_k^2$, where the scaling factor $\alpha(\kappa)$ is greater than one as soon as the dimensionality ratio $\kappa$ is positive. Hence, the LLR is larger than classically assumed. For instance, when $\kappa=0.3$, $\alpha(\kappa)\approx1.5$. In general, we show how to compute the scaling factor by solving a nonlinear system of two equations with two unknowns. Our mathematical arguments are involved and use techniques from approximate message passing theory, non-asymptotic random matrix theory and convex geometry. We also complement our mathematical study by showing that the new limiting distribution is accurate for finite sample sizes. Finally, all the results from this paper extend to some other regression models such as the probit regression model.

Tables

Table 1. Table 1 : Estimates of p-value probabilities with estimated Monte Carlo standard errors in parentheses under i.i.d. Gaussian design.

	Classical	Bartlett-corrected	Adjusted
$\mathbb{P}\left\{\text{p-value}\le 5\%\right\}$	$11.1044\%\left(0.0668\%\right)$	$6.9592\%\left(0.0534\%\right)$	$5.0110\%\left(0.0453\%\right)$
$\mathbb{P}\left\{\text{p-value}\le 1\%\right\}$	$3.6383\%\left(0.038\%\right)$	$1.6975\%\left(0.0261\%\right)$	$0.9944\%\left(0.0186\%\right)$
$\mathbb{P}\left\{\text{p-value}\le 0.5\%\right\}$	$2.2477\%\left(0.0292\%\right)$	$0.9242\%\left(0.0178\%\right)$	$0.4952\%\left(0.0116\%\right)$
$\mathbb{P}\left\{\text{p-value}\le 0.1\%\right\}$	$0.7519\%\left(0.0155\%\right)$	$0.2306\%\left(0.0078\%\right)$	$0.1008\%\left(0.0051\%\right)$
$\mathbb{P}\left\{\text{p-value}\le 0.05\%\right\}$	$0.4669\%\left(0.0112\%\right)$	$0.124\%\left(0.0056\%\right)$	$0.0542\%\left(0.0036\%\right)$
$\mathbb{P}\left\{\text{p-value}\le 0.01\%\right\}$	$0.1575\%\left(0.0064\%\right)$	$0.0342\%\left(0.0027\%\right)$	$0.0104\%\left(0.0014\%\right)$

Table 2. Table 2 : Estimates of p-value probabilities with estimated Monte Carlo standard errors in parentheses under i.i.d. Bernoulli design.

	Adjusted
$\mathbb{P}\left\{\text{p-value}\le 5\%\right\}$	$5.0222\%\left(0.0412\%\right)$
$\mathbb{P}\left\{\text{p-value}\le 1\%\right\}$	$1.0048\%\left(0.0174\%\right)$
$\mathbb{P}\left\{\text{p-value}\le 0.5\%\right\}$	$0.5123\%\left(0.0119\%\right)$
$\mathbb{P}\left\{\text{p-value}\le 0.1\%\right\}$	$0.1108\%\left(0.005\%\right)$
$\mathbb{P}\left\{\text{p-value}\le 0.05\%\right\}$	$0.0521\%\left(0.0033\%\right)$
$\mathbb{P}\left\{\text{p-value}\le 0.01\%\right\}$	$0.0102\%\left(0.0015\%\right)$

Equations

$$E[y_i|{\mathbf{X}}_i]=\mu ({\mathbf{X}}_i^{\top }\mathbf{\beta }),$$

$$\mu (t):=\{\begin{array}{ll}{e}^t/(1+e^t)& \text{in the logistic model},\\ \Phi (t)& \text{in the probit model};\end{array}$$

$$\mu (t)+\mu (-t)=1$$

$$\ell \left(\mathbf{\beta }\right):=-\sum _{i=1}^n\left\{y_i\log \left(\frac{{\mu }_i}{1-{\mu }_i}\right)+\log \left(1-{\mu }_i\right)\right\},{\mu }_i:=\mu ({\mathbf{X}}_i^{\top }\mathbf{\beta }).$$

$$\ell\left(\bm{\beta}\right):=\sum\nolimits_{i=1}^{n}\rho\big{(}-\tilde{y}_{i}\bm{X}_{i}^{\top}\bm{\beta}\big{)},$$

$${\tilde{y}}_i:=\{\begin{array}{ll}1& \text{if }{y}_i=1,\\ -1& \text{if }{y}_i=0,\end{array}\text{and}\rho (t):=\{\begin{array}{ll}\log \left(1+e^t\right)& \text{in the logistic case},\\ -\log \Phi \left(-t\right)& \text{in the probit case}.\end{array}$$

$$H_j:{\beta }_j=0\text{versus}{\beta }_j\ne 0.$$

$$\Lambda_{j}:=\ell\big{(}\hat{\bm{\beta}}_{(-j)}\big{)}-\ell\big{(}\hat{\bm{\beta}}\big{)}.$$

$$\hat{\mathbf{\beta }}=\arg \underset{\mathbf{\beta }\in {\mathbb{R}}^p}{\min }\ell (\mathbf{\beta })\text{and}{\hat{\mathbf{\beta }}}_{(-j)}=\arg \underset{\mathbf{\beta }\in {\mathbb{R}}^p,{\beta }_j=0}{\min }\ell (\mathbf{\beta }).$$

$$2{\Lambda }_j(n)\text{leavevmode}\stackrel{\mathrm{d}}{\to }\text{leavevmode}{\chi }_1^2,\text{as }n\to \infty .$$

$$\mathbb{E}[2{\Lambda }_j]=1+\frac{\alpha }{n}+O\left(\frac{1}{n^2}\right),$$

$$\frac{2{\Lambda }_j}{1+\frac{{\alpha }_n}{n}}$$

$$E\left[\frac{2{\Lambda }_j}{1+\frac{{\alpha }_n}{n}}\right]=1+O\left(\frac{1}{n^2}\right)$$

$${\alpha }_n=\frac{n}{2}\left[\mathrm{Tr}\left({\mathbf{D}}_p^2\right)-\mathrm{Tr}\left({\mathbf{D}}_{p-1}^2\right)\right].$$

$$\underset{n\to \infty }{\lim }\frac{p(n)}{n}=\kappa ,$$

$$\Vert \hat{\mathbf{\beta }}\Vert =\infty \text{if }\kappa >1/2\text{and}\Vert \hat{\mathbf{\beta }}\Vert <\infty \text{if }\kappa <1/2.$$

$${\mathbf{X}}_i^{\top }\mathbf{\beta }>0\text{ if }\text{leavevmode}{\tilde{y}}_i=1\text{and}{\mathbf{X}}_i^{\top }\mathbf{\beta }<0\text{ if }\text{leavevmode}{\tilde{y}}_i=-1$$

$$\left\{\mathbf{X}\mathbf{\beta }\mid \mathbf{\beta }\in {\mathbb{R}}^p\right\}\cap {\mathbb{R}}_{++}^n\ne \varnothing ,$$

$$\left\{\mathbf{Z}\mathbf{\theta }\mid \mathbf{\theta }\in {\mathbb{R}}^p\right\}\cap {\mathbb{R}}_{++}^n\ne \varnothing .$$

$$\left\{\mathbf{Z}\mathbf{\theta }\mid \mathbf{\theta }\in {\mathbb{R}}^p\right\}\cap {\mathbb{R}}_{+}^n\ne \left\{0\right\}$$

$$\delta \left(\left\{\mathbf{Z}\mathbf{\theta }\mid \mathbf{\theta }\in {\mathbb{R}}^p\right\}\right)+\delta \left({\mathbb{R}}_{+}^n\right)>n+o(n).$$

$$\delta (\mathcal{K}):={\mathbb{E}}_{\mathbf{g}\sim \mathcal{N}(0,\mathbf{I})}\left[\Vert {\Pi }_{\mathcal{K}}\left(\mathbf{g}\right){\Vert }^2\right]$$

$$\delta \left(\left\{\mathbf{Z}\mathbf{\theta }\mid \mathbf{\theta }\in {\mathbb{R}}^p\right\}\right)=p\text{and}\delta ({\mathbb{R}}_{+}^n)=n/2,$$

$$p+n/2>n+o(n)\text{or}p/n>1/2+o(1),$$

$$\kappa <1/2.$$

$${\mathsf{prox}}_{b\rho }(z):=\arg \underset{x\in \mathbb{R}}{\min }\left\{b\rho (x)+\frac{1}{2}(x-z{)}^2\right\}$$

$$\Psi (z;b):=b{\rho }^{\prime }({\mathsf{prox}}_{b\rho }(z)),$$

$$\Psi (z;b)+{\mathsf{prox}}_{b\rho }(z)=z.$$

$${\mathsf{prox}}_f(z)+{\mathsf{prox}}_{f^{\ast }}(z)=z,\forall z,$$

$${\tau }^2$$

Full text

Assume that $n$ and $p$ grow large in such a way that $p/n\rightarrow\kappa$ for some constant $\kappa<1/2$. We prove that for a class of logistic models, the LLR converges to a rescaled chi-square, namely, $2\Lambda\leavevmode\nobreak\ \stackrel{{\scriptstyle\mathrm{d}}}{{\rightarrow}}\leavevmode\nobreak\ \alpha(\kappa)\chi_{k}^{2}$, where the scaling factor $\alpha(\kappa)$ is greater than one as soon as the dimensionality ratio $\kappa$ is positive. Hence, the LLR is larger than classically assumed. For instance, when $\kappa=0.3$, $\alpha(\kappa)\approx 1.5$. In general, we show how to compute the scaling factor by solving a nonlinear system of two equations with two unknowns. Our mathematical arguments are involved and use techniques from approximate message passing theory, from non-asymptotic random matrix theory and from convex geometry. We also complement our mathematical study by showing that the new limiting distribution is accurate for finite sample sizes.

Finally, all the results from this paper extend to some other regression models such as the probit regression model.

Keywords. Logistic regression, likelihood-ratio tests, Wilks’ theorem, high-dimensionality, goodness of fit, approximate message passing, concentration inequalities, convex geometry, leave-one-out analysis

1 Introduction

Logistic regression is by far the most widely used tool for relating a binary response to a family of explanatory variables. This model is used to infer the importance of variables and nearly all standard statistical softwares have inbuilt packages for obtaining p-values for assessing the significance of their coefficients. For instance, one can use the snippet of $\mathrm{R}$ code below to fit a logistic regression model from a vector $\verb|y|$ of binary responses and a matrix X of covariates:

fitted <- glm(y ~ X+0, family = ‘binomial’)
pvals  <- summary(fitted)$coefficients[,4]

The vector pvals stores p-values for testing whether a variable belongs to a model or not, and it is well known that the underlying calculations used to produce these p-values can also be used to construct confidence intervals for the regression coefficients. Since logistic models are used hundreds of times every day for inference purposes, it is important to know whether these calculations—e.g. these p-values—are accurate and can be trusted.

1.1 Binary regression

Imagine we have $n$ samples of the form $(y_{i},\bm{X}_{i})$, where $y_{i}\in\{0,1\}$ and $\bm{X}_{i}\in\mathbb{R}^{p}$. In a generalized linear model, one postulates the existence of a link function $\mu(\cdot)$ relating the conditional mean of the response variable to the linear predictor $\bm{X}_{i}^{\top}\bm{\beta}$,

[TABLE]

where $\bm{\beta}=[\beta_{1},\beta_{2},\ldots,\beta_{p}]^{\top}\in\mathbb{R}^{p}$ is an unknown vector of parameters. We focus here on the two most commonly used binary regression models, namely, the logistic and the probit models for which

[TABLE]

here, $\Phi$ is the cumulative distribution function (CDF) of a standard normal random variable. In both cases, the Symmetry Condition

[TABLE]

holds, which says that the two types $y_{i}=0$ and $y_{i}=1$ are treated in a symmetric fashion. Assuming that the observations are independent, the negative log-likelihood function is given by [1, Section 4.1.2]

[TABLE]

Invoking the symmetry condition, a little algebra reveals an equivalent expression

[TABLE]

where

[TABLE]

Throughout we refer to this function $\rho$ as the effective link.

1.2 The likelihood-ratio test and Wilks’ phenomenon

Researchers often wish to determine which covariates are of importance, or more precisely, to test whether the $j{\text{th}}$ variable belongs to the model or not: formally, we wish to test the hypothesis

[TABLE]

Arguably, one of the most commonly deployed techniques for testing $H_{j}$ is the likelihood-ratio test (LRT), which is based on the log-likelihood ratio (LLR) statistic

[TABLE]

Here, $\hat{\bm{\beta}}$ and $\hat{\bm{\beta}}_{(-j)}$ denote respectively the maximum likelihood estimates (MLEs) under the full model and the reduced model on dropping the $j{\text{th}}$ predictor; that is,

[TABLE]

Inference based on such log-likelihood ratio statistics has been studied extensively in prior literature [52, 12, 37]. Arguably, one of the most celebrated results in the large-sample regime is the Wilks’ theorem.

To describe the Wilk’s phenomenon, imagine we have a sequence of observations $(y_{i},\bm{X}_{i})$ where $y_{i}\in\{0,1\}$, $\bm{X}_{i}\in\mathbb{R}^{p}$ with $p$ fixed. Since we are interested in the limit of large samples, we may want to assume that the covariates are i.i.d. drawn from some population with non-degenerate covariance matrix so that the problem is fully $p$-dimensional. As before, we assume a conditional logistic model for the response. In this setting, Wilks’ theorem [52] calculates the asymptotic distribution of $\Lambda_{j}(n)$ when $n$ grows to infinity:

• (Wilks’ phenomenon) Under suitable regularity conditions which, for instance, guarantee that the MLE exists and is unique,111Such conditions would also typically imply asymptotic normality of the MLE. the LLR statistic for testing $H_{j}:\beta_{j}=0$ vs. $\beta_{j}\neq 0$ has asymptotic distribution under the null given by

[TABLE]

This fixed-$p$ large-$n$ asymptotic result, which is a consequence of asymptotic normality properties of the MLE[50, Theorem 5.14], applies to a much broader class of testing problems in parametric models; for instance, it applies to the probit model as well. We refer the readers to [34, Chapter 12] and [50, Chapter 16] for a thorough exposition and details on the regularity conditions under which Wilks’ theorem holds. Finally, there is a well-known extension which states that if we were to drop $k$ variables from the model, then the LLR would converge to a chi-square distribution with $k$ degrees of freedom under the hypothesis that the reduced model is correct.

1.3 Inadequacy of Wilks’ theorem in high dimensions

The chi-square approximation to the distribution of the LLR statistic is used in standard statistical softwares to provide p-values for the single or multiple coefficient likelihood ratio tests. Here, we perform a simple experiment on synthetic data to study the accuracy of the chi-square approximation when $p$ and $n$ are both decently large. Specifically, we set $\bm{\beta}=\bm{0}$ and test $\beta_{1}=0$ vs. $\beta_{1}\neq 0$ using the LRT in a setting where $p=1200$. In each trial, $n=4000$ observations are produced with $y_{i}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}{\mathsf{Bernoulli}}(1/2)$, and $\bm{X}:=[\bm{X}_{1},\cdots,\bm{X}_{n}]^{\top}\in\mathbb{R}^{n\times p}$ is obtained by generating a random matrix composed of i.i.d. $\mathcal{N}(0,1)$ entries. We fit a logistic regression of $\bm{y}$ on $\bm{X}$ using R, and extract the p-values for each coefficient. Figure 1 plots the pooled histogram that aggregates $4.8\times 10^{5}$ p-values in total (400 trials with $1200$ p-values obtained in each trial).

If the $\chi^{2}_{1}$ approximation were true, then we would expect to observe uniformly distributed p-values. The histrogram from Fig. 1 is, however, far from uniform. This is an indication of the inadequacy of Wilks’ theorem when $p$ and $n$ are both large. The same issue was also reported in [11], where the authors observed that this discrepancy is highly problematic since the distribution is skewed towards smaller values. Hence, such p-values cannot be trusted to construct level-$\alpha$ tests and the problem is increasingly severe when we turn attention to smaller p-values as in large-scale multiple testing applications.

1.4 The Bartlett correction?

A natural question that arises immediately is whether the observed discrepancy could be an outcome of a finite-sample effect. It has been repeatedly observed that the chi-square approximation does not yield accurate results with finite sample size. One correction to the LRT that is widely used in finite samples is the Bartlett correction, which dates back to Bartlett [5] and has been extensively studied over the past few decades (e.g. [10, 33, 8, 15, 13]). In the context of testing for a single coefficient in the logistic model, this correction can be described as follows [38]: compute the expectation of the LLR statistic up to terms of order $1/n^{2}$; that is, compute a parameter $\alpha$ such that

[TABLE]

which suggests a corrected LLR statistic

[TABLE]

with $\alpha_{n}$ being an estimator of $\alpha$. With a proper choice of $\alpha_{n}$, one can ensure

[TABLE]

in the classical setting where $p$ is fixed and $n$ diverges. In expectation, this corrected statistic is closer to a $\chi_{1}^{2}$ distribution than the original LLR for finite samples. Notably, the correction factor may in general be a function of the unknown $\bm{\beta}$ and, in that case, must be estimated from the null model via maximum likelihood estimation.

In the context of GLMs, Cordeiro [13] derived a general formula for the Bartlett corrected LLR statistic, see [17, 14] for a detailed survey. In the case where there is no signal ($\bm{\beta}=\bm{0}$), one can compute $\alpha_{n}$ for the logistic regression model following [13] and[38], which yields

[TABLE]

Here, $\bm{D}_{p}$ is the diagonal part of $\bm{X}(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}$ and $\bm{D}_{p-1}$ is that of \bm{X}_{(-j)}\big{(}\bm{X}_{(-j)}^{\top}\bm{X}_{(-j)}\big{)}^{-1}\bm{X}_{(-j)}^{\top} in which $\bm{X}_{(-j)}$ is the design matrix $\bm{X}$ with the $j{\text{th}}$ column removed. Comparing the adjusted LLRs to a $\chi^{2}_{1}$ distribution yields adjusted p-values. In the setting of Fig. 1(a), the histogram of Bartlett corrected p-values is shown in Fig. 1(b). As we see, these p-values are still far from uniform.

If the mismatch is not due to finite sample-size effects, what is the distribution of the LLR in high dimensions? Our main contribution is to provide a very precise answer to this question; below, we derive the high-dimensional asymptotic distribution of the log-likelihood ratios, i.e. in situations where the dimension $p$ is not necessarily negligible compared to the sample size $n$.

2 Main results

2.1 Modelling assumptions

In this paper, we focus on the high-dimensional regime where the sample size is not much larger than the number of parameters to be estimated—a setting which has attracted a flurry of activity in recent years. In particular, we assume that the number $p(n)$ of covariates grows proportionally with the number $n$ of observations; that is,

[TABLE]

where $\kappa>0$ is a fixed constant independent of $n$ and $p(n)$. In fact, we shall also assume $\kappa<1/2$ for both the logistic and the probit models, as the MLE is otherwise at $\infty$; see Section 2.2.

To formalize the notion of high-dimensional asymptotics when both $n$ and $p(n)$ diverge, we consider a sequence of instances $\{\bm{X}(n),\bm{y}(n)\}_{n\geq 0}$ such that for any $n$,

• •

$\bm{X}(n)\in\mathbb{R}^{n\times p(n)}$ has i.i.d. rows $\bm{X}_{i}(n)\sim{\mathcal{N}}(0,\bm{\Sigma})$, where $\bm{\Sigma}\in\mathbb{R}^{p(n)\times p(n)}$ is positive definite;

• •

$y_{i}(n)\,|\,\bm{X}(n)\sim y_{i}(n)\,|\,\bm{X}_{i}(n)\stackrel{{\scriptstyle\mathrm{ind.}}}{{\sim}}{\mathsf{Bernoulli}}\left(\mu(\bm{X}_{i}(n)^{\top}\bm{\beta}(n))\right)$, where $\mu$ satisfies the Symmetry Condition;

• •

we further assume $\bm{\beta}(n)=\bm{0}$. From the Symmetry Condition it follows that $\mu(0)=1/2$, which directly implies that $\bm{y}(n)$ is a vector with i.i.d ${\mathsf{Bernoulli}}(1/2)$ entries.

The MLE is denoted by $\hat{\bm{\beta}}(n)$ and there are $p(n)$ LLR statistics $\Lambda_{j}(n)$ ($1\leq j\leq p(n)$), one for each of the $p(n)$ regression coefficients. In the sequel, the dependency on $n$ shall be suppressed whenever it is clear from the context.

2.2 When does the MLE exist?

Even though we are operating in the regime where $n>p$, the existence of the MLE cannot be guaranteed for all $p$ and $n$. Interestingly, the norm of the MLE undergoes a sharp phase transition in the sense that

[TABLE]

Here, we develop some understanding about this phenomenon. Given that $\rho(t)\geq\rho(-\infty)=0$ for both the logistic and probit models, each summand in (4) is minimized if $\tilde{y}_{i}\bm{X}_{i}^{\top}\bm{\beta}=\infty$, which occurs when $\text{sign}(\bm{X}_{i}^{\top}\bm{\beta})=\text{sign}(\tilde{y}_{i})$ and $\|\bm{\beta}\|=\infty$. As a result, if there exists a nontrivial ray $\bm{\beta}$ such that

[TABLE]

for any $1\leq i\leq n$, then pushing $\|\bm{\beta}\|$ to infinity leads to an optimizer of (4). In other words, the solution to (4) becomes unbounded (the MLE is at $\infty$) whenever there is a hyperplane perfectly separating the two sets of samples $\{i\mid\tilde{y}_{i}=1\}$ and $\{i\mid\tilde{y}_{i}=-1\}$.

Under the assumptions from Section 2.1, $\tilde{y}_{i}$ is independent of $\bm{X}$ and the distribution of $\bm{X}$ is symmetric. Hence, to calculate the chance that there exists a separating hyperplane, we can assume $\tilde{y}_{i}=1$ $(1\leq i\leq n)$ without loss of generality. In this case, the event (12) becomes

[TABLE]

where $\mathbb{R}_{++}^{n}$ is the positive orthant. Write $\bm{X}=\bm{Z}\bm{\Sigma}^{1/2}$ so that $\bm{Z}$ is an $n\times p$ matrix with i.i.d. standard Gaussian entries, and $\bm{\theta}=\bm{\Sigma}^{1/2}\bm{\beta}$. Then the event (13) is equivalent to

[TABLE]

Now the probability that (14) occurs is the same as that

[TABLE]

occurs, where $\mathbb{R}_{+}^{n}$ denotes the non-negative orthant. From the approximate kinematic formula [3, Theorem I] in the literature on convex geometry, the event (15) happens with high probability if and only if the total statistical dimension of the two closed convex cones exceeds the ambient dimension, i.e.

[TABLE]

Here, the statistical dimension of a closed convex cone $\mathcal{K}$ is defined as

[TABLE]

with $\Pi_{\mathcal{K}}\left(\bm{g}\right):=\arg\min_{\bm{z}\in\mathcal{K}}\|\bm{g}-\bm{z}\|$ the Euclidean projection. Recognizing that[3, Proposition 2.4]

[TABLE]

we reduce the condition (16) to

[TABLE]

thus indicating that $\|\hat{\bm{\beta}}\|=\infty$ when $\kappa=\lim p/n>1/2$.

This argument only reveals that $\|\hat{\bm{\beta}}\|=\infty$ in the regime where $\kappa>1/2$. If $\kappa=p/n<1/2$, then $\|\bm{\Sigma}^{1/2}\hat{\bm{\beta}}\|=O(1)$ with high probability, a fact we shall prove in Section 5. In light of these observations we work with the additional condition

[TABLE]

2.3 The high-dimensional limiting distribution of the LLR

In contrast to the classical Wilks’ result, our findings reveal that the LLR statistic follows a rescaled chi-square distribution with a rescaling factor that can be explicitly pinned down through the solution to a system of equations.

2.3.1 A system of equations

We start by setting up the crucial system of equations. Before proceeding, we first recall the proximal operator

[TABLE]

defined for any $b>0$ and convex function $\rho(\cdot)$. As in [18], we introduce the operator

[TABLE]

which is simply the proximal operator of the conjugate $(b\rho)^{*}$ of $b\rho$.222The conjugate $f^{*}$ of a function $f$ is defined as $f^{*}(x)=\sup_{u\in\mathrm{dom}(f)}\{\langle u,x\rangle-f(u)\}$.

To see this, we note that $\Psi$ satisfies the relation [18, Proposition 6.4]

[TABLE]

The claim that $\Psi(\cdot;b)=\mathsf{prox}_{(b\rho)^{*}}(\cdot)$ then follows from the Moreau decomposition

[TABLE]

which holds for a closed convex function $f$ [39, Section 2.5]. Interested readers are referred to [18, Appendix 1] for more properties of $\mathsf{prox}_{b\rho}$ and $\Psi$.

We are now in position to present the system of equations that plays a crucial role in determining the distribution of the LLR statistic in high dimensions:

[TABLE]

where $Z\sim{\mathcal{N}}(0,1)$, and $\Psi^{\prime}\left(\cdot,\cdot\right)$ denotes differentiation with respect to the first variable. The fact that this system of equations would admit a unique solution in $\mathbb{R}_{+}^{2}$ is not obvious a priori. We shall establish this for the logistic and the probit models later in Section 6.

2.3.2 Main result

Theorem 1**.**

Consider a logistic or probit regression model under the assumptions from Section 2.1. If $\kappa\in(0,1/2)$, then for any $1\leq j\leq p$, the log-likelihood ratio statistic $\Lambda_{j}$ as defined in (7) obeys

[TABLE]

*where $(\tau_{\ast},b_{\ast})\in\mathbb{R}_{+}^{2}$ is the unique solution to the system of equations (23) and (24). Furthermore, the LLR statistic obtained by dropping $k$ variables for any fixed $k$ converges to $({\tau_{\ast}^{2}}/{b_{\ast}})\chi_{k}^{2}$. Finally, these results extend to all binary regression models with links obeying the assumptions listed in Section 2.3.3. *

Hence, the limiting distribution is a rescaled chi-square with a rescaling factor ${\tau_{\ast}^{2}}/{b_{\ast}}$ that only depends on the aspect ratio $\kappa$. Fig. 2 illustrates the dependence of the rescaling factor on the limiting aspect ratio $\kappa$ for logistic regression. The figures for the probit model are similar as the rescaling constants actually differ by very small values.

To study the quality of approximation for finite samples, we repeat the same numerical experiments as before but now obtain the p-values by comparing the LLR statistic with the rescaled chi-square suggested by Theorem 1. For a particular run of the experiment ($n=4000,p=1200,\kappa=0.3$), we compute the adjusted LLR statistic $({2b_{\ast}}/{\tau_{\ast}^{2}})\Lambda_{j}$ for each coefficient and obtain the p-values based on the $\chi^{2}_{1}$ distribution. The pooled histogram that aggregates $4.8\times 10^{5}$ p-values in total is shown in Fig. 1(c).

As we clearly see, the p-values are much closer to a uniform distribution now. One can compute the chi-square goodness of fit statistic to test the closeness of the above distribution to uniformity. To this end, we divide the interval $[0,1]$ into $20$ equally spaced bins of width $0.05$ each. For each bin we compute the observed number of times a p-value falls in the bin out of the $4.8\times 10^{5}$ values. Then a chi-square goodness of fit statistic is computed, noting that the expected frequency is $24000$ for each bin. The chi-square statistic in this case is $16.049$, which gives a p-value of $0.654$ in comparison with a $\chi^{2}_{19}$ variable. The same test when performed with the Bartlett corrected p-values (Fig. 1(b)) yields a chi-square statistic $5599$ with a p-value of [math]. 333 Note that the p-values obtained at each trial are not exactly independent. However, they are exchangeable, and weakly dependent (see the proof of Corollary 1 for a formal justification of this fact). Therefore, we expect the goodness of fit test to be an approximately valid procedure in this setting. Thus, our correction gives the desired uniformity in the p-values when the true signal $\bm{\beta}=\bm{0}$.

Practitioners would be concerned about the validity of p-values when they are small—again, think about multiple testing applications. In order to study whether our correction yields valid results for small p-values, we compute the proportion of times the p-values (in all the three cases) lie below $5\%,1\%,0.5\%,0.1\%$ out of the $4.8\times 10^{5}$ times. The results are summarized in Table 1. This further illustrates the deviation from uniformity for the classical and Bartlett corrected p-values, whereas the “adjusted” p-values obtained invoking Theorem 1 are still valid.

2.3.3 Extensions

As noted in Section 1.1, the Symmetry Condition (3) allows to express the negative log-likelihood in the form (4), which makes use of the effective link $\rho(\cdot)$. Theorem 1 applies to any $\rho(\cdot)$ obeying the following properties:

$\rho$ is non-negative, has up to three derivatives, and obeys $\rho(t)\geq t$. 2. 2.

$\rho^{\prime}$ may be unbounded but it should grow sufficiently slowly, in particular, we assume $|\rho^{\prime}(t)|=O(|t|)$ and $\rho^{\prime}(\mathsf{prox}_{c\rho}(Z))$ is a sub-Gaussian random variable for any constant $c>0$ and any $Z\sim{\mathcal{N}}(0,\sigma^{2})$ for some finite $\sigma>0$. 3. 3.

$\rho^{\prime\prime}(t)>0$ for any $t$ which implies that $\rho$ is convex, and $\sup_{t}\rho^{\prime\prime}(t)<\infty$. 4. 4.

$\sup_{t}|\rho^{\prime\prime\prime}(t)|<\infty$. 5. 5.

Given any $\tau>0$ ,the equation (24) has a unique solution in $b$. 6. 6.

The map $\mathcal{V}(\tau^{2})$ as defined in (58) has a fixed point.

It can be checked that the effective links for both the logistic and the probit models (5) obey all of the above. The last two conditions are assumed to ensure existence of a unique solution to the system of equations (23) and (24) as will be seen in Section 6; we shall justify these two conditions for the logistic and the probit models in Section 6.1.

2.4 Reduction to independent covariates

In order to derive the asymptotic distribution of the LLR statistics, it in fact suffices to consider the special case $\bm{\Sigma}=\bm{I}_{p}$.

Lemma 1**.**

Let $\Lambda_{j}(\bm{X})$ be the LLR statistic based on the design matrix $\bm{X}$, where the rows of $\bm{X}$ are i.i.d. ${\mathcal{N}}({\bm{0}},\bm{\Sigma})$ and $\Lambda_{j}(\bm{Z})$ that where the rows are i.i.d. ${\mathcal{N}}({\bm{0}},\bm{I}_{p})$. Then

[TABLE]

• Proof:

Recall from (4) that the LLR statistic for testing the $j$th coefficient can be expressed as

[TABLE]

Write $\bm{Z}^{\prime}=\bm{X}\bm{\Sigma}^{-1/2}$ so that the rows of $\bm{Z}^{\prime}$ are i.i.d. ${\mathcal{N}}({\bm{0}},\bm{I}_{p})$ and set $\bm{\theta}^{\prime}=\bm{\Sigma}^{1/2}\bm{\beta}$. With this reparameterization, we observe that the constraint $\beta_{j}=0$ is equivalent to $\bm{a}_{j}^{\top}\bm{\theta}^{\prime}=0$ for some non-zero vector $\bm{a}_{j}\in\mathbb{R}^{p}$. This gives

[TABLE]

Now let $\bm{Q}$ be an orthogonal matrix mapping $\bm{a}_{j}\in\mathbb{R}^{p}$ into the vector $\|\bm{a}_{j}\|\bm{e}_{j}\in\mathbb{R}^{p}$, i.e. $\bm{Q}\bm{a}_{j}=\|\bm{a}_{j}\|\bm{e}_{j}$. Additionally, set $\bm{Z}=\bm{Z}^{\prime}\bm{Q}$ (the rows of $\bm{Z}$ are still i.i.d. ${\mathcal{N}}({\bm{0}},\bm{I}_{p})$) and $\bm{\theta}=\bm{Q}\bm{\theta}^{\prime}$. Since $\bm{a}_{j}^{\top}\bm{\theta}^{\prime}=0$ occurs if and only if $\theta_{j}=0$, we obtain

[TABLE]

which proves the lemma.

∎

In the remainder of the paper we, therefore, assume $\bm{\Sigma}=\bm{I}_{p}$.

2.5 Proof architecture

This section presents the main steps for proving Theorem 1. We will only prove the theorem for $\{\Lambda_{j}\}$, the LLR statistic obtained by dropping a single variable. The analysis for the LLR statistic obtained on dropping $k$ variables (for some fixed $k$) follows very similar steps and is hence omitted for the sake of conciceness. As discussed before, we are free to work with any configuration of the $y_{i}$’s. For the two steps below, we will adopt two different configurations for convenience of presentation.

2.5.1 Step 1: characterizing the asymptotic distributions of $\hat{\beta}_{j}$

Without loss of generality, we assume here that $y_{i}=1$ (and hence $\tilde{y}_{i}=1$) for all $1\leq i\leq n$ and, therefore, the MLE problem reduces to

[TABLE]

We would first like to characterize the marginal distribution of $\hat{\bm{\beta}}$, which is crucial in understanding the LLR statistic. To this end, our analysis follows by a reduction to the setup of [21, 18, 19, 20], with certain modifications that are called for due to the specific choices of $\rho(\cdot)$ we deal with here. Specifically, consider the linear model

[TABLE]

and prior work [21, 18, 19, 20] investigating the associated M-estimator

[TABLE]

Our problem reduces to (27) on setting $\bm{y}=\bm{w}=\bm{0}$ in (27). When $\rho(\cdot)$ satisfies certain assumptions (e.g. strong convexity), the asymptotic distribution of $\|\hat{\bm{\beta}}\|$ has been studied in a series of works [21, 19, 20] using a leave-one-out analysis and independently in [18] using approximate message passing (AMP) machinery. An outline of their main results is described in Section 2.7. However, the function $\rho(\cdot)$ in our cases has vanishing curvature and, therefore, lacks the essential strong convexity assumption that was utilized in both the aforementioned lines of work. To circumvent this issue, we propose to invoke the AMP machinery as in [18], in conjunction with the following critical additional ingredients:

• •

(Norm Bound Condition) We utilize results from the conic geometry literature (e.g. [3]) to establish that

[TABLE]

with high probability as long as $\kappa<1/2$. This will be elaborated in Theorem 4.

• •

(Likelihood Curvature Condition) We establish some regularity conditions on the Hessian of the log-likelihood function, generalizing the strong convexity condition, which will be detailed in Lemma 4.

• •

(Uniqueness of the Solution to (23) and (24)) We establish that for both the logistic and the probit case, the system of equations (23) and (24) admits a unique solution.

We emphasize that these elements are not straightforward, require significant effort and a number of novel ideas, which form our primary technical contributions for this step.

These ingredients enable the use of the AMP machinery even in the absence of strong convexity on $\rho(\cdot)$, finally leading to the following theorem:

Theorem 2**.**

Under the conditions of Theorem 1,

[TABLE]

This theorem immediately implies that the marginal distribution of $\hat{{\beta}}_{j}$ is normal.

Corollary 1**.**

Under the conditions of Theorem 1, for every $1\leq j\leq p$, it holds that

[TABLE]

• Proof:

From the rotational invariance of our i.i.d. Gaussian design, it can be easily verified that $\hat{\bm{\beta}}/\|\hat{\bm{\beta}}\|$ is uniformly distributed on the unit sphere $\mathbb{S}^{p-1}$ and is independent of $\|\hat{\bm{\beta}}\|$. Therefore, $\hat{\beta}_{j}$ has the same distribution as $\|\hat{\bm{\beta}}\|Z_{j}/\|\bm{Z}\|$, where $\bm{Z}=(Z_{1},\ldots,Z_{p})\sim{\mathcal{N}}(\bm{0},\bm{I}_{p})$ independent of $\|\hat{\bm{\beta}}\|$. Since $\sqrt{p}\|\hat{\bm{\beta}}\|/\|\bm{Z}\|$ converges in probability to $\tau_{\ast}$, we have, by Slutsky’s theorem, that $\sqrt{p}\hat{\beta}_{j}$ converges to ${\mathcal{N}}(0,\tau_{\ast}^{2})$ in distribution.

∎

2.5.2 Step 2: connecting $\Lambda_{j}$ with $\hat{\beta}_{j}$

Now that we have derived the asymptotic distribution of $\hat{\beta}_{j}$, the next step involves a reduction of the LLR statistic to a function of the relevant coordinate of the MLE. Before continuing, we note that the distribution of $\Lambda_{j}$ is the same for all $1\leq j\leq p$ due to exchangeability. As a result, going forward we will only analyze $\Lambda_{1}$ without loss of generality. In addition, we introduce the following convenient notations and assumptions:

• •

the design matrix on dropping the first column is written as $\tilde{\bm{X}}$ and the MLE in the corresponding reduced model as $\tilde{\bm{\beta}}$;

• •

write $\bm{X}=[\bm{X}_{1},\cdots,\bm{X}_{n}]^{\top}\in\mathbb{R}^{n\times p}$ and $\tilde{\bm{X}}=[\tilde{\bm{X}}_{1},\cdots,\tilde{\bm{X}}_{n}]^{\top}\in\mathbb{R}^{n\times(p-1)}$;

• •

without loss of generality, assume that $\tilde{y}_{i}=-1$ for all $i$ in this subsection, and hence the MLEs under the full and the reduced models reduce to

[TABLE]

With the above notations in place, the LLR statistic for testing $\beta_{1}=0$ vs. $\beta_{1}\neq 0$ can be expressed as

[TABLE]

To analyze $\Lambda_{1}$, we invoke Taylor expansion to reach

[TABLE]

where $\gamma_{i}$ lies between $\tilde{\bm{X}}^{\top}_{i}\tilde{\bm{\beta}}$ and $\bm{X}^{\top}_{i}\hat{\bm{\beta}}$. A key observation is that the linear term $Q_{\mathrm{lin}}$ in the above equation vanishes. To see this, note that the first-order optimality conditions for the MLE $\hat{\bm{\beta}}$ is given by

[TABLE]

Replacing $\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{\beta}}$ with \bm{X}_{i}^{\top}\left[\begin{array}[]{c}0\\ \tilde{\bm{\beta}}\end{array}\right] in $Q_{\text{lin}}$ and using the optimality condition, we obtain

[TABLE]

Consequently, $\Lambda_{1}$ simplifies to the following form

[TABLE]

Thus, computing the asymptotic distribution of $\Lambda_{1}$ boils down to analyzing $\bm{X}_{i}^{\top}\hat{\bm{\beta}}-\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{\beta}}$. Our argument is inspired by the leave-one-predictor-out approach developed in [19, 20].

We re-emphasize that our setting is not covered by that of [19, 20], due to the violation of strong convexity and some other technical assumptions. We sidestep this issue by utilizing the Norm Bound Condition and the Likelihood Curvature Condition. In the end, our analysis establishes the equivalence of $\Lambda_{1}$ and $\hat{\beta}_{1}$ up to some explicit multiplicative factors modulo negligible error terms. This is summarized as follows.

Theorem 3**.**

Under the assumptions of Theorem 1,

[TABLE]

Theorem 3 reveals a simple yet surprising connection between the LLR statistic $\Lambda_{1}$ and the MLE $\hat{\bm{\beta}}$. As we shall see in the proof of the theorem, the quadratic term in (36) is $\frac{1}{2}\frac{p}{b_{\ast}}\hat{\beta}_{1}^{2}+o(1)$, while the remaining third-order term of (36) is vanishingly small. Finally, putting Corollary 1 and Theorem 3 together directly establishes Theorem 1.

2.6 Comparisons with the classical regime

We pause to shed some light on the interpretation of the correction factor $\tau_{\ast}^{2}/b_{\ast}$ in Theorem 1 and understand the differences from classical results. Classical theory (e.g. [30, 31]) asserts that when $p$ is fixed and $n$ diverges, the MLE for a fixed design $\bm{X}$ is asymptotically normal, namely,

[TABLE]

where

[TABLE]

is the normalized Fisher information at the true value $\bm{\beta}$. In particular, under the global null and i.i.d. Gaussian design, this converges to

[TABLE]

as $n$ tends to infinity [50, Example 5.40].

The behavior in high dimensions is different. In particular, Corollary 1 states that under the global null, we have

[TABLE]

Comparing the variances in the logistic model, we have that

[TABLE]

Fig. 3 illustrates the behavior of the ratio $\tau_{\ast}^{2}/\kappa$ as a function of $\kappa$. Two observations are immediate:

• •

First, in Fig. 3(a) we have $\tau_{\ast}^{2}\geq 4\kappa$ for all $\kappa\geq 0$. This indicates an inflation in variance or an “extra Gaussian noise” component that appears in high dimensions, as discussed in [18]. The variance of the “extra Gaussian noise” component increases as $\kappa$ grows.

• •

Second, as $\kappa\rightarrow 0$, we have $\tau_{\ast}^{2}/4\kappa\rightarrow 1$ in the logistic model, which indicates that classical theory becomes accurate in this case. In other words, our theory recovers the classical prediction in the regime where $p=o(n)$.

Further, for the testing problem considered here, the LLR statistic in the classical setup can be expressed, through Taylor expansion, as

[TABLE]

where $\tilde{\bm{\beta}}$ is defined in (31). In the high-dimensional setting, we will also establish a quadratic approximation of the form

[TABLE]

In Theorem 7, we shall see that $b_{\ast}$ is the limit of $\frac{1}{n}\mathrm{Tr}(\bm{G}^{-1})$, the Stieltjes transform of the empirical spectral distribution of $\bm{G}$ evaluated at [math]. Thus, this quantity in some sense captures the spread in the eigenvalues of $\bm{G}$ one would expect to happen in high dimensions.

2.7 Prior art

Wilks’ type of phenomenon in the presence of a diverging dimension $p$ has received much attention in the past. For instance, Portnoy [43] investigated simple hypotheses in regular exponential families, and established the asymptotic chi-square approximation for the LLR test statistic as long as $p^{3/2}/n\rightarrow 0$. This phenomenon was later extended in [46] to accommodate the MLE with a quadratic penalization, and in [53] to account for parametric models underlying several random graph models. Going beyond parametric inference, Fan et al. [24, 22] explored extensions to infinite-dimensional non-parametric inference problems, for which the MLE might not even exist or might be difficult to derive. While the classical Wilks’ phenomenon fails to hold in such settings, Fan et al. [24, 22] proposed a generalization of the likelihood ratio statistics based on suitable non-parametric estimators and characterized the asymptotic distributions. Such results have further motivated Boucheron and Massart [9] to investigate the non-asymptotic Wilks’ phenomenon or, more precisely, the concentration behavior of the difference between the excess empirical risk and the true risk, from a statistical learning theory perspective. The Wilks’ phenomenon for penalized empirical likelihood has also been established [48]. However, the precise asymptotic behavior of the LLR statistic in the regime that permits $p$ to grow proportional to $n$ is still beyond reach.

On the other hand, as demonstrated in Section 2.5.1, the MLE here under the global null can be viewed as an M-estimator for a linear regression problem. Questions regarding the behavior of robust linear regression estimators in high dimensions—where $p$ is allowed to grow with $n$—–were raised in Huber [30], and have been extensively studied in subsequent works, e.g. [40, 41, 42, 36]. When it comes to logistic regression, the behavior of the MLE was studied for a diverging number of parameters by [28], which characterized the squared estimation error of the MLE if $(p\log p)/n\rightarrow 0$. In addition, the asymptotic normality properties of the MLE and the penalized MLE for logistic regression have been established by [35] and [23], respectively. A very recent paper by Fan et al. [25] studied the logistic model under the global null $\bm{\beta}=\bm{0}$, and investigated the classical asymptotic normality as given in (38). It was discovered in [25] that the convergence property (38) breaks down even in terms of the marginal distribution, namely,

[TABLE]

as soon as $p$ grows at a rate exceeding $n^{2/3}$. This result, however, does not imply the asymptotic distribution of the likelihood-ratio statistic in this regime. In fact, our theorem implies that LLR statistic $2\Lambda_{j}$ goes to $\chi_{1}^{2}$ (and hence Wilks phenomenon remains valid) when $\kappa=p/n\rightarrow 0$.

The line of work that is most relevant to the present paper was initially started by El Karoui et al. [21]. Focusing on the regime where $p$ is comparable to $n$, the authors uncovered, via a non-rigorous argument, that the asymptotic $\ell_{2}$ error of the MLE could be characterized by a system of nonlinear equations. This seminal result was later made rigorous independently by Donoho et al. [18] under i.i.d. Gaussian design and by El Karoui [19, 20] under more general i.i.d. random design as well as certain assumptions on the error distribution. Both approaches rely on strong convexity on the function $\rho(\cdot)$ that defines the M-estimator, which does not hold in the models considered herein.

2.8 Notations

We adopt the standard notation $f(n)=O\left(g(n)\right)$ or $f(n)\lesssim g(n)$ which means that there exists a constant $c>0$ such that $\left|f(n)\right|\leq c|g(n)|$. Likewise, $f(n)=\Omega\left(g(n)\right)$ or $f(n)\gtrsim g(n)$ means that there exists a constant $c>0$ such that $|f(n)|\geq c\left|g(n)\right|$, $f(n)\asymp g(n)$ means that there exist constants $c_{1},c_{2}>0$ such that $c_{1}|g(n)|\leq|f(n)|\leq c_{2}|g(n)|$, and $f(n)=o(g(n))$ means that $\lim_{n\rightarrow\infty}\frac{f(n)}{g(n)}=0$. Any mention of $C$, $C_{i}$, $c$, $c_{i}$ for $i\in\mathbb{N}$ refers to some positive universal constants whose value may change from line to line. For a square symmetric matrix $\bm{M}$, the minimum eigenvalue is denoted by $\lambda_{\min}(\bm{M})$. Logarithms are base $e$.

3 Numerics

3.1 Non-Gaussian covariates

In this section we first study the sensitivity of our result to the Gaussianity assumption on the design matrix. To this end, we consider a high dimensional binary regression set up with a Bernoulli design matrix. We simulate $n=4000$ i.i.d. observations $(y_{i},\bm{X}_{i})$ with $y_{i}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}{\mathsf{Bernoulli}}(1/2)$, and $\bm{X}_{i}$ generated independent of $y_{i}$, such that each entry takes on values in $\{1,-1\}$ w.p. $1/2$. At each trial, we fit a logistic regression model to the data and obtain the classical, Bartlett corrected and adjusted p-values (using the rescaling factor $\tau_{\ast}^{2}/b_{\ast}$). Figure 4 plots the histograms for the pooled p-values, obtained across $400$ trials.

It is instructive to compare the histograms to that obtained in the Gaussian case (Figure 1). The classical and Bartlett corrected p-values exhibit similar deviations from uniformity as in the Gaussian design case, whereas our adjusted p-values continue to have an approximate uniform distribution. We test for deviations from uniformity using a formal chi-squared goodness of fit test as in Section 2.3.2. For the Bartlett corrected p-values, the chi-squared statistic turns out to be $5885$, with a p-value [math]. For the adjusted p-values,the chi-squared statistic is $24.1024$, with a p-value $0.1922$.444Recall our earlier footnote about the use of a $\chi^{2}$ test.

Once again, the Bartlett correction fails to provide valid p-values whereas the adjusted p-values are consistent with a uniform distribution. These findings indicate that the distribution of the LLR statistic under the i.i.d. Bernoulli design is in agreement to the rescaled $\chi^{2}_{1}$ derived under the Gaussian design in Theorem 1, suggesting that the distribution is not too sensitive to the Gaussianity assumption. Estimates of p-value probabilities for our method are provided in Table 2.

3.2 Quality of approximations for finite sample sizes

In the rest of this section, we report some numerical experiments which study the applicability of our theory in finite sample setups.

Validity of tail approximation The first experiment explores the efficacy of our correction for extremely small p-values. This is particularly important in the context of multiple comparisons, where practitioners care about the validity of exceedingly small p-values. To this end, the empirical cumulative distribution of the adjusted p-values is estimated under a standard Gaussian design with $n=4000$, $p=1200$ and $4.8\times 10^{5}$ p-values. The range $[0.1/p,12/p]$ is divided into points which are equi-spaced with a distance of $1/p$ between any two consecutive points. The estimated empirical CDF at each of these points is represented in blue in Figure 5. The estimated CDF is in near-perfect agreement with the diagonal, suggesting that the adjusted p-values computed using the rescaled chi-square distribution are remarkably close to a uniform, even when we zoom in at very small resolutions as would be the case when applying Bonferroni-style corrections.

Moderate sample sizes The final experiment studies the accuracy of our asymptotic result for moderately large samples. This is especially relevant for applications where the sample sizes are not too large. We repeat our numerical experiments with $n=200$, $p=60$ for i.i.d. Gaussian design, and $4.8\times 10^{5}$ p-values. The empirical CDF for these p-values are estimated and Figure 6 shows that the adjusted p-values are nearly uniformly distributed even for moderate sample sizes such as $n=200$.

4 Preliminaries

This section gathers a few preliminary results that will be useful throughout the paper. We start by collecting some facts regarding i.i.d. Gaussian random matrices.

Lemma 2**.**

Let $\bm{X}=[\bm{X}_{1},\bm{X}_{2},\ldots\bm{X}_{n}]^{\top}$ be an $n\times p$ matrix with i.i.d. standard Gaussian entries. Then

[TABLE]

[TABLE]

• Proof:

This is a straighforward application of [51, Corollary 5.35] and the union bound.

∎

Lemma 3**.**

Suppose $\bm{X}$ is an $n\times p$ matrix with entries i.i.d ${\mathcal{N}}(0,1)$, then there exists a constant $\epsilon_{0}$ such that whenever $0\leq\epsilon\leq\epsilon_{0}$ and $0\leq t\leq\sqrt{1-\epsilon}-\sqrt{p/n}$,

[TABLE]

with probability exceeding $1-2\exp\left(-\left(\frac{\left(1-\epsilon\right)t^{2}}{2}-H\left(\epsilon\right)\right)n\right)$. Here, $H(\epsilon)=-\epsilon\log\epsilon-(1-\epsilon)\log(1-\epsilon)$.

• Proof:

See Appendix A.1.

∎

The above facts are useful in establishing an eigenvalue lower bound on the Hessian of the log-likelihood function. Specifically, recall that

[TABLE]

and the result is this:

Lemma 4** (Likelihood Curvature Condition).**

Suppose that $p/n<1$ and that $\rho^{\prime\prime}(\cdot)\geq 0$. Then there exists a constant $\epsilon_{0}$ such that whenever $0\leq\epsilon\leq\epsilon_{0}$, with probability at least $1-2\exp\left(-n{H}\left(\epsilon\right)\right)-2\exp\left(-{n}/{2}\right)$, the matrix inequality

[TABLE]

holds simultaneously for all $\bm{\beta}\in\mathbb{R}^{p}$.

• Proof:

See Appendix A.2.

∎

The message of Lemma 4 is this: take $\epsilon>0$ to be a sufficiently small constant. Then

[TABLE]

for some non-increasing and positive function $\omega(\cdot)$ independent of $n$. This is a generalization of the strong convexity condition.

5 When is the MLE bounded?

5.1 Phase transition

In Section 2.2, we argued that the MLE is at infinity if we have less than two observations per dimension or $\kappa>1/2$. In fact, a stronger version of the phase transition phenemonon occurs in the sense that

[TABLE]

as soon as $\kappa<1/2$. This is formalized in the following theorem.

Theorem 4** (Norm Bound Condition).**

Fix any small constant $\epsilon>0$, and let $\hat{\bm{\beta}}$ be the MLE for a model with effective link satisfying the conditions from Section 2.3.3.

• (i)

If $p/n\geq 1/2+\epsilon$, then

[TABLE]

with probability exceeding $1-4\exp\left(-\epsilon^{2}n/8\right)$.

• (ii)

There exist universal constants $c_{1},c_{2},C_{2}>0$ such that if $p/n<1/2-c_{1}\epsilon^{3/4}$, then555When $\bm{X}_{i}\sim\mathcal{N}(\bm{0},\bm{\Sigma})$ for a general $\bm{\Sigma}\succ\bm{0}$, one has $\|\bm{\Sigma}^{1/2}\hat{\bm{\beta}}\|\lesssim{1}/{\epsilon^{2}}$ with high probability.**

[TABLE]

with probability at least $1-C_{2}\exp({-c_{2}\epsilon^{2}n})$.

These conclusions clearly continue to hold if $\hat{\bm{\beta}}$ is replaced by $\tilde{\bm{\beta}}$ (the MLE under the restricted model obtained on dropping the first predictor).

The rest of this section is devoted to proving this theorem. As we will see later, the fact that $\|\hat{\bm{\beta}}\|=O(1)$ is crucial for utilizing the AMP machinery in the absence of strong convexity.

5.2 Proof of Theorem 4

As in Section 2.5.1, we assume $\tilde{y}_{i}\equiv 1$ throughout this section, and hence the MLE reduces to

[TABLE]

5.2.1 Proof of Part (i)

Invoking [3, Theorem I] yields that if

[TABLE]

or equivalently, if $p/n\geq 1/2+\epsilon,$ then

[TABLE]

As is seen in Section 2.2, $\|\hat{\bm{\beta}}\|=\infty$ when $\left\{\bm{X}\bm{\beta}\mid\bm{\beta}\in\mathbb{R}^{p}\right\}\cap\mathbb{R}_{+}^{n}\neq\left\{\bm{0}\right\}$, establishing Part (i) of Theorem 4.

5.2.2 Proof of Part (ii)

We now turn to the regime in which $p/n\leq 1/2-O(\epsilon^{3/4})$, where $0<\epsilon<1$ is any fixed constant. Begin by observing that the least singular value of $\bm{X}$ obeys

[TABLE]

with probability at least 1-2\exp\big{(}-\frac{1}{2}\big{(}\frac{3}{4}-\frac{1}{\sqrt{2}}\big{)}^{2}n\big{)} (this follows from Lemma 3 using $\epsilon=0$). Then for any $\bm{\beta}\in\mathbb{R}^{p}$ obeying

[TABLE]

[TABLE]

we must have

[TABLE]

(a) follows since $t\leq\rho(t)$ and (b) is a consequence of (49). Continuing, (48) and (50) give

[TABLE]

This implies the following proposition: if the solution $\hat{\bm{\beta}}$—which necessarily satisfies $\ell_{0}(\hat{\bm{\beta}})\leq\ell_{0}(\bm{0})$—has norm exceeding $\|\hat{\bm{\beta}}\|\geq\frac{4\log 2}{\epsilon^{2}}$, then $\bm{X}\hat{\bm{\beta}}$ must fall within the cone

[TABLE]

Therefore, if one wishes to rule out the possibility of having $\|\hat{\bm{\beta}}\|\geq\frac{4\log 2}{\epsilon^{2}}$, it suffices to show that with high probability,

[TABLE]

This is the content of the remaining proof.

We would like to utilize tools from conic geometry [3] to analyze the probability of the event (52). Note, however, that $\mathcal{A}$ is not convex, while the theory developed in [3] applies only to convex cones. To bypass the non-convexity issue, we proceed in the following three steps:

Generate a set of $N=\exp\left(2\epsilon^{2}p\right)$ closed convex cones $\left\{\mathcal{B}_{i}\mid 1\leq i\leq N\right\}$ such that it forms a cover of $\mathcal{A}$ with probability exceeding $1-\exp\left(-\Omega(\epsilon^{2}p)\right)$. 2. 2.

Show that if $p<\left(\frac{1}{2}-2\sqrt{2}\epsilon^{\frac{3}{4}}-2H(2\sqrt{\epsilon})\right)n$ and if $n$ is sufficiently large, then

[TABLE]

for each $1\leq i\leq N$. 3. 3.

Invoke the union bound to reach

[TABLE]

where we have used the fact that

[TABLE]

Here, the last inequality holds if $\left(\frac{1}{2}-2\sqrt{2}\epsilon^{\frac{3}{4}}-10H(2\sqrt{\epsilon})-\frac{p}{n}\right)^{2}>24\epsilon^{2}$, or equivalently, $\frac{p}{n}<\frac{1}{2}-2\sqrt{2}\epsilon^{\frac{3}{4}}-10H(2\sqrt{\epsilon})-\sqrt{24}\epsilon$.

Taken collectively, these steps establish the following claim: if $\frac{p}{n}<\frac{1}{2}-2\sqrt{2}\epsilon^{\frac{3}{4}}-10H(2\sqrt{\epsilon})-\sqrt{24}\epsilon$, then

[TABLE]

thus establishing Part (ii) of Theorem 4. We defer the complete details of the preceding steps to Appendix D.

6 Asymptotic $\ell_{2}$ error of the

MLE

This section aims to establish Theorem 2, which characterizes precisely the asymptotic squared error of the MLE $\hat{\bm{\beta}}$ under the global null $\bm{\beta}=\bm{0}$. As described in Section 2.5.1, it suffices to assume that $\hat{\bm{\beta}}$ is the solution to the following problem

[TABLE]

In what follows, we derive the asymptotic convergence of $\|\hat{\bm{\beta}}\|$ under the assumptions from our main theorem.

Theorem 5**.**

Under the assumptions of Theorem 1, the solution $\hat{\bm{\beta}}$ to (53) obeys

[TABLE]

Theorem 5 is derived by invoking the AMP machinery [6, 7, 32]. The high-level idea is the following: in order to study $\hat{\bm{\beta}}$, one introduces an iterative algorithm (called AMP) where a sequence of iterates $\hat{\bm{\beta}}^{t}$ is formed at each time $t$. The algorithm is constructed so that the iterates asymptotically converge to the MLE in the sense that

[TABLE]

On the other hand, the asymptotic behavior (asymptotic in $n$) of $\hat{\bm{\beta}}^{t}$ for each $t$ can be described accurately by a scalar sequence $\{\tau_{t}\}$—called state evolution (SE)—following certain update equations [6]. This, in turn, provides a characterization of the $\ell_{2}$ loss of $\hat{\bm{\beta}}$.

Further, in order to prove Theorem 2, one still needs to justify

• (a)

the existence of a solution to the system of equations (23) and (24),

• (b)

and the existence of a fixed point for the iterative map governing the SE sequence updates.

We will elaborate on these steps in the rest of this section.

6.1 State evolution

We begin with the SE sequence $\{\tau_{t}\}$ introduced in [18]. Starting from some initial point $\tau_{0}$, we produce two sequences $\{b_{t}\}$ and $\{\tau_{t}\}$ following a two-step procedure.

• •

For $t=0,1,\ldots$:

• –

Set $b_{t}$ to be the solution in $b$ to

[TABLE]

• –

Set $\tau_{t+1}$ to be

[TABLE]

Suppose that for given any $\tau>0$, the solution in $b$ to (56) with $\tau_{t}=\tau$ exists and is unique, then one can denote the solution as $b(\tau)$, which in turn allows one to write the sequence $\{\tau_{t}\}$ as

[TABLE]

with the variance map

[TABLE]

As a result, if there exists a fixed point $\tau_{\ast}$ obeying $\mathcal{V}(\tau_{\ast}^{2})=\tau_{\ast}^{2}$ and if we start with $\tau_{0}=\tau_{\ast}$, then by induction,

[TABLE]

Notably, $(\tau_{\ast},b_{\ast})$ solves the system of equations (23) and (24). We shall work with this choice of initial condition throughout our proof.

The preceding arguments hold under two conditions: (i) the solution to (24) exists and is unique for any $\tau_{t}>0$; (ii) the variance map (58) admits a fixed point. To verify these two conditions, we make two observations.

• •

Condition (i) holds if one can show that the function

[TABLE]

is strictly monotone for any given $\tau>0$, and that $\lim_{b\rightarrow 0}G(b)<\kappa<\lim_{b\rightarrow\infty}G(b)$.

• •

Since $\mathcal{V}(\cdot)$ is a continuous function, Condition (ii) becomes self-evident once we show that $\mathcal{V}(0)>0$ and that there exists $\tau>0$ obeying $\mathcal{V}(\tau^{2})<\tau^{2}$. The behavior of the variance map is illustrated in Figure 7 for the logistic and probit regression when $\kappa=0.3$. One can in fact observe that the fixed point is unique. For other values of $\kappa$, the variance map shows the same behavior.

In fact, the aforementioned properties can be proved for a certain class of effective links, as summarized in the following lemmas. In particular, they can be shown for the logistic and the probit models.

Lemma 5**.**

Suppose the effective link $\rho$ satisfies the following two properties:

1. (a)

$\rho^{\prime}$* is log-concave.* 2. (b)

For any fixed $\tau>0$ and any fixed $z$, $b\rho^{\prime\prime}(\mathsf{prox}_{b\rho}(\tau z))\rightarrow\infty$ when $b\rightarrow\infty$.

Then for any $\tau>0$, the function $G(b)$ defined in (59) is an increasing function in $b$ ($b>0$), and the equation

[TABLE]

has a unique positive solution.

• Proof:

See Appendix B.

∎

Lemma 6**.**

Suppose that $0<\kappa<1/2$ and that $\rho=\log(1+e^{t})$ or $\rho=-\log\Phi(-t)$. Then

• (i)

$\mathcal{V}(0)>0$;

• (ii)

$\mathcal{V}(\tau^{2})<\tau^{2}$* for some sufficiently large $\tau^{2}$.*

• Proof:

See Appendix C and the supplemental material [47].

∎

Remark 1**.**

A byproduct of the proof is that the following relations hold for any constant $0<\kappa<1/2$:

• •

In the logistic case,

[TABLE]

• •

*In the probit case, *

[TABLE]

Remark 2**.**

Lemma 6 is proved for the two special effective link functions, the logistic and the probit cases. However, the proof sheds light on general conditions on the effective link that suffice for the lemma to hold. Such general sufficient conditions are also discussed in the supplemental material [47].

6.2 AMP recursion

In this section, we construct the AMP trajectory tracked by two sequences $\{\hat{\bm{\beta}}^{t}(n)\in\mathbb{R}^{p}\}$ and $\{\bm{\eta}^{t}(n)\in\mathbb{R}^{n}\}$ for $t\geq 0$. Going forward we suppress the dependence on $n$ to simplify presentation. Picking $\hat{\bm{\beta}}^{0}$ such that

[TABLE]

and taking $\bm{\eta}^{-1}=\bm{0}$ and $b_{-1}=0$, the AMP path is obtained via Algorithm 1, which is adapted from the algorithm in [18, Section 2.2].

As asserted by [18], the SE sequence $\{\tau_{t}\}$ introduced in Section 6.1 proves useful as it offers a formal procedure for predicting operating characteristics of the AMP iterates at any fixed iteration. In particular it assigns predictions to two types of observables: observables which are functions of the $\hat{\bm{\beta}}^{t}$ sequence and those which are functions of $\bm{\eta}^{t}$. Repeating identical argument as in [18, Theorem 3.4], we obtain

[TABLE]

6.3 AMP converges to the MLE

We are now in position to show that the AMP iterates $\{\hat{\bm{\beta}}^{t}\}$ converge to the MLE in the large $n$ and $t$ limit. Before continuing, we state below two properties that are satisfied under our assumptions.

• •

The MLE $\hat{\bm{\beta}}$ obeys

[TABLE]

almost surely.

• •

And there exists some non-increasing continuous function $0<\omega\left(\cdot\right)<1$ independent of $n$ such that

[TABLE]

In fact, the norm bound (65) follows from Theorem 4 together with Borel-Cantelli, while the likelihood curvature condition (66) is an immediate consequence of Lemma 4. With this in place, we have:

Theorem 6**.**

Suppose (65) and (66) hold. Let $(\tau_{\ast},b_{\ast})$ be a solution to the system (23) and (24), and assume that $\lim_{n\rightarrow\infty}\|\hat{\bm{\beta}}^{0}\|^{2}=\tau_{\ast}^{2}$. Then the AMP trajectory as defined in Algorithm 1 obeys

[TABLE]

Taken collectively, Theorem 6 and Eqn. (64) imply that

[TABLE]

thus establishing Theorem 5. In addition, an upshot of these theorems is a uniqueness result:

Corollary 2**.**

The solution to the system of equations (23) and (24) is unique.

• Proof:

When the AMP trajectory $\hat{\bm{\beta}}^{t}$ is started with the initial condition from Theorem 6, $\lim_{n\rightarrow\infty}\|\hat{\bm{\beta}}\|^{2}=_{\text{a.s.}}\tau_{\ast}^{2}$. This holds for any $\tau_{\ast}$ such that $(\tau_{\ast},b_{\ast})$ is a solution to (23) and (24). However, since the MLE problem is strongly convex and hence admits a unique solution $\hat{\bm{\beta}}$, this implies that $\tau_{\ast}$ must be unique, which together with the monotonicity of $G(\cdot)$ (cf. (59)) implies that $b_{\ast}$ is unique as well.

∎

• Proof of Theorem 6:

To begin with, repeating the arguments in [18, Lemma 6.9] we reach

[TABLE]

To show that the AMP iterates converge to the MLE, we shall analyze the log-likelihood function. Recall from Taylor’s theorem that

[TABLE]

holds for some $0<\lambda<1$. To deal with the quadratic term, we would like to control the Hessian of the likelihood at a point between $\hat{\bm{\beta}}$ and $\hat{\bm{\beta}}^{t}$. Invoking the likelihood curvature condition (66), one has

[TABLE]

with high probability. Apply Cauchy-Schwarz to yield that with exponentially high probability,

[TABLE]

where the last inequality follows since $0<\omega(\cdot)<1$ and $\omega(\cdot)$ is non-decreasing.

It remains to control $\|\nabla\ell(\hat{\bm{\beta}}^{t})\|$. The identity $\Psi(z;b_{*})=z-\mathsf{prox}_{b_{*}\rho}(z)$ and (61) give

[TABLE]

In addition, substituting $\Psi\left(z;b\right)=b\rho^{\prime}(\mathsf{prox}_{\rho b}(z))$ into (63) yields

[TABLE]

We are now ready to bound $\|\nabla\ell(\hat{\bm{\beta}}^{t})\|$. Recalling that

[TABLE]

and that $\sup_{z}\rho^{\prime\prime}(z)<\infty$, we have

[TABLE]

This establishes that with probability at least $1-c_{1}e^{-c_{2}n}$,

[TABLE]

Using (43) together with Borel-Cantelli yields $\lim_{n\rightarrow\infty}\|\bm{X}\|/\sqrt{n}<\infty$ almost surely. Further, it follows from (64) that $\lim_{n\rightarrow\infty}\|\hat{\bm{\beta}}^{t}\|$ is finite almost surely as $\tau_{\ast}<\infty.$ These taken together with (65), (68) and (69) yield

[TABLE]

as claimed.

∎

7 Likelihood ratio analysis

This section presents the analytical details for Section 2.5.2, which relates the log-likelihood ratio statistic $\Lambda_{i}$ with $\hat{\beta}_{i}$. Recall from (36) that the LLR statistic for testing $\beta_{1}=0$ vs. $\beta_{1}\neq 0$ is given by

[TABLE]

where

[TABLE]

and $\gamma_{i}$ lies between $\bm{X}_{i}^{\top}\hat{\bm{\beta}}$ and $\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{\beta}}$. The asymptotic distribution of $\Lambda_{1}$ claimed in Theorem 3 immediately follows from the result below, whose proof is the subject of the rest of this section.

Theorem 7**.**

Let $(\tau_{\ast},b_{\ast})$ be the unique solution to the system of equations (23) and (24), and define

[TABLE]

*Suppose $p/n\rightarrow\kappa\in(0,1/2)$ . Then *

• (a)

the log-likelihood ratio statistic obeys

[TABLE]

• (b)

and the scalar $\tilde{\alpha}$ converges,

[TABLE]

7.1 More notations and preliminaries

Before proceeding, we introduce some notations that will be used throughout. For any matrix $\bm{X}$, denote by $X_{ij}$ and $\bm{X}_{\cdot j}$ its $(i,j)$-th entry and $j$th column, respectively. We denote an analogue $\bm{r}=\{r_{i}\}_{1\leq i\leq n}$ (resp. $\tilde{\bm{r}}=\{\tilde{r}_{i}\}_{1\leq i\leq n}$) of residuals in the full (resp. reduced) model by

[TABLE]

As in (75), set

[TABLE]

where $\gamma_{i}^{*}$ is between $\bm{X}_{i}^{\top}\hat{\bm{\beta}}$ and $\bm{X}_{i}^{\top}\tilde{\bm{b}}$, and $\tilde{\bm{b}}$ is to be defined later in Section 7.2. Further, as in (76), introduce the Gram matrices

[TABLE]

Let $\tilde{\bm{G}}_{(i)}$ denote the version of $\tilde{\bm{G}}$ without the term corresponding to the $i{\text{th}}$ observation, that is,

[TABLE]

Additionally, let $\hat{\bm{\beta}}_{[-i]}$ be the MLE when the $i{\text{th}}$ observation is dropped and let $\bm{G}_{[-i]}$ be the corresponding Gram matrix,

[TABLE]

Further, let $\tilde{\bm{\beta}}_{[-i]}$ be the MLE when the first predictor and $i$th observation are removed, i.e.

[TABLE]

Below $\tilde{\bm{G}}_{[-i]}$ is the corresponding version of $\tilde{\bm{G}}$,

[TABLE]

For these different versions of $\bm{G}$, their least eigenvalues are all bounded away from 0, as asserted by the following lemma.

Lemma 7**.**

There exist some absolute constants $\lambda_{\mathrm{lb}},C,c>0$ such that

[TABLE]

Moreover, the same result holds for $\tilde{\bm{G}}$, $\bm{G}_{\hat{\bm{\beta}},\tilde{\bm{b}}}$, $\tilde{\bm{G}}_{(i)}$, $\bm{G}_{[-i]}$ and $\tilde{\bm{G}}_{[-i]}$ for all $i\in[n]$.

• Proof:

This result follows directly from Lemma 2, Lemma 4, and Theorem 4.

∎

Throughout the rest of this section, we restrict ourselves (for any given $n$) to the following event:

[TABLE]

By Lemma 7, $\mathcal{A}_{n}$ arises with exponentially high probability, i.e.

[TABLE]

7.2 A surrogate for the MLE

In view of (74), the main step in controlling $\Lambda_{1}$ consists of characterizing the differences $\bm{X}\hat{\bm{\beta}}-\tilde{\bm{X}}\tilde{\bm{\beta}}$ or \hat{\bm{\beta}}-\left[\begin{array}[]{c}0\\ \tilde{\bm{\beta}}\\ \end{array}\right]. Since the definition of $\hat{\bm{\beta}}$ is implicit and not amenable to direct analysis, we approximate $\hat{\bm{\beta}}$ by a more amenable surrogate $\tilde{\bm{b}}$, an idea introduced in [21, 20, 19]. We collect some properties of the surrogate which will prove valuable in the subsequent analysis.

To begin with, our surrogate is

[TABLE]

where $\tilde{\bm{G}}$ is defined in (81),

[TABLE]

and $\tilde{b}_{1}$ is a scalar to be specified later. This vector is constructed in the hope that

[TABLE]

where $\hat{\bm{\beta}}_{2:p}$ contains the $2{\text{nd}}$ through $p{\text{th}}$ components of $\hat{\bm{\beta}}$.

Before specifying $\tilde{b}_{1}$, we shall first shed some insights into the remaining terms in $\tilde{\bm{b}}$. By definition,

[TABLE]

Employing the first-order approximation of $\nabla\ell(\cdot)$ gives

[TABLE]

Suppose $\hat{\bm{\beta}}_{2:p}$ is well approximated by $\tilde{\bm{\beta}}$. Then all but the first coordinates of $\nabla\ell(\tilde{\bm{\beta}})$ and $\nabla\ell\left(\begin{bmatrix}0\\ \hat{\bm{\beta}}\end{bmatrix}\right)$ are also very close to each other. Therefore, taking the $2{\text{nd}}$ through $p{\text{th}}$ components of (90) and approximating them by zero give

[TABLE]

This together with a little algebra yields

[TABLE]

which coincides with (89). In fact, for all but the $1{\text{st}}$ entries, $\tilde{\bm{b}}$ is constructed by moving $\tilde{\bm{\beta}}$ one-step in the direction which takes it closest to $\hat{\bm{\beta}}$.

Next, we come to discussing the scalar $\tilde{b}_{1}$. Introduce the projection matrix

[TABLE]

and define $\tilde{b}_{1}$ as

[TABLE]

where $\tilde{\bm{r}}$ comes from (79). In fact, the expression $\tilde{b}_{1}$ is obtained through similar (but slightly more complicated) first-order approximation as for $\tilde{\bm{b}}_{2:p}$, in order to make sure that $b_{1}\approx\hat{\beta}_{1}$; see [21, Pages 14560-14561] for a detailed description.

We now formally justify that the surrogate $\tilde{\bm{b}}$ and the MLE $\hat{\bm{\beta}}$ are close to each other.

Theorem 8**.**

The MLE $\hat{\bm{\beta}}$ and the surrote $\tilde{\bm{b}}$ (87) obey

[TABLE]

[TABLE]

and

[TABLE]

with probability tending to one as $n\rightarrow\infty$.

• Proof:

See Section 7.4. ∎

The global accuracy (93) immediately leads to a coordinate-wise approximation result between $\hat{\beta}_{1}$ and $\tilde{b}_{1}$.

Corollary 3**.**

With probability tending to one as $n\rightarrow\infty$,

[TABLE]

Another consequence from Theorem 8 is that the value $\bm{X}_{i}^{\top}\hat{\bm{\beta}}$ in the full model and its counterpart $\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{\beta}}$ in the reduced model are uniformly close.

Corollary 4**.**

The values $\bm{X}_{i}^{\top}\hat{\bm{\beta}}$ and $\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{\beta}}$ are uniformly close in the sense that

[TABLE]

holds with probability approaching one as $n\rightarrow\infty$.

• Proof:

Note that

[TABLE]

The second term in the right-hans side is upper bounded by $n^{-1/2+o(1)}$ with probability $1-o(1)$ according to Theorem 8. Invoking Lemma 2 and Theorem 8 and applying Cauchy-Schwarz inequality yield that the first term is $O(n^{-1/2+o(1)})$ with probability $1-o(1)$. This establishes the claim.

∎

7.3 Analysis of the likelihood-ratio statistic

We are now positioned to use our surrogate $\tilde{\bm{b}}$ to analyze the likelihood-ratio statistic. In this subsection we establish Theorem 7(a). The proof for Theorem 7(b) is deferred to Appendix I.

Recall from (36) that

[TABLE]

To begin with, Corollary 4 together with the assumption $\sup_{z}\rho^{\prime\prime\prime}(z)<\infty$ implies that

[TABLE]

with probability $1-o(1)$. Hence, $I_{3}$ converges to zero in probability.

Reorganize the quadratic term as follows:

[TABLE]

We control each of the three terms in the right-hand side of (98).

• •

Since $\sup_{z}\rho^{\prime\prime}(z)<\infty$, the first term in (98) is bounded by

[TABLE]

with probability $1-o(1)$, by an application of Theorem 8 and Lemma 2.

• •

From the definition of $\tilde{\bm{b}}$, the second term can be upper bounded by

[TABLE]

with probability $1-o(1)$, where the last line follows from a combination of Theorem 8, Lemma 2 and the following lemma.

Lemma 8**.**

Let $\tilde{\bm{G}}$ and $\bm{w}$ be as defined in (81) and (88), respectively. Then

[TABLE]

• Proof:

See Appendix E. ∎

• •

The third term in (98) can be decomposed as

[TABLE]

for some $\tilde{\gamma}_{i}$ between $\bm{X}_{i}^{\top}\hat{\bm{\beta}}$ and $\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{\beta}}$. From Theorem 8 and Corollary 4, the first term in (100) is $O(n^{-1/2+o(1)})$ with probability $1-o(1)$. Hence, the only remaining term is the second.

In summary, we have

[TABLE]

where $\bm{v}:=\tilde{b}_{1}\begin{bmatrix}1\\ -\tilde{\bm{G}}^{-1}\bm{w}\end{bmatrix}$ according to (87). On simplification, the quadratic form reduces to

[TABLE]

recalling the definitions (81), (88), and (91). Hence, the log-likelihood ratio $2\Lambda_{1}$ simplifies to $n\tilde{b}_{1}^{2}\xi+o_{P}(1)$ on $\mathcal{A}_{n}$.

Finally, rewrite $\bm{v}^{\top}\bm{X}^{\top}\bm{D}_{\tilde{\bm{\beta}}}\bm{X}\bm{v}$ as $n(\tilde{b}_{1}^{2}-\hat{\beta}_{1}^{2})\xi+n\hat{\beta}_{1}^{2}\xi$. To analyze the first term, note that

[TABLE]

with probability $1-o(1)$ in view of Theorem 8 and Corollary 3. It remains to analyze $\xi$. Recognize that $\bm{X}_{\cdot 1}$ is independent of $\bm{D}_{\tilde{\bm{\beta}}}^{1/2}\bm{H}\bm{D}_{\tilde{\bm{\beta}}}^{1/2}$. Applying the Hanson-Wright inequality [27, 44] and the Sherman-Morrison-Woodbury formula (e.g. [26]) leads to the following lemma:

Lemma 9**.**

Let $\tilde{\alpha}=\frac{1}{n}\mathrm{Tr}(\tilde{\bm{G}}^{-1})$, where $\tilde{\bm{G}}=\frac{1}{n}\tilde{\bm{X}}^{\top}\bm{D}_{\tilde{\bm{\beta}}}\tilde{\bm{X}}$. Then one has

[TABLE]

with probability approaching one as $n\rightarrow\infty$.

• Proof:

See Appendix F. ∎

In addition, if one can show that $\tilde{\alpha}$ is bounded away from zero with probability $1-o(1)$, then it is seen from Lemma 9 that

[TABLE]

To justify the above claim, we observe that since $\rho^{\prime\prime}$ is bounded, $\lambda_{\max}(\tilde{\bm{G}})\lesssim\lambda_{\max}(\tilde{\bm{X}}^{\top}\tilde{\bm{X}})/n\lesssim 1$ with exponentially high probability (Lemma 2). This yields

[TABLE]

with probability $1-o(1)$. On the other hand, on $\mathcal{A}_{n}$ one has

[TABLE]

Hence, it follows that $\xi=\Omega(1)$ with probability $1-o(1)$. Putting this together with (102) gives the approximation

[TABLE]

Taken collectively (101), (104) and (105) yields the desired result

[TABLE]

7.4 Proof of Theorem 8

This subsection outlines the main steps for the proof of Theorem 8. To begin with, we shall express the difference $\hat{\bm{\beta}}-\tilde{\bm{b}}$ in terms of the gradient of the negative log-likelihood function. Note that $\nabla\ell(\hat{\bm{\beta}})=\bm{0}$, and hence

[TABLE]

where $\gamma_{i}^{*}$ is between $\bm{X}_{i}^{\top}\hat{\bm{\beta}}$ and $\bm{X}_{i}^{\top}\tilde{\bm{b}}$. Recalling the notation introduced in (81), this can be rearranged as

[TABLE]

Hence, on $\mathcal{A}_{n}$, this yields

[TABLE]

The next step involves expressing $\nabla\ell(\tilde{\bm{b}})$ in terms of the difference $\tilde{\bm{b}}-\begin{bmatrix}0\\ \tilde{\bm{\beta}}\end{bmatrix}$.

Lemma 10**.**

On the event $\mathcal{A}_{n}$ (85), the negative log-likelihood evaluated at the surrogate $\tilde{\bm{b}}$ obeys

[TABLE]

where $\gamma_{i}^{*}$ is some quantity between $\bm{X}_{i}^{\top}\tilde{\bm{b}}$ and $\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{\beta}}$.

• Proof:

The proof follows exactly the same argument as in the proof of [20, Proposition 3.11], and is thus omitted.

∎

The point of expressing $\nabla\ell(\tilde{\bm{b}})$ in this way is that the difference \tilde{\bm{b}}-\left[\begin{array}[]{c}0\\ \tilde{\bm{\beta}}\end{array}\right] is known explicitly from the definition of $\tilde{\bm{b}}$. Invoking Lemma 10 and the definition (87) allows one to further upper bound (106) as

[TABLE]

with probability at least $1-\exp(-\Omega(n))$. The last inequality here comes from our assumption that $\sup_{z}|\rho^{\prime\prime\prime}(z)|<\infty$ together with Lemmas 2 and 8.

In order to bound (110), we first make use of the definition of $\tilde{\bm{b}}$ to reach

[TABLE]

The following lemma provides an upper bound on $\sup_{i}|X_{i1}-\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{G}}^{-1}\bm{w}|$.

Lemma 11**.**

With $\tilde{\bm{G}}$ and $\bm{w}$ as defined in (81) and (88),

[TABLE]

• Proof:

See Appendix G.

∎

In view of Lemma 11, the second term in the right-hand side of (111) is bounded above by $n^{o(1)}$ with high probability. Thus, in both the bounds (110) and (111), it only remains to analyze the term $\tilde{b}_{1}$. To this end, we control the numerator and the denominator of $\tilde{b}_{1}$ separately.

• •

Recall from the definition (92) that the numerator of $\tilde{b}_{1}$ is given by $\bm{X}_{\cdot 1}^{\top}\tilde{\bm{r}}$ and that $\tilde{\bm{r}}$ is independent of $\bm{X}_{\cdot 1}$. Thus, conditional on $\tilde{\bm{X}}$, the quantity $\bm{X}_{\cdot 1}^{\top}\tilde{\bm{r}}$ is distributed as a Gaussian with mean zero and variance

[TABLE]

Since $|\rho^{\prime}(x)|=O(|x|)$, the variance is bounded by

[TABLE]

with probability at least $1-\exp(-\Omega(n)))$, a consequence from Theorem 4 and Lemma 2. Therefore, with probability $1-o(1)$, we have

[TABLE]

• •

We now move on to the denominator of $\tilde{b}_{1}$ in (92). In the discussion following Lemma 9 we showed $\frac{1}{n}\bm{X}_{\cdot 1}^{\top}\bm{D}_{\tilde{\bm{\beta}}}^{1/2}\bm{H}\bm{D}_{\tilde{\bm{\beta}}}^{1/2}\bm{X}_{\cdot 1}=\Omega(1)$ with probability $1-o(1)$.

Putting the above bounds together, we conclude

[TABLE]

Substitution into (110) and (111) yields

[TABLE]

with probability $1-o(1)$ as claimed.

8 Discussion

In this paper, we derived the high-dimensional asymptotic distribution of the LLR under our modelling assumptions. In particular, we showed that the LLR is inflated vis a vis the classical Wilks’ approximation and that this inflation grows as the dimensionality $\kappa$ increases. This inflation is typical of high-dimensional problems, and one immediate practical consequence is that it explains why classically computed p-values are completely off since they tend to be far too small under the null hypothesis. In contrast, we have shown in our simulations that our new limiting distribution yields reasonably accurate p-values in finite samples. Having said this, our work raises a few important questions that we have not answered and we conclude this paper with a couple of them.

• •

We expect that our results continue to hold when the covariates are not normally distributed, see Section 3 for some numerical evidence in this direction. To be more precise, we expect the same limiting distribution to hold when the variables are simply sub-Gaussian. If this were true, then this would imply that our rescaled chi-square has a form of universal validity.

• •

The major limitation of our work is arguably the fact that our limiting distribution holds under the global null; that is, under the assumption that all the regression coefficients vanish. It is unclear to us how the distribution would change in the case where the coefficients are not all zero. In particular, would the limiting distribution depend upon the unknown values of these coefficients? Are there assumptions under which it would not? Suppose for instance that we model the regression coefficients as i.i.d. samples from the mixture model

[TABLE]

where $0<\epsilon<1$ is a mixture parameter, $\delta_{0}$ is a point mass at zero and $\Pi^{\star}$ is a distribution with vanishing mass at zero. Then what would we need to know about $\epsilon$ and $\Pi^{\star}$ to compute the asymptotic distribution of the LLR under the null?

Acknowledgements

E. C. was partially supported by the Office of Naval Research under grant N00014-16-1-2712, and by the Math + X Award from the Simons Foundation. Y. C. and P. S. are grateful to Andrea Montanari for his help in understanding AMP and [18]. Y. C. thanks Kaizheng Wang and Cong Ma for helpful discussion about [20], and P. S. thanks Subhabrata Sen for several helpful discussions regarding this project. E. C. would like to thank Iain Johnstone for a helpful discussion as well.

Appendix A Proofs for Eigenvalue Bounds

A.1 Proof of Lemma 3

Fix $\epsilon\geq 0$ sufficiently small. For any given $S\subseteq[n]$ obeying $|S|=(1-\epsilon)n$ and $0\leq t\leq\sqrt{1-\epsilon}-\sqrt{p/n}$ it follows from [51, Corollary 5.35] that

[TABLE]

holds with probability at most $2\exp\left(-\frac{t^{2}|S|}{2}\right)=2\exp\left(-\frac{\left(1-\epsilon\right)t^{2}n}{2}\right)$. Taking the union bound over all possible subsets $S$ of size $(1-\epsilon)n$ gives

[TABLE]

where the last line is a consequence of the inequality ${n\choose(1-\epsilon)n}\leq e^{nH(\epsilon)}$ [16, Example 11.1.3].

A.2 Proof of Lemma 4

Define

[TABLE]

for any $B>0$ and any $\bm{\beta}$. Then

[TABLE]

If one also has $|S_{B}\left(\bm{\beta}\right)|\geq(1-\epsilon)n$ (for $\epsilon\geq 0$ sufficiently small), then this together with Lemma 3 implies that

[TABLE]

with probability at least $1-2\exp\left(-\left(\frac{\left(1-\epsilon\right)t^{2}}{2}-H\left(\epsilon\right)\right)n\right)$.

Thus if we can ensure that with high probability, $|S_{B}\left(\bm{\beta}\right)|\geq(1-\epsilon)n$ holds simultaneously for all $\bm{\beta}$, then we are done. From Lemma 2 we see that $\frac{1}{n}\left\|\bm{X}^{\top}\bm{X}\right\|\leq 9$ with probability exceeding $1-2\exp\left(-n/2\right)$. On this event,

[TABLE]

On the other hand, the definition of $S_{B}(\bm{\beta})$ gives

[TABLE]

Taken together, (116) and (117) yield

[TABLE]

with probability at least $1-2\exp(-n/2)$. Therefore, with probability $1-2\exp(-n/2)$, $\left|S_{3/\sqrt{\epsilon}}(\bm{\beta})\right|\geq\left(1-\epsilon\right)n$ holds simultaneously for all $\bm{\beta}$. Putting the above results together and setting $t=2\sqrt{\frac{H(\epsilon)}{1-\epsilon}}$ give

[TABLE]

simultaneously for all $\bm{\beta}$ with probability at least $1-2\exp\left(-nH\left(\epsilon\right)\right)-2\exp\left(-{n}/2\right)$.

Appendix B Proof of Lemma 5

Applying an integration by parts leads to

[TABLE]

with $\phi(z)=\frac{1}{\sqrt{2\pi}}\exp(-z^{2}/2)$. This reveals that

[TABLE]

where the second identity comes from [18, Proposition 6.4], and the last identity holds since $\phi^{\prime}(z)=-\phi^{\prime}(-z)$.

Next, we claim that

• (a)

The function $h\left(z\right):=\frac{\rho^{\prime}\left(z\right)}{1+b\rho^{\prime\prime}\left(z\right)}$ is increasing in $z$;

• (b)

$\mathrm{prox}_{b\rho}(z)$ is increasing in $z$.

These two claims imply that

[TABLE]

which combined with the fact $\phi^{\prime}(z)<0$ for $z>0$ reveals

[TABLE]

In other words, the integrand in (118) is positive, which allows one to conclude that $G^{\prime}(b)>0$.

We then move on to justify (a) and (b). For the first, the derivative of $h$ is given by

[TABLE]

Since $\rho^{\prime}$ is log concave, this directly yields $(\rho^{\prime\prime})^{2}-\rho^{\prime}\rho^{\prime\prime\prime}>0$. As $\rho^{\prime\prime}>0$ and $b\geq 0$, the above implies $h^{\prime}(z)>0$ for all $z$. The second claim follows from $\frac{\partial\mathrm{prox}_{b\rho}(z)}{\partial z}\geq\frac{1}{1+b\|\rho^{\prime\prime}\|_{\infty}}>0$ (cf. [18, Equation (56)]).

It remains to analyze the behavior of $G$ in the limits when $b\rightarrow 0$ and $b\rightarrow\infty$. From [18, Proposition 6.4], $G(b)$ can also be expressed as

[TABLE]

Since $\rho^{\prime\prime}$ is bounded and the integrand is at most $1$, the dominated convergence theorem gives

[TABLE]

When $b\rightarrow\infty$, $b\rho^{\prime\prime}(\mathsf{prox}_{b\rho}(\tau z))\rightarrow\infty$ for a fixed $z$. Again by applying the dominated convergence theorem,

[TABLE]

It follows that $\lim_{b\rightarrow 0}G(b)<\kappa<\lim_{b\rightarrow\infty}G(b)$ and, therefore, $G(b)=\kappa$ has a unique positive solution.

Remark 3**.**

Finally, we show that the logistic and the probit effective links obey the assumptions of Lemma 5. We work with a fixed $\tau>0$.

• •

A direct computation shows that $\rho^{\prime}$ is log-concave for the logistic model. For the probit, it is well-known that the reciprocal of the hazard function (also known as Mills’ ratio) is strictly log-convex **[4]**.

• •

To check the other condition, recall that the proximal mapping operator satisfies

[TABLE]

For a fixed $z$, we claim that if $b\rightarrow\infty$, $\mathsf{prox}_{b\rho}(\tau z)\rightarrow-\infty$. To prove this claim, we start by assuming that this is not true. Then either $\mathsf{prox}_{b\rho}(\tau z)$ is bounded or diverges to $\infty$. If it is bounded, the LHS above diverges to $\infty$ while the RHS is fixed, which is a contradiction. Similarly if $\mathsf{prox}_{b\rho}(\tau z)$ diverges to $\infty,$ the left-hand side of (119) diverges to $\infty$ while the right-hand side is fixed, which cannot be true as well. Further, when $b\rightarrow\infty$, we must have $\mathsf{prox}_{b\rho}(\tau z)\rightarrow-\infty$, $b\rho^{\prime}(\mathsf{prox}_{b\rho}(\tau z))\rightarrow\infty$, such that the difference of these two is $\tau z$. Observe that for the logistic, $\rho^{\prime\prime}(x)=\rho^{\prime}(x)(1-\rho^{\prime}(x))$ and for the probit, $\rho^{\prime\prime}(x)=\rho^{\prime}(x)(\rho^{\prime}(x)-x)$ **[45]**. Hence, combining the asymptotic behavior of $\mathsf{prox}_{b\rho}(\tau z)$ and $b\rho^{\prime}(\mathsf{prox}_{b\rho}(\tau z))$, we obtain that $b\rho^{\prime\prime}(\mathsf{prox}_{b\rho}(\tau z))$ diverges to $\infty$ in both models when $b\rightarrow\infty$.

Appendix C Proof of Lemma 6

C.1 Proof of Part (i)

Recall from [18, Proposition 6.4] that

[TABLE]

If we denote $c:=\mathsf{prox}_{b\rho}(0)$, then $b(0)$ is given by the following relation:

[TABLE]

as $\rho^{\prime\prime}(c)>0$ for any given $c>0$. In addition, since $\rho^{\prime}(c)>0$, we have

[TABLE]

where (a) comes from (20).

C.2 Proof of Part (ii)

We defer the proof of this part to the supplemental materials [47].

Appendix D Proof of Part (ii) of Theorem 4

As discussed in Section 5.2.2, it suffices to (1) construct a set $\left\{\mathcal{B}_{i}\mid 1\leq i\leq N\right\}$ that forms a cover of the cone $\mathcal{A}$ defined in (51), and (2) upper bound $\mathbb{P}\left\{\left\{\bm{X}\bm{\beta}\mid\bm{\beta}\in\mathbb{R}^{p}\right\}\cap\mathcal{B}_{i}\neq\left\{\bm{0}\right\}\right\}$. In what follows, we elaborate on these two steps.

• •

Step 1. Generate $N=\exp\left(2\epsilon^{2}p\right)$ i.i.d. points $\bm{z}^{(i)}\sim\mathcal{N}(\bm{0},\frac{1}{p}\bm{I}_{p})$, $1\leq i\leq N$, and construct a collection of convex cones

[TABLE]

In words, $\mathcal{C}_{i}$ consists of all directions that have nontrivial positive correlation with $\bm{z}^{(i)}$. With high probability, this collection $\left\{\mathcal{C}_{i}\mid 1\leq i\leq N\right\}$ forms a cover of $\mathbb{R}^{p}$, a fact which is an immediate consequence of the following lemma.

Lemma 12**.**

Consider any given constant $0<\epsilon<1$, and let $N=\exp\left(2\epsilon^{2}p\right)$. Then there exist some positive universal constants $c_{5},C_{5}>0$ such that with probability exceeding $1-C_{5}\exp\left(-c_{5}\epsilon^{2}p\right)$,

[TABLE]

holds simultaneously for all $\bm{x}\in\mathbb{R}^{p}$.

With our family $\left\{\mathcal{C}_{i}\mid 1\leq i\leq N\right\}$ we can introduce

[TABLE]

which in turn forms a cover of the nonconvex cone $\mathcal{A}$ defined in (51). To justify this, note that for any $\bm{u}\in\mathcal{A}$, one can find $i\in\{1,\cdots,N\}$ obeying $\bm{u}\in\mathcal{C}_{i}$, or equivalently, $\left\langle\bm{u},\frac{\bm{z}^{(i)}}{\|\bm{z}^{(i)}\|}\right\rangle\geq\epsilon\|\bm{u}\|$, with high probability. Combined with the membership to $\mathcal{A}$ this gives

[TABLE]

indicating that $\bm{u}$ is contained within some $\mathcal{B}_{i}$.

• •

Step 2. We now move on to control $\mathbb{P}\left\{\left\{\bm{X}\bm{\beta}\mid\bm{\beta}\in\mathbb{R}^{p}\right\}\cap\mathcal{B}_{i}\neq\left\{\bm{0}\right\}\right\}$. If the statistical dimensions of the two cones obey $\delta\left(\mathcal{B}_{i}\right)<n-\delta\left(\left\{\bm{X}\bm{\beta}\mid\bm{\beta}\in\mathbb{R}^{p}\right\}\right)=n-p$, then an application of [3, Theorem I] gives

[TABLE]

It then comes down to upper bounding $\delta(\mathcal{B}_{i})$, which is the content of the following lemma.

Lemma 13**.**

Fix $\epsilon>0$. When $n$ is sufficiently large, the statistical dimension of the convex cone $\mathcal{B}_{i}$ defined in (121) obeys

[TABLE]

where $H(x):=-x\log x-(1-x)\log(1-x)$.

Substitution into (122) gives

[TABLE]

Finally, we prove Lemmas 12-13 in the next subsections. These are the only remaining parts for the proof of Theorem 4.

D.1 Proof of Lemma 12

To begin with, it is seen that all $\|\bm{z}^{(i)}\|$ concentrates around $1$. Specifically, apply [29, Proposition 1] to get

[TABLE]

and set $t=3\epsilon^{2}p$ to reach

[TABLE]

Taking the union bound we obtain

[TABLE]

Next, we note that it suffices to prove Lemma 12 for all unit vectors $\bm{x}$. The following lemma provides a bound on $\left\langle\bm{z}^{(i)},\bm{x}\right\rangle$ for any fixed unit vector $\bm{x}\in\mathbb{R}^{p}$.

Lemma 14**.**

Consider any fixed unit vector $\bm{x}\in\mathbb{R}^{p}$ and any given constant $0<\epsilon<1$, and set $N=\exp\left(2\epsilon^{2}p\right)$. There exist positive universal constants $c_{5},c_{6},C_{6}>0$ such that

[TABLE]

Recognizing that Lemma 12 is a uniform result, we need to extend Lemma 14 to all $\bm{x}$ simultaneously, which we achieve via the standard covering argument. Specifically, one can find a set $\mathcal{C}:=\left\{\bm{x}^{(j)}\in\mathbb{R}^{p}\mid 1\leq j\leq K\right\}$ of unit vectors with cardinality $K=\left(1+2p^{2}\right)^{p}$ to form a cover of the unit ball of resolution $p^{-2}$ [51, Lemma 5.2]; that is, for any unit vector $\bm{x}\in\mathbb{R}^{p}$, there exists a $\bm{x}^{(j)}\in\mathcal{C}$ such that

[TABLE]

Apply Lemma 14 and take the union bound to arrive at

[TABLE]

with probability exceeding $1-K\exp\left\{-2\exp\left(\left(1-o(1)\right)\frac{7}{4}\epsilon^{2}p\right)\right\}\geq 1-\exp\left\{-2\left(1-o\left(1\right)\right)\exp\left(\left(1-o(1)\right)\frac{7}{4}\epsilon^{2}p\right)\right\}$. This guarantees that for each $\bm{x}^{(j)}$, one can find at least one $\bm{z}^{(i)}$ obeying

[TABLE]

This result together with (125) yields that with probability exceeding $1-C\exp\left(-c\epsilon^{2}p\right)$, for some universal constants $C,c>0$.

[TABLE]

holds simultaneously for all unit vectors $\bm{x}\in\mathbb{R}^{p}$. Since $\epsilon>0$ can be an arbitrary constant, this concludes the proof.

• Proof of Lemma 14:

Without loss of generality, it suffices to consider $\bm{x}=\bm{e}_{1}=[1,0,\cdots,0]^{\top}$. For any $t>0$ and any constant $\zeta>0$, it comes from [2, Theorem A.1.4] that

[TABLE]

Setting $t=1-\Phi\left(\zeta\sqrt{p}\right)$ gives

[TABLE]

Recall that for any $t>1$, one has $(t^{-1}-t^{-3})\phi(t)\leq 1-\Phi(t)\leq t^{-1}\phi(t)$ which implies that

[TABLE]

Taking $\zeta=\frac{1}{2}\epsilon$, we arrive at

[TABLE]

This justifies that

[TABLE]

as claimed.

∎

D.2 Proof of Lemma 13

First of all, recall from the definition (17) that

[TABLE]

where $\bm{g}\sim\mathcal{N}\left(\bm{0},\bm{I}_{n}\right)$, and $\mathcal{D}_{i}$ is a superset of $\mathcal{B}_{i}$ defined by

[TABLE]

Recall from the triangle inequality that

[TABLE]

Since $\bm{0}\in\mathcal{D}_{i}$, this implies that

[TABLE]

revealing that

[TABLE]

In what follows, it suffices to look at the set of $\bm{u}$’s within $\mathcal{D}_{i}$ obeying $\|\bm{u}\|\leq 2\|\bm{g}\|$, which verify

[TABLE]

It is seen that

[TABLE]

Regarding the first term of (130), we first recognize that

[TABLE]

where the last inequality follows from the constraint (129). As a consequence,

[TABLE] 2. 2.

Next, we turn to the second term of (130), which can be bounded by

[TABLE]

where the last inequality follows from the constraint (129).

Putting the above results together, we have

[TABLE]

for any $\bm{u}\in\mathcal{D}_{i}$ obeying $\|\bm{u}\|\leq 2\|\bm{g}\|$, whence

[TABLE]

Finally, it follows from [29, Proposition 1] that for any $t>2\sqrt{\epsilon}n$,

[TABLE]

which together with the union bound gives

[TABLE]

This gives

[TABLE]

for any given $\epsilon>0$ with the proviso that $n$ is sufficiently large. This combined with (131) yields

[TABLE]

as claimed.

Appendix E Proof of Lemma 8

Throughout, we shall restrict ourselves on the event $\mathcal{A}_{n}$ as defined in (85), on which $\tilde{\bm{G}}\succeq\lambda_{\mathrm{lb}}\bm{I}$. Recalling the definitions of $\tilde{\bm{G}}$ and $\bm{w}$ from (81) and (88), we see that

[TABLE]

If we let the singular value decomposition of $\frac{1}{\sqrt{n}}\bm{D}_{\tilde{\bm{\beta}}}^{1/2}\tilde{\bm{X}}$ be $\bm{U}\bm{\Sigma}\bm{V}^{\top}$, then a little algebra gives $\bm{\Sigma}\succeq\sqrt{\lambda_{\mathrm{lb}}}\bm{I}$ and

[TABLE]

Substituting this into (133) and using the fact $\|\bm{X}_{\cdot 1}\|^{2}\lesssim n$ with high probability (by Lemma 2), we obtain

[TABLE]

with probability at least $1-\exp(-\Omega(n))$.

Appendix F Proof of Lemma 9

Throughout this and the subsequent sections, we consider $H_{n}$ and $K_{n}$ to be two diverging sequences with the following properties:

[TABLE]

for any constants $c_{i}>0$, $i=1,2$ and any $\epsilon>0$. This lemma is an analogue of [20, Proposition 3.18]. We modify and adapt the proof ideas to establish the result in our setup. Throughout we shall restrict ourselves to the event $\mathcal{A}_{n}$, on which $\tilde{\bm{G}}\succeq\lambda_{\mathrm{lb}}\bm{I}$.

Due to independence between $\bm{X}_{\cdot 1}$ and $\{\bm{D}_{\tilde{\bm{\beta}}},\bm{H}\}$, one can invoke the Hanson-Wright inequality [44, Theorem 1.1] to yield

[TABLE]

where $\|.\|_{\mathrm{F}}$ denotes the Frobenius norm. Choose t=C^{2}\big{\|}\bm{D}_{\tilde{\bm{\beta}}}^{1/2}\bm{H}\bm{D}_{\tilde{\bm{\beta}}}^{1/2}\big{\|}H_{n}/\sqrt{n} with $C>0$ a sufficiently large constant, and take $H_{n}$ to be as in (134). Substitution into the above inequality and unconditioning give

[TABLE]

for some universal constants $C,c>0$.

We are left to analyzing \mathrm{Tr}\big{(}\bm{D}_{\tilde{\bm{\beta}}}^{1/2}\bm{H}\bm{D}_{\tilde{\bm{\beta}}}^{1/2}\big{)}. Recall from the definition (91) of $\bm{H}$ that

[TABLE]

and, hence,

[TABLE]

This requires us to analyze $\tilde{\bm{G}}^{-1}$ carefully. To this end, recall that the matrix $\tilde{\bm{G}}_{(i)}$ defined in (82) obeys

[TABLE]

Invoking Sherman-Morrison-Woodbury formula (e.g. [26]), we have

[TABLE]

It follows that

[TABLE]

which implies that

[TABLE]

The relations (136) and (138) taken collectively reveal that

[TABLE]

We shall show that the trace above is close to $\mathrm{Tr}(\bm{I}-\bm{H})$ up to some factors. For this purpose we analyze the latter quantity in two different ways. To begin with, observe that

[TABLE]

On the other hand, it directly follows from the definition of $\bm{H}$ and (138) that the $i{\text{th}}$ diagonal entry of $\bm{H}$ is given by

[TABLE]

Applying this relation, we can compute $\mathrm{Tr}(\bm{I}-\bm{H})$ analytically as follows:

[TABLE]

where $\tilde{\alpha}:=\frac{1}{n}\mathrm{Tr}\left(\tilde{\bm{G}}^{-1}\right)$.

Observe that the first quantity in the right-hand side above is simply \tilde{\alpha}\mathrm{Tr}\big{(}\bm{D}_{\tilde{\bm{\beta}}}^{1/2}\bm{H}\bm{D}_{\tilde{\bm{\beta}}}^{1/2}\big{)}. For simplicity, denote

[TABLE]

Note that $\tilde{\bm{G}}_{(i)}\succ\bm{0}$ on $\mathcal{A}_{n}$ and that $\rho^{\prime\prime}>0$. Hence the denominator in the second term in (142) is greater than $1$ for all $i$. Comparing (140) and (142), we deduce that

[TABLE]

on $\mathcal{A}_{n}$. It thus suffices to control $\sup_{i}|\eta_{i}|$. The above bounds together with Lemma (86) and the proposition below complete the proof.

Proposition 1**.**

Let $\eta_{i}$ be as defined in (143). Then there exist universal constants $C_{1},C_{2},C_{3}>0$ such that

[TABLE]

where $K_{n},H_{n}$ are diverging sequences as specified in (134)

• Proof of Proposition 1:

Fix any index $i$. Recall that $\tilde{\bm{\beta}}_{[-i]}$ is the MLE when the $1{\text{st}}$ predictor and $i{\text{th}}$ observation are removed. Also recall the definition of $\tilde{\bm{G}}_{[-i]}$ in (84). The proof essentially follows three steps. First, note that $\tilde{\bm{X}}_{i}$ and $\tilde{\bm{G}}_{[-i]}$ are independent. Hence, an application of the Hanson-Wright inequality [44] yields that

[TABLE]

We choose t=C^{2}\big{\|}\tilde{\bm{G}}_{[-i]}^{-1}\big{\|}H_{n}/\sqrt{n}, where $C>0$ is a sufficiently large constant. Now marginalizing gives

[TABLE]

where $C^{\prime}>0$ is a sufficiently large constant. On $\mathcal{A}_{n}$, the spectral norm \big{\|}\tilde{\bm{G}}_{(i)}^{-1}\big{\|} is bounded above by $\lambda_{\mathrm{lb}}$ for all $i$. Invoking (86) we obtain that there exist universal constants $C_{1},C_{2},C_{3}>0$ such that

[TABLE]

The next step consists of showing that \mathrm{Tr}\big{(}\tilde{\bm{G}}_{[-i]}^{-1}\big{)} (resp. $\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{G}}_{[-i]}^{-1}\tilde{\bm{X}}_{i}$) and \mathrm{Tr}\big{(}\tilde{\bm{G}}_{(i)}^{-1}\big{)} (resp. $\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{G}}_{(i)}^{-1}\tilde{\bm{X}}_{i}$) are uniformly close across all $i$. This is established in the following lemma.

Lemma 15**.**

Let $\tilde{\bm{G}}_{(i)}$ and $\tilde{\bm{G}}_{[-i]}$ be defined as in (82) and (84), respectively. Then there exist universal constants $C_{1},C_{2},C_{3},C_{4},c_{2},c_{3}>0$ such that

[TABLE]

[TABLE]

where $K_{n},H_{n}$ are diverging sequences as defined in (134).

This together with (145) yields that

[TABLE]

The final ingredient is to establish that \frac{1}{n}\mathrm{Tr}\big{(}\tilde{\bm{G}}_{(i)}^{-1}\big{)} and \frac{1}{n}\mathrm{Tr}\big{(}\tilde{\bm{G}}^{-1}\big{)} are uniformly close across $i$.

Lemma 16**.**

*Let $\tilde{\bm{G}}$ and $\tilde{\bm{G}}_{(i)}$ be as defined in (81) and (82), respectively. Then one has *

[TABLE]

This completes the proof.

∎

• Proof of Lemma 15:

For two invertible matrices $\bm{A}$ and $\bm{B}$ of the same dimensions, the difference of their inverses can be written as

[TABLE]

Applying this identity, we have

[TABLE]

From the definition of these matrices, it follows directly that

[TABLE]

As $\rho^{\prime\prime\prime}$ is bounded, by the mean-value theorem, it suffices to control the differences $\bm{X}_{j}^{\top}\tilde{\bm{\beta}}_{[-i]}-\tilde{\bm{X}}_{j}^{\top}\tilde{\bm{\beta}}$ uniformly across all $j$. This is established in the following lemma, the proof of which is deferred to Appendix H.

Lemma 17**.**

Let $\hat{\bm{\beta}}$ be the full model MLE and $\hat{\bm{\beta}}_{[-i]}$ be the MLE when the $i{\text{th}}$ observation is dropped. Let $q_{i}$ be as described in Lemma 18 and $K_{n},H_{n}$ be as in (134). Then there exist universal constants $C_{1},C_{2},C_{3},C_{4},c_{2},c_{3}>0$ such that

[TABLE]

[TABLE]

Invoking this lemma, we see that the spectral norm of (150) is bounded above by some constant times

[TABLE]

with high probability as specified in (151). From Lemma 2, the spectral norm here is bounded by some constant with probability at least $1-c_{1}\exp(-c_{2}n)$. These observations together with (86) and the fact that on $\mathcal{A}_{n}$ the minimum eigenvalues of $\tilde{\bm{G}}_{(i)}$ and $\tilde{\bm{G}}_{[-i]}$ are bounded by $\lambda_{\mathrm{lb}}$ yield that

[TABLE]

This is true for any $i$. Hence, taking the union bound we obtain

[TABLE]

In order to establish the first result, note that

[TABLE]

To obtain the second result, note that

[TABLE]

Therefore, combining (153) and Lemma 2 gives the desired result.

∎

• Proof of Lemma 16:

We restrict ourselves to the event $\mathcal{A}_{n}$ throughout. Recalling (137), one has

[TABLE]

In addition, on $\mathcal{A}_{n}$ we have

[TABLE]

Combining these results and recognizing that $\rho^{\prime\prime}>0$, we get

[TABLE]

as claimed.

∎

Appendix G Proof of Lemma 11

Again, we restrict ourselves to the event $\mathcal{A}_{n}$ on which $\tilde{\bm{G}}\succeq\lambda_{\mathrm{lb}}\bm{I}$. Note that

[TABLE]

Note that $\{\tilde{\bm{G}},\tilde{\bm{X}}\}$ and $\bm{X}_{\cdot 1}$ are independent. Conditional on $\tilde{\bm{X}}$, the left-hand side is Gaussian with mean zero and variance $\frac{1}{n^{2}}\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{G}}^{-1}\tilde{\bm{X}}^{\top}\bm{D}_{\tilde{\bm{\beta}}}^{2}\tilde{\bm{X}}\tilde{\bm{G}}^{-1}\tilde{\bm{X}}_{i}$. The variance is bounded above by

[TABLE]

In turn, Lemma 2 asserts that $n^{-1}\|\tilde{\bm{X}}_{i}\|^{2}$ is bounded by a constant with high probability. As a result, applying Gaussian concentration results [49, Theorem 2.1.12] gives

[TABLE]

with probability exceeding $1-C\exp\left(-cH_{n}^{2}\right)$, where $C,c>0$ are universal constants.

In addition, $\sup_{i}|X_{i1}|\lesssim H_{n}$ holds with probability exceeding $1-C\exp\left(-cH_{n}^{2}\right)$. Putting the above results together, applying the triangle inequality $|X_{i1}-\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{G}}^{-1}\bm{w}|\leq|X_{i1}|+|\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{G}}^{-1}\bm{w}|$, and taking the union bound, we obtain

[TABLE]

Appendix H Proof of Lemma 17

The goal of this section is to prove Lemma 17, which relates the full-model MLE $\hat{\bm{\beta}}$ and the MLE $\hat{\bm{\beta}}_{[-i]}$. To this end, we establish the key lemma below.

Lemma 18**.**

Suppose $\hat{\bm{\beta}}_{[-i]}$ denote the MLE when the $i{\text{th}}$ observation is dropped. Further let $\bm{G}_{[-i]}$ be as in (83), and define $q_{i}$ and $\hat{\bm{b}}$ as follows:

[TABLE]

Suppose $K_{n},H_{n}$ are diverging sequences as in (134). Then there exist universal constants $C_{1},C_{2},C_{3}>0$ such that

[TABLE]

[TABLE]

The proof ideas are inspired by the leave-one-observation-out approach of [20]. We however emphasize once more that the adaptation of these ideas to our setup is not straightforward and crucially hinges on Theorem 4, Lemma 7 and properties of the effective link function.

• Proof of Lemma 18:

Invoking techniques similar to that for establishing Lemma 7, it can be shown that

[TABLE]

with probability at least $1-\exp(\Omega(n))$, where $\gamma_{i}^{*}$ is between $\bm{X}_{i}^{\top}\hat{\bm{b}}$ and $\bm{X}_{i}^{\top}\hat{\bm{\beta}}$. Denote by $\mathcal{B}_{n}$ the event where (159) holds. Throughout this proof, we work on the event $\mathcal{C}_{n}:=\mathcal{A}_{n}\cap\mathcal{B}_{n}$, which has probability $1-\exp\left(-\Omega(n)\right)$. As in (106) then,

[TABLE]

Next, we simplify (160). To this end, recall the defining relation of the proximal operator

[TABLE]

which together with the definitions of $\hat{\bm{b}}$ and $q_{i}$ gives

[TABLE]

Now, let $\ell_{[-i]}$ denote the negative log-likelihood function when the $i{\text{th}}$ observation is dropped, and hence \nabla\ell_{[-i]}\big{(}\hat{\bm{\beta}}_{[-i]}\big{)}=\bm{0}. Expressing $\nabla\ell(\hat{\bm{b}})$ as \nabla\ell(\hat{\bm{b}})-\nabla\ell_{[-i]}\big{(}\hat{\bm{\beta}}_{[-i]}\big{)}, applying the mean value theorem, and using the analysis similar to that in [20, Proposition 3.4], we obtain

[TABLE]

where $\gamma_{j}^{*}$ is between $\bm{X}_{j}^{\top}\hat{\bm{b}}$ and $\bm{X}_{j}^{\top}\hat{\bm{\beta}}_{[-i]}$. Combining (160) and (162) leads to the upper bound

[TABLE]

We need to control each term in the right-hand side. To start with, the first term is bounded by a universal constant with probability $1-\exp(-\Omega(n))$ (Lemma 2). For the second term, since $\gamma_{j}^{*}$ is between $\bm{X}_{j}^{\top}\hat{\bm{b}}$ and $\bm{X}_{j}^{\top}\hat{\bm{\beta}}_{[-i]}$ and $\|\rho^{\prime\prime\prime}\|_{\infty}<\infty$, we get

[TABLE]

Given that $\{\bm{X}_{j},\bm{G}_{[-i]}\}$ and $\bm{X}_{i}$ are independent for all $j\neq i$, conditional on $\{\bm{X}_{j},\bm{G}_{[-i]}\}$ one has

[TABLE]

In addition, the variance satisfies

[TABLE]

with probability at least $1-\exp(-\Omega(n))$. Applying standard Gaussian concentration results [49, Theorem 2.1.12], we obtain

[TABLE]

By the union bound

[TABLE]

Consequently,

[TABLE]

In addition, the third term in the right-hand side of (163) can be upper bounded as well since

[TABLE]

with high probability.

It remains to bound $\left|\rho^{\prime}\left(\mathsf{prox}_{q_{i}\rho}(\bm{X}_{i}^{\top}\hat{\bm{\beta}}_{[-i]})\right)\right|$. To do this, we begin by considering $\rho^{\prime}(\mathsf{prox}_{c\rho}(Z))$ for any constant $c>0$ (rather than a random variable $q_{i}$). Recall that for any constant $c>0$ and any $Z\sim{\mathcal{N}}(0,\sigma^{2})$ with finite variance, the random variable $\rho^{\prime}(\mathsf{prox}_{c\rho}(Z))$ is sub-Gaussian. Conditional on $\hat{\bm{\beta}}_{[-i]}$, one has \bm{X}_{i}^{\top}\hat{\bm{\beta}}_{[-i]}\sim{\mathcal{N}}\big{(}0,\|\hat{\bm{\beta}}_{[-i]}\|^{2}\big{)}. This yields

[TABLE]

for some constants $C_{1},C_{2},C_{3},C_{4},C_{5}>0$ since$\|\hat{\bm{\beta}}_{[-i]}\|$ is bounded with high probability (see Theorem 4).

Note that $\frac{\partial\mathsf{prox}_{b\rho}(z)}{\partial b}\leq 0$ by [18, Proposition 6.3]. Hence, in order to move over from the above concentration result established for a fixed constant $c$ to the random variables $q_{i}$, it suffices to establish a uniform lower bound for $q_{i}$ with high probability. Observe that for each $i$,

[TABLE]

with probability $1-\exp(-\Omega(n))$, where $C^{*}$ is some universal constant. On this event, one has

[TABLE]

This taken collectively with (Proof of Lemma 18:) yields

[TABLE]

This controls the last term.

To summarize, if $\{K_{n}\}$ and $\{H_{n}\}$ are diverging sequences satisfying the assumptions in (134), combining (163) and the bounds for each term in the right-hand side finally gives (157). On the other hand, combining (169) and (174) yields (158).

∎

With the help of Lemma 18 we are ready to prove Lemma 17. Indeed, observe that

[TABLE]

and hence by combining Lemma 2 and Lemma 18, we establish the first claim (151). The second claim (152) follows directly from Lemmas 2, 18 and (161).

Appendix I Proof of Theorem 7(b)

This section proves that the random sequence \tilde{\alpha}=\mathrm{Tr}\big{(}\tilde{\bm{G}}^{-1}\big{)}/n converges in probability to the constant $b_{\ast}$ defined by the system of equations (23) and (24). To begin with, we claim that $\tilde{\alpha}$ is close to a set of auxiliary random variables $\{\tilde{q}_{i}\}$ defined below.

Lemma 19**.**

Define $\tilde{q}_{i}$ to be

[TABLE]

where $\tilde{\bm{G}}_{[-i]}$ is defined in (84). Then there exist universal constants $C_{1},C_{2},C_{3},C_{4},c_{2},c_{3}>0$ such that

[TABLE]

where $K_{n},H_{n}$ are as in (134).

• Proof:

This result follows directly from Proposition 1 and equation (146).

∎

A consequence is that $\mathsf{prox}_{\tilde{q}_{i}\rho}\left(\bm{X}_{i}^{\top}\hat{\bm{\beta}}_{[-i]}\right)$ becomes close to $\mathsf{prox}_{\tilde{\alpha}\rho}\left(\bm{X}_{i}^{\top}\hat{\bm{\beta}}_{[-i]}\right)$.

Lemma 20**.**

Let $\tilde{q}_{i}$ and $\tilde{\alpha}$ be as defined earlier. Then one has

[TABLE]

where $K_{n},H_{n}$ are as in (134).

The key idea behind studying $\mathsf{prox}_{\tilde{\alpha}\rho}\left(\bm{X}_{i}^{\top}\hat{\bm{\beta}}_{[-i]}\right)$ is that it is connected to a random function $\delta_{n}(\cdot)$ defined below, which happens to be closely related to the equation (24). In fact, we will show that $\delta_{n}(\tilde{\alpha})$ converges in probability to 0; the proof relies on the connection between $\mathsf{prox}_{\tilde{\alpha}\rho}\left(\bm{X}_{i}^{\top}\hat{\bm{\beta}}_{[-i]}\right)$ and the auxiliary quantity $\mathsf{prox}_{\tilde{q}_{i}\rho}\left(\bm{X}_{i}^{\top}\hat{\bm{\beta}}_{[-i]}\right)$. The formal results is this:

Proposition 2**.**

For any index $i$, let $\hat{\bm{\beta}}_{[-i]}$ be the MLE obtained on dropping the $i{\text{th}}$ observation. Define $\delta_{n}(x)$ to be the random function

[TABLE]

Then one has $\delta_{n}(\tilde{\alpha})\leavevmode\nobreak\ \stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}\leavevmode\nobreak\ 0$.

Furthermore, the random function $\delta_{n}(x)$ converges to a deterministic function $\Delta(x)$ defined by

[TABLE]

where $Z\sim{\mathcal{N}}(0,1)$, and $\tau_{\ast}$ is such that $(\tau_{\ast},b_{\ast})$ is the unique solution to (23) and (24).

Proposition 3**.**

With $\Delta(x)$ as in (177), $\Delta(\tilde{\alpha})\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}0.$

In fact, one can easily verify that

[TABLE]

and hence by Lemma 5, the solution to $\Delta(x)=0$ is exactly $b_{\ast}$. As a result, putting the above claims together, we show that $\tilde{\alpha}$ converges in probability to $b_{\ast}$.

It remains to formally prove the preceding lemmas and propositions, which is the goal of the rest of this section.

• Proof of Lemma 20:

By [18, Proposition 6.3], one has

[TABLE]

which yields

[TABLE]

where $q_{\tilde{\alpha},i}$ is between $\tilde{q}_{i}$ and $\tilde{\alpha}$. Here, the last inequality holds since $q_{\tilde{\alpha},i},\rho^{\prime\prime}\geq 0$.

In addition, just as in the proof of Lemma 18, one can show that $q_{i}$ is bounded below by some constant $C^{*}>0$ with probability $1-\exp(-\Omega(n))$. Since $q_{\tilde{\alpha},i}\geq\min\{\tilde{q}_{i},\tilde{\alpha}\}$, on the event $\sup_{i}|\tilde{q}_{i}-\tilde{\alpha}|\leq C_{1}K_{n}^{2}H_{n}/\sqrt{n}$, which happens with high probability (Lemma 19), $q_{\tilde{\alpha},i}\geq C_{\alpha}$ for some universal constant $C_{\alpha}>0$. Hence, by an argument similar to that establishing (174), we have

[TABLE]

This together with (179) and Lemma 19 concludes the proof.

∎

• Proof of Proposition 2:

To begin with, recall from (140) and (141) that on $\mathcal{A}_{n}$,

[TABLE]

Using the fact that \big{|}\frac{1}{1+x}-\frac{1}{1+y}\big{|}\leq|x-y| for $x,y\geq 0$, we obtain

[TABLE]

with high probability (Proposition 1). This combined with (180) yields

[TABLE]

The above bound concerns $\frac{1}{n}\sum_{i=1}^{n}\frac{1}{1+\rho^{\prime\prime}(\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{\beta}})\tilde{\alpha}}$, and it remains to relate it to $\frac{1}{n}\sum_{i=1}^{n}\frac{1}{1+\rho^{\prime\prime}\left(\mathsf{prox}_{\tilde{\alpha}\rho}\left(\tilde{\bm{X}}_{i}^{\top}\tilde{\bm{\beta}}\right)\right)\tilde{\alpha}}$. To this end, we first get from the uniform boundedness of $\rho^{\prime\prime\prime}$ and Lemma 17 that

[TABLE]

Note that

[TABLE]

By the bound (181), an application of Lemma 20, and the fact that $\tilde{\alpha}\leq p/(n\lambda_{\mathrm{lb}})$ (on $\mathcal{A}_{n}$), we obtain

[TABLE]

This establishes that $\delta_{n}(\tilde{\alpha})\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}0$.

∎

• Proof of Proposition 3:

Note that since $0<\alpha\leq p/(n\lambda_{\mathrm{lb}}):=B$ on $\mathcal{A}_{n}$, it suffices to show that

[TABLE]

We do this by following three steps. Below, $M>0$ is some sufficiently large constant.

First we truncate the random function $\delta_{n}(x)$ and define

[TABLE]

The first step is to show that $\sup_{x\in[0,B]}\left|\tilde{\delta}_{n}(x)-\delta_{n}(x)\right|\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}0.$ We stress that this truncation does not arise in [20], and we keep track of the truncation throughout the rest of the proof. 2. 2.

Show that \sup_{x\in[0,B]}\left|\tilde{\delta}_{n}(x)-\operatorname{\mathbb{E}}\big{[}\tilde{\delta}_{n}(x)\big{]}\right|\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}0. 3. 3.

Show that \sup_{x\in[0,B]}\left|\operatorname{\mathbb{E}}\big{[}\tilde{\delta}_{n}(x)\big{]}-\Delta(x)\right|\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}0.

To argue about the first step, observe that $|\frac{1}{1+y}-\frac{1}{1+z}|\leq|y-z|$ for any $y,z>0$ and that \big{|}\frac{\partial\mathsf{prox}_{c\rho}(x)}{\partial x}\big{|}\leq 1 [18, Proposition 6.3]. Then

[TABLE]

For a sufficiently large constant $M>0$, $\mathbb{P}(\|\hat{\bm{\beta}}_{[-i]}\|\geq M)\leq\exp(-\Omega(n))$ by Theorem 4. Hence, for any $\epsilon>0$,

[TABLE]

establishing Step 1.

To argue about Step 2, note that for any $x$ and $z$,

[TABLE]

where $\tilde{z}\in(x,z).$ Setting

[TABLE]

then for any $\epsilon,\eta>0$ we have

[TABLE]

where $C_{1}(\epsilon)$ is some function independent of $n$. The inequality (183) is an analogue of [20, Lemma 3.24]. We remark that the truncation is particularly important here in guaranteeing that $\operatorname{\mathbb{E}}[Y_{n}]<\infty$.

Set

[TABLE]

and observe that

[TABLE]

A similarly inequality applies to $\Delta(x)$ in which $\bm{X}_{i}^{\top}\hat{\bm{\beta}}_{[-i]}\bm{1}_{\{\|\hat{\bm{\beta}}_{[-i]}\|\leq M\}}$ is replaced by $\tau_{\ast}Z$. In either case,

[TABLE]

for any $\eta>0$.

For any $\epsilon^{\prime}>0$, set $K=\max\{C_{1}(\epsilon^{\prime}),C_{2}\}$. Next, divide $[0,B]$ into finitely many segments

[TABLE]

such that the length of each segment is $\eta/K$ for any $\eta>0$. Then for every $x\in[0,B]$, there exists $l$ such that $|x-x_{l}|\leq\eta/K$. As a result, for any $x\in[0,B]$,

[TABLE]

Now fix $\delta>0,\epsilon>0$. Applying the above inequality gives

[TABLE]

Choose $\eta<\min\{\epsilon/2,\delta\},K=\max\{C_{1}(\frac{\delta-\eta}{2}),C_{2}\}$. Then the first term in the right-hand side is at most $\epsilon/2$. Furthermore, suppose one can establish that for any fixed $x$,

[TABLE]

Since the second term in the right-hand side is a supremum over finitely many points, there exists an integer $N$ such that for all $n\geq N$, the second term is less than or equal to $\epsilon/2$. Hence, for all $n\geq N$, right-hand side is at most $\epsilon$, which proves Step $2$. Thus, it remains to prove (186) for any fixed $x\geq 0$. We do this after the analysis of Step 3.

We argue about Step 3 in a similar fashion and letting $K=\max\{C_{2},C_{3}\}$, divide $(0,B]$ into segments of length $\eta/K$. For every $x\in(0,B]$ there exists $x_{l}$ such that $|x-x_{l}|\leq\eta/K$. Then

[TABLE]

[TABLE]

Hence, it suffices to show that for any fixed $x$, $|G_{n}(x)-\Delta(x)|\rightarrow 0.$ To this end, observe that for any fixed $x$,

[TABLE]

Additionally,

[TABLE]

where $\tilde{Z}\sim{\mathcal{N}}(0,1)$. Since $\|\hat{\bm{\beta}}_{[-1]}\|\bm{1}_{\{\|\hat{\bm{\beta}}_{[-1]}\|\leq M\}}\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}\tau_{\ast}$, by Slutsky’s theorem, $\bm{X}_{1}^{\top}\hat{\bm{\beta}}_{[-1]}\bm{1}_{\{\|\hat{\bm{\beta}}_{[-1]}\|\leq M\}}$ converges weakly to $\tau_{\ast}\tilde{Z}$. Since $t\mapsto 1/(1+x\rho"(\mathsf{prox}_{x\rho}(t)))$ is bounded, one directly gets

[TABLE]

for every $x$.

Finally, we establish (186). To this end, we will prove instead that for any given $x$,

[TABLE]

Define

[TABLE]

Then $\delta_{n}(x)-G_{n}(x)=\frac{1}{n}\sum_{i=1}^{n}f(M_{i})$. Hence, for any $x\in[0,B]$,

[TABLE]

The first term in the right-hand side is at most $1/n$. Hence, it suffices to show that \operatorname{\mathbb{E}}\big{[}f(M_{1})f(M_{2})\big{]}\rightarrow 0. Let $\hat{\bm{\beta}}_{[-12]}$ be the MLE when the $1{\text{st}}$ and $2{\text{nd}}$ observations are dropped. Define

[TABLE]

By an application of Lemma 18,

[TABLE]

Also, by the triangle inequality,

[TABLE]

Invoking Lemma 2, (187), and an argument similar to that leading to (174) and (168), we obtain

[TABLE]

The event $\{\|\hat{\bm{\beta}}_{[-1]}\|\leq M\}\cap\{\|\hat{\bm{\beta}}_{[-12]}\|\leq M\}$ occurs with probability at least $1-C\exp(-cn)$. Hence, one obtains

[TABLE]

A similar statement continues to hold with $\bm{X}_{1}$ replaced by $\bm{X}_{2}$ and $\hat{\bm{\beta}}_{[-1]}$ replaced by $\hat{\bm{\beta}}_{[-2]}$. Some simple computation yields that $\|f^{\prime}\|_{\infty}$ is bounded by some constant times $|x|$. By the mean value theorem and the fact that $\|f\|_{\infty}\leq 1,$

[TABLE]

Consequently,

[TABLE]

As $\|f\|_{\infty}\leq 1$, this implies convergence in $\mathrm{L}_{1}$. Thus, it simply suffices to show that

[TABLE]

Denote the design matrix on dropping the first and second row as $\bm{X}_{[-12]}$. Note that conditional on $\bm{X}_{[-12]}$, $\bm{X}_{1}^{\top}\hat{\bm{\beta}}_{[-12]}\bm{1}_{\{\|\hat{\bm{\beta}}_{[-12]}\|\leq M\}}$ and $\bm{X}_{2}^{\top}\hat{\bm{\beta}}_{[-12]}\bm{1}_{\{}\|\hat{\bm{\beta}}_{[-12]}\|\leq M\}$ are independent and have distribution

[TABLE]

Using this and by arguments similar to [20, Lemma 3.23], one can show that

[TABLE]

On repeated application of the multivariate inversion theorem for obtaining densities from characteristic functions, we get that

[TABLE]

Since $f$ is centered, this completes the proof.

∎