The existence of maximum likelihood estimate in high-dimensional binary response generalized linear models

TL;DR

This paper investigates the conditions under which maximum likelihood estimates exist in high-dimensional binary response generalized linear models, revealing a phase transition phenomenon influenced by covariate distribution.

Contribution

It extends previous work by establishing a phase transition for a broad class of models beyond Gaussian covariates, using geometric and probabilistic tools.

Findings

Existence of MLE exhibits a phase transition depending on model parameters.

The phase transition is influenced by covariate distribution, not just Gaussianity.

Simulation studies support theoretical phase transition results.

Abstract

Motivated by recent works on the high-dimensional logistic regression, we establish that the existence of the maximum likelihood estimate exhibits a phase transition for a wide range of generalized linear models with binary outcome and elliptical covariates. This extends a previous result of Cand\`es and Sur who proved the phase transition for the logistic regression with Gaussian covariates. Our result reveals a rich structure in the phase transition phenomenon, which is simply overlooked by Gaussianity. The main tools for deriving the result are data separation, convex geometry and stochastic approximation. We also conduct simulation studies to corroborate our theoretical findings, and explore other features of the problem.

Tables

Table 1. Table 1 . Summary of the theoretical h MLE subscript ℎ MLE h_{\text{MLE}} and the simulations of the phase transition for the maximum likelihood estimate with the logit link function. h 0.5 subscript ℎ 0.5 h_{0.5} is the κ 𝜅 \kappa such that the proportion of times that the maximum likelihood estimate exists is 0.5 0.5 0.5 , and the number inside the bracket is the width of the uncertainty interval (a slice of the uncertainty band ) for a given γ 0 subscript 𝛾 0 \gamma_{0} . MIW is the mean width of the uncertainty intervals across γ 0 ∈ ( 0 , 10 ] subscript 𝛾 0 0 10 \gamma_{0}\in(0,10] ; MD is the mean difference between h MLE subscript ℎ MLE h_{\text{MLE}} and h 0.5 subscript ℎ 0.5 h_{0.5} .

	${\gamma }_0=1$		${\gamma }_0=9$		Overall
Distribution	$h_{0.5}$	$h_{\text{MLE}}$	$h_{0.5}$	$h_{\text{MLE}}$	MIW	MD
$Gamma(0.5,{\theta }_{0.5})$	0.435 (0.075)	0.4238	0.310 (0.095)	0.310	0.101	0.0077
$Gamma(1,{\theta }_1)$	0.450 (0.050)	0.458	0.260 (0.080)	0.246	0.071	0.0186
$Gamma(2,{\theta }_2)$	0.450 (0.050)	0.447	0.205 (0.070)	0.191	0.057	0.0160
$Pareto(2.5,x_m\left(2.5\right))$	0.455 (0.045)	0.458	0.165 (0.065)	0.172	0.055	0.0045
$Pareto(3.5,x_m\left(3.5\right))$	0.440 (0.045)	0.439	0.120 (0.060)	0.137	0.057	0.0095
$Pareto(4.5,x_m\left(4.5\right))$	0.435 (0.055)	0.430	0.110 (0.060)	0.128	0.055	0.0095
half-normal	0.450 (0.050)	0.448	0.240 (0.075)	0.211	0.064	0.0215
log-normal	0.380 (0.135)	0.497	0.355 (0.150)	0.450	0.152	0.1199

Equations

$$\mathbb{P}(y_i=1|{\mathbf{x}}_i)=\sigma ({\beta }_0+{\mathbf{x}}_i^T\mathbf{\beta }).$$

$$\ell ({\beta }_0,\mathbf{\beta }):=\sum _{y_i=1}\log (\sigma ({\beta }_0+{\mathbf{x}}_i^T\mathbf{\beta }))+\sum _{y_i=-1}\log (1-\sigma ({\beta }_0+{\mathbf{x}}_i^T\mathbf{\beta })).$$

$$\mathbf{X}\stackrel{(d)}{=}\mathbf{\mu }+R\mathbf{A}\mathbf{U},$$

$${\mathbf{\mu }}_{1|2}:={\mathbf{\mu }}_1+{\mathbf{\Sigma }}_{12}{\mathbf{\Sigma }}_{22}^{-1}({\mathbf{x}}_2-{\mathbf{\mu }}_2),{\mathbf{\Sigma }}_{1|2}:={\mathbf{\Sigma }}_{11}-{\mathbf{\Sigma }}_{12}{\mathbf{\Sigma }}_{22}^{-1}{\mathbf{\Sigma }}_{21},$$

$$F_{R_{1|2}}(r)=\frac{{\int }_{d_{{\mathbf{\Sigma }}_{22}}({\mathbf{x}}_2,{\mathbf{\mu }}_2)}^{\sqrt{r^2+d_{{\mathbf{\Sigma }}_{22}}^2({\mathbf{x}}_2,{\mathbf{\mu }}_2)}}{\left(s^2-d_{{\mathbf{\Sigma }}_{22}}^2({\mathbf{x}}_2,{\mathbf{\mu }}_2)\right)}^{p_1/2-1}{s}^{-p+2}d{F}_R(s)}{{\int }_{d_{{\mathbf{\Sigma }}_{22}}({\mathbf{x}}_2,{\mathbf{\mu }}_2)}^{\infty }{\left(s^2-d_{{\mathbf{\Sigma }}_{22}}^2({\mathbf{x}}_2,{\mathbf{\mu }}_2)\right)}^{p_1/2-1}{s}^{-p+2}d{F}_R(s)},$$

$$(Y^{(p)},X^{(p)})\sim F_{{\alpha }_0,{\beta }_0,{\gamma }_0}\text{mbox}if(Y^{(p)},X^{(p)})\stackrel{(d)}{=}(V^{(p)},V^{(p)}{U}^{(p)}),$$

$$\sum _{k=1}^{\infty }{m}_{2k}^{-\frac{1}{2k}}=\infty ,$$

$$p_{+}(x)=\sigma \left({\beta }_0+\frac{{\gamma }_0}{{\alpha }_0}x\right)\text{mbox}and{p}_{-}(x):=1-p_{+}(x).$$

$$G_{p,+}(x)={\int }_{z\le x}{p}_{+}(z)f_{X^{(p)}}(z)dz\text{mbox}and{G}_{p,-}(x)={\int }_{z\le x}{p}_{-}(z)f_{X^{(p)}}(z)dz,$$

$${\overline{G}}_{p,+}(x)={\int }_{z>x}{p}_{+}(z)f_{X^{(p)}}(z)dz\text{mbox}and{\overline{G}}_{p,-}(x)={\int }_{z>x}{p}_{-}(z)f_{X^{(p)}}(z)dz.$$

$$h_{\text{MLE}}({\alpha }_0,{\beta }_0,{\gamma }_0):$$

$$=\underset{{\lambda }_0,{\lambda }_1\in \mathbb{R}}{\min }\mathbb{E}({\lambda }_{0}Y+{\lambda }_{1}X-Z{)}_{+}^2,$$

$$\mathbb{E}[p_{\pm }(X^{(p)})(G_{p,\mp }(X^{(p)})+{\overline{G}}_{p,\pm }(X^{(p)}){)}^{n-1}]=o\left(\frac{1}{n}\right),$$

$$\kappa >h_{\text{MLE}}({\alpha }_0,{\beta }_0,{\gamma }_0)\Rightarrow \underset{n,p\to \infty }{\lim }\mathbb{P}(\text{mbox}maximumlikelihoodestimateexists)=0.$$

$$\kappa <h_{\text{MLE}}({\alpha }_0,{\beta }_0,{\gamma }_0)\Rightarrow \underset{n,p\to \infty }{\lim }\mathbb{P}(\text{mbox}maximumlikelihoodestimateexists)=1.$$

$$G_{p,\mp }(x)+{\overline{G}}_{p,\pm }(x)<1-ϵ\text{mbox}forallx.$$

$$G_{p,-}(x)+{\overline{G}}_{p,+}(x)=\frac{1}{2}+{\int }_x^{\infty }\frac{e^z-1}{1+e^z}{f}_X(z)dz<1.$$

$$\underset{b_0,\mathbf{b}}{\max }$$

$$y_i(b_0+{\mathbf{x}}_i^T\mathbf{b})\ge 0,i=1,\dots ,n$$

$$\mathbb{E}{R}^2=p{\alpha }_0^2+1\text{mbox}and\mathbf{\beta }=(\tilde{\mathbf{W}}/||\tilde{\mathbf{W}}|{|}_2+1/p)\cdot {\gamma }_0/{\alpha }_0,$$

$$(y_i,y_i{\mathbf{x}}_i)\stackrel{(d)}{=}(Y^{(p)},X^{(p)},X_2,\dots ,X_p),$$

$$(X_2,\dots ,X_p|X^{(p)}=x)\sim {\mathcal{E}}_{p-1}(0,{\mathbf{I}}_{p-1},F_{R_{-1}}),$$

$$\mathcal{C}(\mathcal{W})=\{\mathbf{w}+\mathbf{u}:\mathbf{w}\in \mathcal{W},\mathbf{u}\ge 0\},$$

$$\mathcal{L}:=\mathrm{span}({\mathbf{X}}_2,\dots ,{\mathbf{X}}_p)\text{mbox}and\mathcal{W}:=\mathrm{span}({\mathbf{Y}}^{(p)},{\mathbf{X}}^{(p)}).$$

$$0\le \mathbb{P}(\text{mbox}noMLE)-\mathbb{P}(\mathcal{L}\cap \mathcal{C}(\mathcal{W})\ne \{0\})\le \mathbb{P}(\text{mbox}NoMLESingle).$$

$$p-1+\delta (\mathcal{C}(\mathcal{W}))>n+a_{\epsilon }\sqrt{n}⟹\mathbb{P}(\mathcal{L}\cap \mathcal{C}(\mathcal{W}))\ge 1-\epsilon ,$$

$$p-1+\delta (\mathcal{C}(\mathcal{W}))<n-a_{\epsilon }\sqrt{n}⟹\mathbb{P}(\mathcal{L}\cap \mathcal{C}(\mathcal{W}))\le \epsilon ,$$

$$\delta (\mathcal{C}(\mathcal{W}))=n-\mathbb{E}\left(\underset{\mathbf{w}\in \mathcal{W}}{\min }|(\mathbf{w}-\mathbf{Z}{)}_{+}{|}^2\right).$$

$$Q_{p,n}:=\underset{{\lambda }_0,{\lambda }_1\in \mathbb{R}}{\min }\frac{1}{n}|({\lambda }_0{\mathbf{Y}}^{(p)}+{\lambda }_1{\mathbf{X}}^{(p)}-\mathbf{Z}{)}_{+}{|}^2.$$

$$\underset{\mathbf{\lambda }\in S}{\min }G(\mathbf{\lambda }),G(\mathbf{\lambda }):=\mathbb{E}g(\mathbf{\lambda },\mathbf{\xi }),$$

$$\underset{\mathbf{\lambda }\in S}{\min }{\hat{G}}_n(\mathbf{\lambda };{\mathbf{\xi }}_{1,n},\dots ,{\mathbf{\xi }}_{n,n}),{\hat{G}}_n(\mathbf{\lambda };{\mathbf{\xi }}_{1,n},\dots ,{\mathbf{\xi }}_{n,n}):=\frac{1}{n}\sum _{i=1}^{n}g(\mathbf{\lambda },{\mathbf{\xi }}_{i,n}),$$

$${\hat{v}}_n\to v\text{mbox}inprobability.$$

$$\mathbb{P}(\mathcal{L}\cap \mathcal{C}(\mathcal{W})\ne \{0\})\ge \mathbb{P}(p/n>\mathbb{E}(Q_{p.n}|\mathbf{X},\mathbf{Y})+a_n{n}^{-1/2})-{\epsilon }_n,$$

Full text

1. Introduction

In this paper, we are concerned with the maximum likelihood estimate of generalized linear models [32, 34] with binary outcome. More precisely, we consider $n$ independent and identically distributed observations $(\boldsymbol{x}_{i},y_{i})$, $i=1,\ldots,n$, where the binary outcome $y_{i}\in\{-1,1\}$ is connected to the covariates $\boldsymbol{x}_{i}\in\mathbb{R}^{p}$ by the probability model

[TABLE]

Here $\sigma:\mathbb{R}\rightarrow[0,1]$ is the inverse link function, and $(\beta_{0},\boldsymbol{\beta})\in\mathbb{R}^{p+1}$ are unknown parameters of the model. Popular choices for $\sigma$ are

• •

logit link function [7]: $\sigma(t)=e^{t}/(1+e^{t})$.

• •

probit link function [9]: $\sigma(t)=\Phi(t)$, where $\Phi(t):=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{t}\exp(-s^{2}/2)ds$.

• •

cloglog link function [21]: $\sigma(t)=1-e^{-e^{t}}$.

The maximum likelihood estimate of $(\beta_{0},\boldsymbol{\beta})$ is any maximizer of the log-likelihood

[TABLE]

In contrast to the maximum likelihood estimate of linear models, that of generalized linear models does not always exist. This phenomenon is closely related to the separability of observed data, see Section 2.1 for a review. Classical theory deals with this issue when the number of covariates $p$ is fixed, and the number of observed data $n$ tends to infinity. In the era of data deluge, we are often in a situation where the number of covariates $p$ and the number of observations $n$ are comparable in size. The problems of interest are in high-dimensional asymptotics, in which case the number of parameters $p$ and the number of observations $n$ both tend to infinity, at the same rate.

In a series of papers [12, 39, 40], Sur, Chen and Candès developed a theory for the logistic regression with Gaussian covariates in high-dimensional regimes. They studied the asymptotic properties of the maximum likelihood estimate when $p/n\rightarrow\kappa$, with applications in hypothesis testing. Candès and Sur [12] proved a phase transition for the existence of the maximum likelihood estimate in high-dimensional logistic regression with Gaussian covariates. This extends an earlier result of Cover [14] in the context of information theory. Formally, there exists a threshold $h_{\text{MLE}}$, depending on the parameters of the model, such that

• •

if $\kappa>h_{\text{MLE}}$, then $\mathbb{P}({\mbox{maximum likelihood estimate exists}})\rightarrow 0$ as $n,p\rightarrow\infty$.

• •

if $\kappa<h_{\text{MLE}}$, then $\mathbb{P}({\mbox{maximum likelihood estimate exists}})\rightarrow 1$ as $n,p\rightarrow\infty$.

This phenomenon is referred to as the phase transition for the existence of the maximum likelihood estimate. The latter is crucial to justify the use of large sample approximations to numerous measures of goodness-of-fit, and derive the limiting distribution of the likelihood ratio, as mentioned in [12]. But they only studied the existence of the maximum likelihood estimate for Gaussian covariates. This is rarely the case in reality. For instance, the covariates are often heavy-tail distributed in financial problems where $p$ and $n$ are large.

The purpose of this paper is to further generalize the results of [12], proving the phase transition for a large class of generalized linear models with elliptical covariates. Here we consider a large number of covariates sampled from elliptical distributions, and predict whether one can expect the maximum likelihood estimate to be found or not. Elliptical symmetry is a natural generalization of multivariate normality. The contribution of this paper is twofold.

• •

Theoretical justifications. We give a universal threshold on $p/n$ for the existence of the maximum likelihood estimate in the binary classification. Here the word ’universal’ refers to a wide class of link functions and covariate distributions. Our work aims to explore to which extent the phase transition occurs in terms of link functions and covariates, including the logit link and Gaussian covariates as a special case. We notice that the projection limit assumption (Assumption 2.6) is essential to our result, which is disguised for the special choice of Gaussian covariates. This is analogous to regularity assumptions in high-dimensional signal processing [6, 18]. Without this assumption, the phase transition formula may fail (for example, log-normal covariates).

• •

Novel techniques. While the high level idea is the same as [12], we bring a few techniques into this field. First, we give a checkable condition to the projection limit assumption. Second, we use a stochastic approximation approach (Theorem 4.3) to prove the phase transition formula. Compared to the bare-hands analysis in [12], our argument is more general and is easily adapted to other problems of interest. Finally, we provide fairly general conditions (for example, (2.7)) under which the phase transition occurs. These conditions reveal additional structures which are masked by Gaussianity.

In a recent paper of Montanari, Ruan, Sohn and Yan [33], they considered the max-margin classifier of the random feature problem in the high-dimensional regime. They provided a phase transition threshold for the existence of the max-margin classifier, and further studied the limiting max-margin classifier. But similar to [12], they assume that the covariates are Gaussian. De Loera and Hogan [31] considered the maximum likelihood estimate of a multi-class logistic regression in a different way. They explored a condition on the number of observed data and the number of classes such that the maximum likelihood estimate exists. We hope that our work will trigger further research towards a theory for multi-class classification models with non-Gaussian covariates.

The rest of the paper is organized as follows. In Section 2, we provide the background and state the main result, Theorem 2.7. In Section 3, we perform simulations to corroborate our theoretical findings. The proof of Theorem 2.7 is given in Section 4. In Section 5, we conclude with some insights and directions for future work.

2. Background and Main Result

In this section we provide the background on the existence of the maximum likelihood estimate in binary response generalized linear models, and the properties of elliptical distributions. Then we present the main result, Theorem 2.7.

2.1. Existence of the maximum likelihood estimate and Data Geometry

Often the maximum likelihood estimate in the logit model, implemented in many statistical packages, runs smoothly. But sometimes it fails, even when the number of covariates $p$ is much smaller than the sample size $n$. One reason for this undesirable phenomenon is that the maximum likelihood estimate does not exist. It is a classical problem in statistics to characterize the existence and uniqueness of the maximum likelihood estimate in generalized linear models.

Historically, Haberman [22] and Weddenburn [41] provided general criteria for the maximum likelihood estimate to exist. Silvapulle [38], and Albert and Anderson [1] gave conditions for the existence of the maximum likelihood estimate in logistic regression via data geometry. More precisely, they classified the data into the following three categories:

• •

The data points $(\boldsymbol{x}_{i},y_{i})$ are said to be completely separated if there exists $\boldsymbol{b}\in\mathbb{R}^{p}$ such that $y_{i}\boldsymbol{x}_{i}^{T}\boldsymbol{b}>0$ for all $i$.

• •

The data points $(\boldsymbol{x}_{i},y_{i})$ are said to be quasi-completely separated if for each $\boldsymbol{b}\neq\boldsymbol{0}$, $y_{i}\boldsymbol{x}_{i}^{T}\boldsymbol{b}\geq 0$ for all $i$, and equality holds for some $i$.

• •

The data points $(\boldsymbol{x}_{i},y_{i})$ are said to overlap if for each $\boldsymbol{b}\neq\boldsymbol{0}$, there exists one $i$ such that $y_{i}\boldsymbol{x}_{i}^{T}\boldsymbol{b}>0$, and another $i$ such that $y_{i}\boldsymbol{x}_{i}^{T}\boldsymbol{b}<0$.

In [1], it was proved that the maximum likelihood estimate exists in logistic regression if and only if the data points overlap. See also [36] for a generalization. Later Lesaffre and Kaufmann [29] proposed a necessary and sufficient condition for the existence of the maximum likelihood estimate in probit regression, which coincides with that derived in [1] for logistic regression. In fact, their result holds for a general class of generalized linear models.

Theorem 2.1**.**

Consider the generalized linear model defined by (1.1), and assume that $\sigma(\cdot)$ and $1-\sigma(\cdot)$ are log-concave. Then the maximum likelihood estimate exists if and only if the data points overlap.

It is easily seen that the logit, probit and cloglog links all satisfy the log-concavity. Despite this nice characterization, it is not clear how to check these criteria efficiently. See also [13, 16, 28] for algorithmic aspects for detecting separation/overlaps.

2.2. Elliptical Distributions

Elliptical distributions are natural generalizations of multivariate normal, which preserve spherical symmetry. In the sequel, $\mathbb{S}^{p-1}$ denotes the unit sphere in $\mathbb{R}^{p}$. The following definition of elliptical distributions is due to Kelker [25], and Cambanis, Huang and Simons [11].

Definition 2.2**.**

A random vector $\boldsymbol{X}:=(X_{1},\ldots,X_{p})\in\mathbb{R}^{p}$ is elliptically contoured, or simply elliptical if

[TABLE]

where $\boldsymbol{\mu}\in\mathbb{R}^{p}$, $\boldsymbol{A}\in\mathbb{R}^{p\times r}$, $\boldsymbol{U}$ is uniformly distributed on $\mathbb{S}^{r-1}$ for some $r>0$, and $R$ is a non-negative random variable independent of $\boldsymbol{U}$. Write $\boldsymbol{X}\sim\mathcal{E}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},F_{R})$, where $\boldsymbol{\Sigma}=\boldsymbol{A}\boldsymbol{A}^{T}$, and $F_{R}$ is the cumulative distribution function of $R$.

If $\boldsymbol{\mu}=\boldsymbol{0}$, $r=p$, and $\boldsymbol{A}$ is orthogonal, then $\boldsymbol{X}\sim\mathcal{E}_{p}(\boldsymbol{0},\boldsymbol{I}_{p},F_{R})$ is said to be spherically symmetric. For $\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{p})$, the random variable $R$ is chi-distributed with degree of freedom $p$. Below we list a few useful properties of elliptical distributions, see [10, 11, 19] for further development.

1. (1)

The random vector $\boldsymbol{X}\sim\mathcal{E}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},F_{R})$ has finite moments of order $k>0$ if and only if $\mathbb{E}R^{k}<\infty$. If the first two moments exist, then $\mathbb{E}\boldsymbol{X}=\boldsymbol{\mu}$ and $\operatorname{Var}\boldsymbol{X}=\frac{\mathbb{E}R^{2}}{r}\boldsymbol{\Sigma}$. 2. (2)

The distribution of $\boldsymbol{X}\sim\mathcal{E}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},F_{R})$ is absolutely continuous with respect to Lebesgue measure on $\mathbb{R}^{p}$ if and only if $r=p=\operatorname{rank}(\boldsymbol{\Sigma})$, and the distribution of $R$ is absolutely continuous with respect to Lebesgue measure on $\mathbb{R}_{+}$. 3. (3)

The marginal and conditional distributions of $\boldsymbol{X}\sim\mathcal{E}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},F_{R})$ are also elliptical. For the sake of simplicity, assume that $r=p=\operatorname{rank}(\boldsymbol{\Sigma})$. Let $\boldsymbol{X}:=(\boldsymbol{X}_{1},\boldsymbol{X_{2}})$, with $\boldsymbol{X}_{1}\in\mathbb{R}^{p_{1}}$ and $\boldsymbol{X}_{2}\in\mathbb{R}^{p_{2}}$ $(p_{1}+p_{2}=p)$. Let $\boldsymbol{\mu}:=\begin{pmatrix}\boldsymbol{\mu}_{1}\\ \boldsymbol{\mu}_{2}\end{pmatrix}$ and $\boldsymbol{\Sigma}:=\begin{pmatrix}\boldsymbol{\Sigma}_{11}&\boldsymbol{\Sigma}_{12}\\ \boldsymbol{\Sigma}_{21}&\boldsymbol{\Sigma}_{22}\end{pmatrix}$, with $\boldsymbol{\mu}_{1}\in\mathbb{R}^{p_{1}}$, $\boldsymbol{\mu}_{2}\in\mathbb{R}^{p_{2}}$, $\boldsymbol{\Sigma}_{11}\in\mathcal{M}_{p_{1}}(\mathbb{R})$, $\boldsymbol{\Sigma}_{12}=\boldsymbol{\Sigma}_{21}^{T}\in\mathcal{M}_{p_{1},p_{2}}(\mathbb{R})$, and $\boldsymbol{\Sigma}_{22}\in\mathcal{M}_{p_{2}}(\mathbb{R})$. Then $\boldsymbol{X}_{1}\sim\mathcal{E}_{p_{1}}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},F_{R})$, and $(\boldsymbol{X}_{1}|\boldsymbol{X}_{2}=\boldsymbol{x}_{2})\sim\mathcal{E}_{p_{1}}(\boldsymbol{\mu}_{1|2},\boldsymbol{\Sigma}_{1|2},F_{R_{1|2}})$ where

[TABLE]

and

[TABLE]

with $d_{\boldsymbol{\Sigma}_{22}}(\boldsymbol{x}_{2},\boldsymbol{\mu}_{2}):=\sqrt{(\boldsymbol{x}_{2}-\boldsymbol{\mu}_{2})^{T}\boldsymbol{\Sigma}^{-1}_{22}(\boldsymbol{x}_{2}-\boldsymbol{\mu}_{2})}$ the Mahalanobis distance between $\boldsymbol{x}_{2}$ and $\boldsymbol{\mu}_{2}$ in the metric associated with $\boldsymbol{\Sigma}_{22}$.

2.3. Main Result

Before stating the main result, we make a few assumptions on the link function, the covariate distribution, and model parameters. As seen in Section 2.1, the existence of the maximum likelihood estimate of binary response generalized linear models can be translated into data geometry. Thus, we need the assumptions on the link function in Theorem 2.1.

Assumption 2.3** (link function).**

For $\sigma:\mathbb{R}\rightarrow[0,1]$, both $\sigma(\cdot)$ and $1-\sigma(\cdot)$ are log-concave.

The phase transition for the existence of the maximum likelihood estimate is expected to occur not only with Gaussian covariates, but with a broad range of covariate distributions. Here we consider elliptical distributions. To exclude singularity, we assume that the covariate distribution is absolutely continuous with respect to Lebesgue measure on $\mathbb{R}^{p}$.

Assumption 2.4** (covariate distribution).**

The covariates $\boldsymbol{x}_{i}\sim\mathcal{E}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},F_{R})$ are of full rank. That is, $r=p=\operatorname{rank}(\boldsymbol{\Sigma})$ and $F_{R}$ is absolutely continuous on $\mathbb{R}_{+}$.

To get a meaningful result in diverging dimension, we consider a sequence of problems with the intercept $\beta_{0}$ fixed, and $\operatorname{Var}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta})\rightarrow\gamma_{0}^{2}$. Recall that for $\boldsymbol{x}_{i}\sim\mathcal{E}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},F_{R})$, we have $\operatorname{Var}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta})=\frac{\mathbb{E}R^{2}}{p}\boldsymbol{\beta}^{T}\boldsymbol{\Sigma}\boldsymbol{\beta}$. This leads to the following assumption on model parameters.

Assumption 2.5** (parameter scaling).**

As $p\rightarrow\infty$, $\mathbb{E}R^{2}/p\rightarrow\alpha_{0}^{2}$ and $|\boldsymbol{\Sigma}^{1/2}\boldsymbol{\beta}|\rightarrow\gamma_{0}/\alpha_{0}$.

In the remaining of the paper, we define

[TABLE]

where $U^{(p)}$ is distributed as any component of $\boldsymbol{x}_{i}\sim\mathcal{E}_{p}(\boldsymbol{0},\boldsymbol{I}_{p},F_{R})$, and $\mathbb{P}(V^{(p)}=1|U^{(p)})=1-\mathbb{P}(V^{(p)}=-1|U^{(p)})=\sigma(\beta_{0}+\gamma_{0}U^{(p)}/\alpha_{0})$. The parametric distribution $F_{\alpha_{0},\beta_{0},\gamma_{0}}$ is an analog of $F_{\beta_{0},\gamma_{0}}$ introduced in [12]. Here the superscript $(p)$ emphasizes that the distribution of $(Y^{(p)},X^{(p)})$ or $(V^{(p)},U^{(p)})$ may depend on $p$. The dependence in the Gaussian case is easily ignored since for $\boldsymbol{x}_{i}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{p})$, $U^{(p)}$ is distributed as standard normal independent of $p$. For the general elliptical covariates, we make the following technical assumption.

Assumption 2.6** (projection limit).**

The projection $U^{(p)}$ converges in distribution to $U$. That way, $(Y^{(p)},X^{(p)})$ converges in distribution to $(Y,X)$.

Now we give a sufficient condition for Assumption 2.6 to hold. Let $m_{p,k}$ be the $k^{th}$ moment of $U^{(p)}$. If for each $k\geq 1$, $m_{p,k}\rightarrow m_{k}$ as $p\rightarrow\infty$ and

[TABLE]

then $U^{(p)}$ converges in distribution to $U$ whose distribution is entirely characterized by the moments $(m_{k};\,k\geq 1)$. The condition (2.3) is referred to as the Carleman’s condition. See Lin [30, Theorem 1] for a list of equivalent conditions.

Define

[TABLE]

Denote $f_{X}$ as the density of $X$. Also define

[TABLE]

So $G_{p,+}(x)+G_{p,-}(x)$ is the cumulative distribution function of $X^{(p)}$, and $G_{p,\pm}(x)+\overline{G}_{p,\pm}(x)=\mathbb{E}p_{\pm}(X^{(p)})$. The main result is stated as follows.

Theorem 2.7**.**

Let $(Y^{(p)},X^{(p)})\sim F_{\alpha_{0},\beta_{0},\gamma_{0}}$, and $(Y,X)$ be the limit in distribution under Assumption 2.6. Let $Z\sim\mathcal{N}(0,1)$ be independent of $(Y,X)$. Define

[TABLE]

where $x_{+}:=\max\{x,0\}$. If Assumptions 2.3-2.5 are satisfied and $\sup_{p}\mathbb{E}[(X^{(p)})^{8}]$ $<\infty$, and

[TABLE]

then we have

[TABLE]

The proof is deferred to Section 4. The assumption $\sup_{p}\mathbb{E}[(X^{(p)})^{8}]<\infty$, which is purely technical, is used to prove the law of large numbers for triangular arrays. As suggested by simulations in Section 3, the phase transition exists even without this moment condition. The assumption (2.7) is used to prove that the probability the data points can be separated via a univariate model is small. A sufficient condition for (2.7) to hold is that $G_{p,\mp}+\overline{G}_{p,\pm}$ is bounded away from $1$. That is, there exists $\epsilon>0$ independent of $p$ such that

[TABLE]

So the term $\mathbb{E}[p_{\pm}(X^{(p)})(G_{\mp}(X^{(p)})+\overline{G}_{\pm}(X^{(p)}))^{n-1}]$ is exponentially small in $n$. To illustrate, we check the condition (2.8) for the logistic regression with Gaussian covariates. To simplify the discussion, we take $\beta_{0}=0$, and $\alpha_{0}=\gamma_{0}=1$. In this case, $p_{+}(x)=1-p_{-}(x)=e^{x}/(1+e^{x})$ and $X^{(p)}=X\sim\mathcal{N}(0,1)$. Fixing $x\geq 0$, we get

[TABLE]

Similarly, we can prove that $G_{p,-}(x)+\overline{G}_{p,+}(x)<1$ for $x<0$. It is also easily checked that all examples in Section 3 except the log-normal distribution satisfy the sufficient conditions (2.3)-(2.8).

3. Empirical Results

In this section we perform experiments to verify the phase transition of the maximum likelihood estimate by $(i)$ computing $h_{\text{MLE}}$ defined by (2.7); $(ii)$ checking whether the data is separated by the linear programming as in [12]:

[TABLE]

Note that the maximum likelihood estimate of the binary response generalized linear model exists if the linear programming (3.1) only has the trivial solution. We compare the theoretical phase transition curve with the empirical observations under several simulation designs described as below.

First, we argue that $\boldsymbol{A}=\boldsymbol{I}_{p}$ suffices to validate our theory. Consider $\tilde{\boldsymbol{X}}=\boldsymbol{A}^{-1}\boldsymbol{X}$ and $\tilde{\boldsymbol{\beta}}=\boldsymbol{A}^{\text{T}}\boldsymbol{\beta}$ so that $\tilde{\boldsymbol{X}}\tilde{\boldsymbol{\beta}}=\boldsymbol{X}\boldsymbol{\beta}$. The conditional distribution of the response $Y$ given $\tilde{\boldsymbol{X}}\tilde{\boldsymbol{\beta}}$ is the same as that given $\boldsymbol{X}\boldsymbol{\beta}$. If $(\tilde{\boldsymbol{X}}_{1},Y_{1}),\ldots,(\tilde{\boldsymbol{X}}_{n},Y_{n})$ are linearly separated by a hyperplane $\boldsymbol{\beta}^{*}$, then $(\boldsymbol{X}_{1},Y_{1}),\ldots,(\boldsymbol{X}_{n},Y_{n})$ are linearly separated by a hyperplane $(\boldsymbol{A}^{\text{T}})^{-1}\boldsymbol{\beta}^{*}$. The other way around also holds. Note that $\tilde{\boldsymbol{X}}\sim\mathcal{E}_{p}(\boldsymbol{A}^{-1}\boldsymbol{\mu},\boldsymbol{I}_{p},F_{R})$, so it suffices to consider $\boldsymbol{A}=\boldsymbol{I}_{p}$.

For the elliptical covariates, we set $\boldsymbol{\mu}=\mathbb{0}$, $\boldsymbol{A}=\boldsymbol{I}_{p}$ and consider different distributions for the non-negative variable $R$, including the chi distribution with degree of freedom $p$, the Gamma distributions, the Pareto distributions, the half-normal distribution and the log-normal distribution. When generating the binary response by (1.1), we choose the logit function, the cloglog function and the probit function for the link functino. We simply take $\beta_{0}=0$. See [12] for results with $\beta_{0}\neq 0$. To ensure Assumption 2.5, let

[TABLE]

where $\widetilde{\bf{W}}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{p})$.

We fix $n=1000$, $\alpha_{0}=1$, and vary $\gamma_{0}\in\{0.01,0.02,\ldots,10.00\}$ and $\kappa=p/n\in\{0.005,0.01,\ldots,0.6\}$. The parameter $\alpha_{0}$ is simply set as 1 since we observed that it does not affect the phase transition much. A large $\alpha_{0}$ might slightly shift the phase transition curve (the reds curve in Figure 1) to the right, and enlarge the uncertain band (the green bands). Once the data is generated, we solve the problem (3.1) by checking whether a non-trivial solution exists. We repeat the procedure for $100$ times, and get a heat map which indicates the proportion of times that the maximum likelihood estimate exists for each pair $(\gamma_{0},\kappa)$. See Figure 1 for the chi case with degree of freedom $p$. Results for other designs can be found in Appendix A.

3.1. Multivariate Gaussian covariates with different link functions

We consider the multivariate Gaussian covariates, which have been studied for the logit link in [12]. In our setup, $R$ is sampled from a chi distribution with degree of freedom $p$ and the link function is one of {logit, cloglog, probit}. Figure 1 displays the phase transition for the existence of the maximum likelihood estimate for different link functions. There are green bands in these figures, which indicates that the maximum likelihood estimate exists indefinitely when $(\gamma_{0},\kappa)$ falls in this band with the given sample size. This region is referred to as the uncertainty band. Observe that for the multivariate Gaussian covariates, as expected, the $h_{\text{MLE}}$ curves lie in the uncertainty bands for different link functions.

We also use the same setup to investigate the situation when $\boldsymbol{A}\neq\boldsymbol{I}_{p}$. In this case, we generate $\boldsymbol{A}$ such that each entry $A_{ij}$ is i.i.d sampled from the standard Gaussian distribution. To make sure it is full rank, we let $\boldsymbol{A}\leftarrow\frac{1}{10000}\|\boldsymbol{A}\|_{2}\cdot\boldsymbol{I}_{p}+\boldsymbol{A}$. The result is deferred to Figure 3 (Middle). We find that the result is very similar to that in Figure 1, which corroborates our argument that it suffices to use $\boldsymbol{A}=\boldsymbol{I}_{p}$ to validate our theory.

3.2. Gamma-distributed $R$

In [12], $U^{(p)}$ defined by (2.2) simplifies to the Gaussian distribution, which does not depend on $p$. This is key to their proof of the phase transition for the existence of the maximum likelihood estimate. However, when we go beyond the chi distribution for $R$, $U^{(p)}$ depends on $p$. We observe that Assumption 2.6 is satisfied for Gamma distributions, and the resulting theoretical phase transition curves agree with the simulations. More precisely, assume $R\sim\text{Gamma}(k,\theta)$ where $k$ is the shape parameter and $\theta$ is the scale parameter. The second moment condition (3.2) gives $\theta=\sqrt{\mathbb{E}R^{2}/(k^{2}+k)}$. When $k=0.5$, we get $\theta_{0.5}=\sqrt{4(p+1)/3}$ which corresponds to $\chi$ distribution with degree of freedom $2$ if $\theta_{0.5}$ is an integer; when $k=1$, it is the Exponential distribution with $\theta_{1}=\sqrt{(p+1)/2}$; when $k=2$, it is a Gamma distribution with $\theta_{2}=\sqrt{(p+1)/6}$. Figure 2 implies that $h_{\text{MLE}}$ defined by (2.7) converges quickly as $p$ increases. Table 1 indicates that all the theoretical phase transition curves align with the corresponding middle curves of the uncertainty bands.

3.3. The moment condition and the tail behavior of $R$

First we explore a case where the eighth moment of $X^{(p)}$ does not exist as required by Theorem 2.7. To this end, we sample $R$ from the Pareto distribution of type I. We specify the shape parameter $\alpha$, and set the scale parameter $x_{m}=\sqrt{(\alpha-2)/\alpha\cdot\mathbb{E}R^{2}}$. Recall that for the Pareto distribution of type I, the fourth moment exists when $\alpha>4$, and the third moment exists when $\alpha>3$. From Table 1, we see that the simulation results match the theoretical $h_{\text{MLE}}$ well. This suggests that the moment condition in Theorem 2.7 may be further relaxed.

Subsequently, we study how the tail behavior of the $R$ distribution influences the phase transition curve $h_{\text{MLE}}$. In the previous empirical studies, we consider the chi distributions with degree of freedom $p$ and Gamma distributions which have sub-exponential tails; the Pareto distributions have polynomial tails. We also investigate the half-normal distribution with a sub-Gaussian tail, and the log-normal distribution with another heavy tail. To ensure (3.2), we set the scale parameter $\sigma^{2}=\mathbb{E}R^{2}=p+1$ for the half-normal distribution, and $\sigma^{2}=0.5\cdot\log\mathbb{E}R^{2}=\log\sqrt{p+1}$ for the log-normal distribution. From Table 1, we observe that the theoretical $h_{\text{MLE}}$ successfully predicts the phase transition in the simulations of the half-normal distribution, but fails for the log-normal distribution. This is due to the fact that the log-normal distribution is not uniquely characterized by its moments, and the sufficient condition (2.3) does not hold. See Appendix A for more detailed results for the exploration.

4. Proof of Theorem 2.7

4.1. Roadmap to the proof of Theorem 2.7

Elliptical covariates Assume that the covariates $\boldsymbol{x}_{i}\sim\mathcal{E}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},F_{R})$ with $r=p=\operatorname{rank}(\boldsymbol{\Sigma})$. It is easily seen that $\boldsymbol{x}_{i}=\boldsymbol{\mu}+\boldsymbol{\Sigma}^{1/2}\boldsymbol{z}_{i}$ with $\boldsymbol{z}_{i}\sim\mathcal{E}_{p}(\boldsymbol{0},\boldsymbol{I}_{p},F_{R})$. Recall from Section 2.1 that for GLMs satisfying Assumption 2.3, the MLE does not exist if and only if there is $(b_{0},\boldsymbol{b})\neq\boldsymbol{0}$ such that $y_{i}(b_{0}+\boldsymbol{x}_{i}^{T}\boldsymbol{b})\geq 0$ for all $i$. This is equivalent to the existence of $(\widetilde{b}_{0},\widetilde{\boldsymbol{b}})\neq\boldsymbol{0}$ such that $y_{i}(\widetilde{b}_{0}+\boldsymbol{z}_{i}^{T}\widetilde{\boldsymbol{b}})\geq 0$ for all $i$. Without loss of generality, we assume $\boldsymbol{x}_{i}\sim\mathcal{E}_{p}(\boldsymbol{0},\boldsymbol{I}_{p},F_{R})$ in the sequel.

We are in a situation where the covariates $\boldsymbol{x}_{i}:=(x_{i1},\ldots,x_{ip})$ is spherically symmetric and $\operatorname{Var}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta})=|\boldsymbol{\beta}|^{2}\mathbb{E}R^{2}/p\rightarrow\gamma_{0}^{2}$. By rotational invariance, we assume that all the signal is in the first coordinate. That is, $\mathbb{P}(y_{i}=1|\boldsymbol{x}_{i})=\sigma\left(\beta_{0}+\frac{\gamma_{0}}{\alpha_{0}}x_{i1}\right)$. The results in Section 2.2 show that

[TABLE]

where $(Y^{(p)},X^{(p)})\sim F_{\alpha_{0},\beta_{0},\gamma_{0}}$, and

[TABLE]

with $F_{R_{-1}}(r)=\frac{\int_{|x|}^{\sqrt{r^{2}+x^{2}}}(s^{2}-x^{2})^{(p-3)/2}s^{-p+2}dF_{R}(s)}{\int_{|x|}^{\infty}(s^{2}-x^{2})^{(p-3)/2}s^{-p+2}dF_{R}(s)}$.

Now we want to express $\mathbb{P}(\mbox{no MLE})$ via conic geometry. For a fixed space $\mathcal{W}\in\mathbb{R}^{n}$, let

[TABLE]

be the convex cone generated by $\mathcal{W}$. The following proposition is read from [12, Propositions 1 $\&$ 2], and we include the proof in Section 4.2 for completeness.

Proposition 4.1**.**

Let the $n$-dimensional vectors $(\boldsymbol{Y}^{(p)},\boldsymbol{X}^{(p)},\boldsymbol{X}_{2},\ldots,\boldsymbol{X}_{p})$ be $n$ i.i.d. copies of $(Y^{(p)},X^{(p)},X_{2},\ldots,X_{p})$ distributed as in (4.1). Let

[TABLE]

Let $\{\mbox{No MLE Single}\}$ be the event that the data points can be completely or quasi-completely separated by the intercept and the first coordinate only, i.e. $\mathcal{W}\cap\mathbb{R}_{+}^{n}\neq\{\boldsymbol{0}\}$. Then

[TABLE]

By Proposition 4.1, the existence of the MLE boils down to whether $\mathcal{L}$ intersects $\mathcal{C}(\mathcal{W})$ in a non-trivial way. It remains to prove the following: $(1)$ The probability $\mathbb{P}(\mbox{No MLE Single})$ is relatively small. $(2)$ The probability $\mathbb{P}(\mathcal{L}\cap\mathcal{C}(\mathcal{W})\neq\{\boldsymbol{0}\})$ exhibits a phase transition through the ratio $\kappa:=p/n$, and $h_{\text{MLE}}(\alpha_{0},\beta_{0},\gamma_{0})$ defined by (2.7).

Separation of data in a univariate model We aim to prove that the probability $\mathbb{P}(\mbox{No MLE Single})$ is small. In [12], a sketch of proof is given for the logistic regression with Gaussian covariates. In Section 4.3, we give a rigorous proof of this result in the setting of Theorem 2.7. The main difficulty comes from the fact that though the probability the data can be separated via any fixed $t_{0}\in\mathbb{R}$ is exponentially small, there are uncountably many such $t_{0}$ and the union bound does not give a good estimate.

Proposition 4.2**.**

Under the assumptions in Theorem 2.7, the event $\{\mbox{No MLE}$ $\mbox{Single}\}$ occurs with small probability. That is, $\mathbb{P}(\mbox{No MLE Single})=o(1)$.

Convex geometry and phase transition We want to prove the phase transition of $\mathbb{P}(\mathcal{L}\cap\mathcal{C}(\mathcal{W}))$ through the interplay between $\kappa$ and $h_{\text{MLE}}(\alpha_{0},\beta_{0},\gamma_{0})$. The key is to understand when a random subspace $\mathcal{L}$ with uniform orientation intersects $\mathcal{C}(\mathcal{W})$ in a non-trivial way.

For any fixed subspace $\mathcal{W}\in\mathbb{R}^{n}$, the approximate kinematic formula [2, Theorem I] shows that for any $\varepsilon\in(0,1)$, there exists $a_{\varepsilon}>0$ such that

[TABLE]

Here $\delta(\mathcal{C})$ is the statistical dimension of the convex cone $\mathcal{C}$ defined by $\delta(\mathcal{C}):=n-\mathbb{E}|\boldsymbol{Z}-\Pi_{\mathcal{C}}(\boldsymbol{Z})|^{2}$, where $\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{n})$ and $\Pi_{\mathcal{C}}$ is the projection onto $\mathcal{C}$. The following identity is given in [12, Lemma 3]:

[TABLE]

Theorem 2.7 can be derived from the formulas (4.1)-(4.4) and the following theorem.

Theorem 4.3**.**

Let $(\boldsymbol{Y}^{(p)},\boldsymbol{X}^{(p)})$ be $n$ i.i.d. samples from $F_{\alpha_{0},\beta_{0},\gamma_{0}}$ satisfying Assumption 2.6, and $\sup_{p}\mathbb{E}[(X^{(p)})^{8}]<\infty$. Let

[TABLE]

Then $Q_{p,n}$ converges in probability to $\min_{\lambda_{0},\lambda_{1}\in\mathbb{R}}\mathbb{E}(\lambda_{0}Y+\lambda_{1}X-Z)_{+}^{2}$ as $n,p\rightarrow\infty$.

In [12], the authors proved Theorem 4.3 in the setting of the logistic regression by a bare-hands argument. One can adapt their argument to prove Theorem 4.3, with possibly more complications. However, the statement of Theorem 4.3 suggests it be a form of stochastic approximation. Here we show how this result follows systematically from stochastic approximation, which is of independent interest.

Stochastic approximation We sketch a proof of Theorem 4.3 via a stochastic approximation. In the stochastic approximation literature [26, 35], people seek to approximate the optimization problem

[TABLE]

where $S\subset\mathbb{R}^{k}$ for some $k>0$ and $\boldsymbol{\xi}$ is a generic random vector, by a sequence of stochastic optimization problems $\min_{\boldsymbol{\lambda}\in S}\widehat{G}_{n}(\boldsymbol{\lambda};\boldsymbol{\xi}_{1},\ldots,\boldsymbol{\xi}_{n})$, $\widehat{G}_{n}(\boldsymbol{\lambda};\boldsymbol{\xi}_{1},\ldots,\boldsymbol{\xi}_{n}):=\frac{1}{n}\sum_{i=1}^{n}g(\boldsymbol{\lambda},\boldsymbol{\xi}_{i})$ where $(\boldsymbol{\xi}_{i};\,i\geq 1)$ are i.i.d. copies of $\boldsymbol{\xi}$.

A deep connection between stochastic approximation and convergence of random closed sets was established by Attouch and Wets [5] via the epi-convergence of functions. A sequence of lower semi-continuous functions $f_{n}:\mathbb{R}^{k}\rightarrow(\infty,\infty]$ is said to epi-converges to $f$ if for each $x\in\mathbb{R}^{k}$,

• •

$\lim\inf f_{n}(x_{n})\rightarrow f(x)$ if $x_{n}\rightarrow x$,

• •

$\lim f_{n}(x_{n})\rightarrow f(x)$ for at least one sequence $x_{n}\rightarrow x$.

See [3, 4, 17, 23, 27, 37] for further development on epi-convergence.

Here we consider a sequence of stochastic optimization problems with triangular arrays

[TABLE]

where $(\boldsymbol{\xi}_{i,n};\,1\leq i\leq n)$ are i.i.d. copies of $\boldsymbol{\xi}_{n}$, and $\boldsymbol{\xi}_{n}$ converges in distribution to $\boldsymbol{\xi}$. Let $v$, $\widehat{v}_{n}$, and $\operatorname{argmin}G$, $\operatorname{argmin}\widehat{G}_{n}$ be optimal values, and optimal solutions to the problems (4.5)-(4.6). Note that $\operatorname{argmin}\,G$ and $\operatorname{argmin}\,\widehat{G}_{n}$ are set-valued. The following result gives asymptotic inference of $\widehat{v}_{n}$ as $n\rightarrow\infty$. The proof will be given in Section 4.4.

Lemma 4.4**.**

Assume that $g(\cdot,\cdot)$ is measurable and bounded from below, and $g(\cdot,\boldsymbol{\xi})$ is convex. Assume that $\sup_{n}\mathbb{E}g^{4}(\boldsymbol{\lambda},\boldsymbol{\xi}_{n})<\infty$ and $\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi}_{n})\rightarrow\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi})$ for all $\boldsymbol{\lambda}$. Further assume that $\operatorname{argmin}\,\widehat{G}_{n}$, $n\geq 1$ are non-empty and bounded in probability. Then $\operatorname{argmin}\,G\neq\emptyset$, and

[TABLE]

To prove Theorem 4.3, we need to show that the set of minimizers $\operatorname{argmin}\,\widehat{G}_{n}$ is non-empty and bounded in probability. The argument is similar in spirit to [12], and we give the proof for ease of reference. In Section 4.4, we prove that under the assumptions in Theorem 2.7:

• •

The problem (4.5) has a unique minimizer $\boldsymbol{\lambda}_{0}$.

• •

For any minimizer $\widehat{\boldsymbol{\lambda}}_{n}\in\operatorname{argmin}\,\widehat{G}_{n}$, $|\widehat{\boldsymbol{\lambda}}_{n}-\boldsymbol{\lambda}_{0}|=\mathcal{O}_{P}(1)$. Here we use the Chebyshev inequality instead of the concentration inequality for sub-exponential variables in [12].

Proof of Theorem 2.7 The proof goes along the same line as [12], with two modifications. Assume that $\kappa>h_{\text{MLE}}(\alpha_{0},\beta_{0},\gamma_{0})$.

• •

Given $(Y^{(p)},X^{(p)})$, the random vector $(X_{2},\ldots,X_{p})$ is also elliptical whose distribution is given by (4.1). By the geometric characterization (4.1), we get

[TABLE]

for some $\varepsilon_{n}\rightarrow 0$.

• •

The random variable $Q_{p,n}$ is uniformly integrable, since $Q_{p,n}\leq|\boldsymbol{Z}_{+}|^{2}/n$ which is sub-exponential. This implies that $\mathbb{E}(Q_{p,n}|\boldsymbol{X},\boldsymbol{Y})$ converges in probability to $h_{\text{MLE}}(\alpha_{0},\beta_{0},\gamma_{0})$.

Thus, $\mathbb{P}(\mathcal{L}\cap\mathcal{C}(\mathcal{W})\neq\{\boldsymbol{0}\})\rightarrow 1$ and $\mathbb{P}(\mbox{MLE exists})\rightarrow 0$. Similarly, we can prove if $\kappa<h_{\text{MLE}}(\alpha_{0},\beta_{0},\gamma_{0})$, then $\mathbb{P}(\mbox{MLE exists})\rightarrow 1$.

4.2. Proof of Proposition 4.1

We aim to prove that

[TABLE]

from which the result follows. If $\{\mbox{No MLE Single}\}$ occurs, there is no MLE. Assume that $\{\mbox{No MLE Single}\}$ does not occur. If

[TABLE]

then there is no MLE if and only if there is a non-zero vector $(b_{0},\ldots,b_{p})$ such that $b_{0}\boldsymbol{Y}^{(p)}+b_{1}\boldsymbol{X}^{(p)}+\cdots+b_{p}\boldsymbol{X}_{p}=\boldsymbol{u}$, $\boldsymbol{u}\geq 0,\boldsymbol{u}\neq 0$. By assumption, $b_{0}\boldsymbol{Y}^{(p)}+b_{1}\boldsymbol{X}^{(p)}\neq u$ so $b_{2}\boldsymbol{X}_{2}+\cdots b_{p}\boldsymbol{X}_{p}$ is a non-zero element of $\mathcal{C}(\mathcal{W})$. This leads to (4.7). Note that there is no MLE if and only if there is a non-zero vector $(b_{0},\ldots,b_{p})$ such that $b_{0}\boldsymbol{Y}^{(p)}+b_{1}\boldsymbol{X}^{(p)}+\cdots+b_{p}\boldsymbol{X}_{p}\geq 0$. The identity in law (9) implies that the equality occurs with probability [math], which proves (4.8).

4.3. Proof of Proposition 4.2

Let $(X_{1},\ldots,X_{n})$ be i.i.d. samples with density $f_{X^{(p)}}$. It is well known that the distribution of the order statistics $(X_{(1)},\ldots,X_{(n)})$ is given by $n!\,\prod_{i=1}^{n}f_{X^{(p)}}(x_{i})$ for $x_{1}<\cdots<x_{n}$. Note that there exists $t\in\mathbb{R}$ separating $X_{(1)}<\cdots<X_{(n)}$ if and only if for some $k\in\{0,\ldots,n\}$, the responses corresponding to $X_{(1)},\ldots,X_{(k)}$ is of the same sign, and those corresponding to $X_{(k+1)},\ldots,X_{(n)}$ is of the opposite sign. Consequently,

[TABLE]

Moreover, $\int_{x_{1}<\cdots<x_{n}}n!\,\prod_{i=1}^{n}f_{X^{(p)}}(x_{i})\prod_{i=1}^{n}p_{\pm}(x_{i})dx_{1}\ldots dx_{n}$ is equal to

[TABLE]

For $1\leq k\leq n-1$,

[TABLE]

Combining the above identities yields

[TABLE]

Finally, the condition (2.7) together with (4.9) lead to the desired result.

4.4. Proof of Theorem 4.3

We start with the proof of Lemma 4.4.

Proof of Lemma 4.4.

By law of large numbers of triangular arrays, the condition $\sup_{i}\mathbb{E}g^{4}(\boldsymbol{\lambda},\boldsymbol{\xi}_{i,n})<\infty$ implies that $\frac{1}{n}\sum_{i=1}^{n}g(\boldsymbol{\lambda},\boldsymbol{\xi}_{i,n})-\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi}_{n})\rightarrow 0$ a.s. It follows from [3, Theorem 2.3] that

[TABLE]

It is well known that if a sequence of convex functions converge pointwise, then they converge uniformly on compact sets. Since $\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi}_{n})\rightarrow\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi})$ for all $\boldsymbol{\lambda}$ and $g(\cdot,\boldsymbol{\xi})$ is convex, the convergence is uniform on compact sets. This implies the epi-convergence. Therefore,

[TABLE]

Combining with [17, Proposition 3.3] yields the desired result. ∎

Now we are ready to prove Theorem 4.3. We specialize to $\boldsymbol{\lambda}=(\lambda_{0},\lambda_{1})$, $\boldsymbol{\xi}_{n}=(Y^{(p)},X^{(p)},Z)$ with $(Y^{(p)},X^{(p)})\sim F_{\alpha_{0},\beta_{0},\gamma_{0}}$, $\boldsymbol{\xi}=(Y,X,Z)$, and $Z\sim\mathcal{N}(0,1)$ independent of $\{(Y^{(p)},X^{(p)}),(Y,X)\}$, and

[TABLE]

It is clear that the function $g$ defined by (4.10) is measurable and non-negative, and $g(\cdot,\boldsymbol{\xi})$ is convex. It follows from $\sup_{p}\mathbb{E}[(X^{(p)})^{8}]<\infty$ that $\sup_{n}g^{4}(\boldsymbol{\lambda},\boldsymbol{\xi}_{n})<\infty$. By Assumption 2.5, $\mathbb{E}[(X^{(p)})^{2}]$ converges to $\alpha_{0}^{2}$, and by Assumption 2.6, $(Y^{(p)},X^{(p)})$ converges in distribution to $(Y,X)$. Now by [8, Lemma 8.3], we get $\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi}_{n})\rightarrow\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi})$ for all $\boldsymbol{\lambda}$.

Let $\boldsymbol{\lambda}_{\min}$ be a minimum of $\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi})$. By convexity of $\boldsymbol{\lambda}\rightarrow\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi})$, there exists $r>0$ such that $\mathbb{E}g(\boldsymbol{\lambda}_{\min},\boldsymbol{\xi})<\min_{r\leq|\boldsymbol{\lambda}|\leq r+1}\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi})$, and $\min_{\boldsymbol{\lambda}}\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi})=\min_{|\boldsymbol{\lambda}|\leq r+1}\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi})$. Note that $\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi}_{n})$ converges uniformly to $\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi})$ on $\{\boldsymbol{\lambda}:r\leq|\boldsymbol{\lambda}|\leq r+1\}$. So for $p$ large enough, $\mathbb{E}g(\boldsymbol{\lambda}_{\min},\boldsymbol{\xi}_{n})<\min_{r\leq|\boldsymbol{\lambda}|\leq r+1}\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi}_{n})$ and $\min_{\boldsymbol{\lambda}}\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi}_{n})=\min_{r\leq|\boldsymbol{\lambda}|\leq r+1}\mathbb{E}g(\boldsymbol{\lambda},\boldsymbol{\xi}_{n})$. Now by [24, Theorem 2.1], we have as $p\rightarrow\infty$,

[TABLE]

By Lemma 4.4, it suffices to prove that the set of minimizers $\operatorname{argmin}\widehat{G}_{n}$ is non-empty and bounded in probability. We aim to show that under the assumptions in Theorem 2.7, the problem $(13)$ has a unique minimizer $\boldsymbol{\lambda}_{0}$, and for any minimizer $\widehat{\boldsymbol{\lambda}}_{n}\in\operatorname{argmin}\widehat{G}_{n}$, $|\widehat{\boldsymbol{\lambda}}_{n}-\boldsymbol{\lambda}_{0}|=\mathcal{O}_{P}(1)$.

We prove these statements in the next two lemmas. It is easily seen that $G(\boldsymbol{\lambda})$ and $\widehat{G}_{n}(\boldsymbol{\lambda})$ are convex. The following lemma shows that the function $G$ is strongly convex, which was stated in [12] without proof. Here we give a complete proof.

Lemma 4.5**.**

Under the assumptions in Theorem 2.7, the function $\boldsymbol{\lambda}\rightarrow G(\boldsymbol{\lambda})$ with $g(\cdot)$ defined by (4.10) is strongly convex. That is, there exists $\alpha_{1}>\alpha_{0}>0$ such that

[TABLE]

Proof.

Elementary analysis shows that

[TABLE]

where $\Phi(\cdot)$ is the CDF of standard normal, $U$ is defined in Assumption 2.6, and $\mathbb{P}(V=1|U)=1-\mathbb{P}(V=-1|U)=p_{+}(U)$. The r.h.s. of (4.11) is clear. By Cauchy-Schwarz inequality, $\det\nabla^{2}G(\boldsymbol{\lambda})\geq 0$. If $\det\nabla^{2}G(\boldsymbol{\lambda})=0$, then $U$ is constant almost surely which violates the non-degeneracy of $U$. Thus, $G(\boldsymbol{\lambda})$ is strictly convex.

Note that

[TABLE]

and the decomposition holds for other terms. So

[TABLE]

Without loss of generality, consider $\lambda_{1},\lambda_{2}>0$. For $\lambda_{0},\lambda_{1}$ sufficiently large,

• •

if $\lambda_{0}/\lambda_{1}$ is large, then $\mathbb{E}[\Phi(\lambda_{0}+\lambda_{1}U)p_{+}(U)]$ can be approximated by $\mathbb{E}p_{+}(U)$.

• •

if $\lambda_{0}/\lambda_{1}$ is small, then $\mathbb{E}[\Phi(\lambda_{0}+\lambda_{1}U)p_{+}(U)]$ can be approximated by $\mathbb{E}[1(U>-\lambda_{0}/\lambda_{1}+\eta)p_{+}(U)]$, where $\eta$ is a fixed small value.

The approximation also holds for other terms. By the strict convexity and the approximation, we can show that there exist $N_{+}>0$ such that for $\lambda_{0},\lambda_{1}>N_{+}$,

[TABLE]

for some $\epsilon_{+}>0$. Similarly, there exist $N_{-}>0$ such that for $\lambda_{0},\lambda_{1}>N_{-}$, we get a bound $\epsilon_{-}\boldsymbol{I}_{2}$ for the second term on the r.h.s. of (4.4). Thus, for $\lambda_{1},\lambda_{2}>\max(N_{+},N_{-})$, $\nabla^{2}G(\boldsymbol{\lambda})\succeq\min(\epsilon_{+},\epsilon_{-})\boldsymbol{I}_{2}$. By continuity of $\det\nabla^{2}G(\boldsymbol{\lambda})$, we get $\nabla^{2}G(\boldsymbol{\lambda})\succeq\epsilon^{\prime}\boldsymbol{I}_{2}$ for $\lambda_{0},\lambda_{1}<\max(N_{+},N_{-})$. It suffices to take $\alpha_{0}=\min(\epsilon_{+},\epsilon_{-},\epsilon^{\prime})$ to conclude. ∎

Finally, we prove that the set of minimizers $\operatorname{argmin}\,\widehat{G}_{n}$ is bounded in probability. The argument can be used to show that $|\widehat{\boldsymbol{\lambda}}_{n}-\boldsymbol{\lambda}_{0}|=\mathcal{O}_{P}(n^{-1/4})$, where $\widehat{\boldsymbol{\lambda}}_{n}$ is any minimizer of $\widehat{G}_{n}$. The proof is adapted from [12, Lemma 4], which we include for completeness.

Lemma 4.6**.**

Under the assumptions in Theorem 2.7, we have $|\widehat{\boldsymbol{\lambda}}_{n}-\boldsymbol{\lambda}_{0}|=\mathcal{O}_{P}(1)$, where $\widehat{\boldsymbol{\lambda}}_{n}\in\operatorname{argmin}\,\widehat{G}_{n}$.

Proof.

For any $\boldsymbol{\lambda}\in\mathbb{R}^{2}$, the strong convexity (4.11) gives that

[TABLE]

Fix $x\geq 1$. For any $\boldsymbol{\lambda}$ on the circle $C(x):=\{\boldsymbol{\lambda}:|\boldsymbol{\lambda}-\boldsymbol{\lambda}_{0}|=x\}$, we have

[TABLE]

Fix $z=G(\boldsymbol{\lambda}_{0})+y$, and consider the event

[TABLE]

By convexity of $\widehat{G}_{n}$, when $E$ occurs, $\widehat{\boldsymbol{\lambda}}_{n}$ must lie in the circle. Hence, $|\widehat{\boldsymbol{\lambda}}_{n}-\boldsymbol{\lambda}_{0}|\leq x$.

Next we prove that the event $E$ occurs with high probability. Fix $d$ equi-spaced point $\{\boldsymbol{\lambda}_{i}\}_{i=1}^{d}$ on the set $C(x)$. Take any point $\boldsymbol{\lambda}$ on the circle, and let $\boldsymbol{\lambda}_{i}$ be its closest point. So $|\boldsymbol{\lambda}-\boldsymbol{\lambda}_{i}|\leq\pi x/d$. By convexity of $\widehat{G}_{n}$,

[TABLE]

Define the event

[TABLE]

By Chebyshev inequality and union bound, we get

[TABLE]

where $\sigma^{2}:=\sup_{n}\operatorname{Var}[g(\boldsymbol{\lambda},\boldsymbol{\xi}_{n})]<\infty$. As $|\nabla^{2}G|$ is bounded and $\nabla G(\boldsymbol{\lambda}_{0})=0$,

[TABLE]

Now for $n$ sufficiently large, on $B$ we have $|\nabla\widehat{G}_{n}(\boldsymbol{\lambda}_{i})||\boldsymbol{\lambda}-\boldsymbol{\lambda}_{i}|\leq\pi(1+\alpha_{1})x^{2}/d$. Choose $d>6\pi(1+\alpha_{1})/\alpha_{0}$ so that $|\nabla\widehat{G}_{n}(\boldsymbol{\lambda}_{i})||\boldsymbol{\lambda}-\boldsymbol{\lambda}_{i}|\leq y$. Then it follows from (4.14) that on $B$,

[TABLE]

Observe that

[TABLE]

since $G(\boldsymbol{\lambda}_{i})\geq G(\boldsymbol{\lambda}_{0})+3y$ by (4.13). Consequently,

[TABLE]

where $\mathbb{P}(B^{c})$ is controlled by (4.15), and the last two terms can also be bounded by Chebyshev inequality. ∎

5. Conclusion

In this paper, we proved a phase transition for the existence of the maximum likelihood estimate in high-dimensional generalized linear models with elliptical covariates. We derived an explicit formula for the phase transition boundary, depending on the regression coefficients and the scaling parameter of the covariate distribution. Our result extends a previous one in [12], and elucidates a rich structure in the phase transition phenomenon. We believe that the phase transition also holds for multinomial response models such as the Poisson regression and the log-linear regression. See [15, 20] for further discussions. We hope that this work will trigger further research towards a theory of hypothesis testing for generalized linear models with non-Gaussian covariates.

Acknowledgements

Wenpin Tang gratefully acknowledges financial support through a startup grant at Columbia University.

Appendix A Empirical Results

To study whether the MLE exists for a given pair of $(\gamma_{0},\kappa)$ with some distribution for $R$, we generate the data by the mechanism described in Section 3. We fix the sample size at $n=1000$, and vary $\gamma_{0}\in\{0.01,0.02,\ldots,10.00\}$ and $\kappa=p/n\in\{0.005,0.01,\ldots,0.6\}$. We generate $100$ such datasets, and for each dataset we solve the linear programming (3.1) by checking whether there is a nonzero solution using the package CVXOPT111https://cvxopt.org/documentation/index.html in python.

On the other hand, the optimization problem (2.7) does not have a closed-form solution. Here we solve numerically this convex optimization. Using the same data generation mechanism but with a sample size $n=4000$, we compute $h_{\text{MLE}}$ by CVXOPT as well. We repeat the procedure for $100$ times and take the average of these replicates as the reported $h_{\text{MLE}}$. Here we take $\gamma_{0}\in\{0.5,1.0,\ldots,10.00\}$ and $\kappa=p/n\in\{0.02,0.04,\ldots,0.6\}$.