High-dimensional single-index models:
link estimation and marginal inference

We review key technical aspects of statistical inference for each coordinate $\beta_{j}$ in the proportionally high-dimensional regime (2). We maintain $\bm{\beta}^{\top}\bm{X}_{i}$ of constant order by considering the setting $\bm{\beta}^{\top}\bm{\Sigma}\bm{\beta}=\Theta(1)$. We define $\bm{\Theta}=\bm{\Sigma}^{-1}$ as the precision matrix for the distribution of $\bm{X}_{i}$ and set $\tau_{j}^{2}:=\Theta_{jj}^{-1}>0$. An unbiased estimator of $\tau_{j}$ can be constructed using nodewise regression (cf. Section 5.1 in Zhao et al. (2022)). For simplicity, we assume $\tau_{j}$ is known, following prior studies.

In the high-dimensional regime (2), statistical inference must address two components: the asymptotic distribution and the inferential parameters of an estimator. We review the asymptotic distribution of the MLE $\hat{\bm{\beta}}^{\mathrm{m}}$ for logistic regression. According to Zhao et al. (2022), for all $j\in\{1,\ldots,p\}$ such that $\sqrt{p}\tau_{j}\beta_{j}=O(1)$ as $n\to\infty$, the estimator achieves the following asymptotic normality:

Here, we define $\mu_{*}\in\mathbb{R}$ and $\sigma_{*}>0$ as the asymptotic bias and variance, respectively, ensuring the convergence (3). It is crucial to note that both the estimator $\hat{\beta}_{j}^{\mathrm{m}}$ and the target $\beta_{j}$ scale as $O_{p}(1/\sqrt{p})$ here.

To perform statistical inference based on (3), it is necessary to estimate the inferential parameters $\mu_{*}$ and $\sigma_{*}$. Several studies including El Karoui et al. (2013); Thrampoulidis et al. (2018); Sur and Candès (2019); Loureiro et al. (2021) theoretically characterize these parameters as solutions to a system of nonlinear equations that depend on the data-generating process and the loss function. Additionally, various approaches have been developed to practically solve the equations by determining their hyperparameter $\bm{\beta}^{\top}\bm{\Sigma}\bm{\beta}$ under different conditions. Specifically, Sur and Candès (2019) introduces ProbeFrontier for estimating $\bm{\beta}^{\top}\bm{\Sigma}\bm{\beta}$ based on the asymptotic existence/non-existence boundary of the maximum likelihood estimator (MLE) in logistic regression. SLOE, proposed by Yadlowsky et al. (2021), enhances this estimation using a leave-one-out technique. Moreover, Sawaya et al. (2023) takes a different approach to estimate $\bm{\beta}^{\top}\bm{\Sigma}\bm{\beta}$ for generalized linear models.

For single-index models, Bellec (2022) introduces observable adjustments that estimate the inferential parameters directly under the identification condition $\bm{\beta}^{\top}\bm{\Sigma}\bm{\beta}=1$ irrespective of link misspecification, bypassing the system of equations. In our study, we develop an estimator for the single-index model satisfying the asymptotic normality (3), with corresponding estimators for the inferential parameters using observable adjustments.

Research into the asymptotic behavior of statistical models in high-dimensional settings, where both $n$ and $p$ diverge proportionally, has gained momentum in recent years. Notable areas of exploration include (regularized) linear regression models (Donoho et al., 2009; Bayati and Montanari, 2011b; Krzakala et al., 2012; Bayati et al., 2013; Thrampoulidis et al., 2018; Mousavi et al., 2018; Takahashi and Kabashima, 2018; Miolane and Montanari, 2021; Guo and Cheng, 2022; Hastie et al., 2022; Li and Sur, 2023), robust estimation (El Karoui et al., 2013; Donoho and Montanari, 2016), generalized linear models (Rangan, 2011; Sur et al., 2019; Sur and Candès, 2019; Salehi et al., 2019; Barbier et al., 2019; Zhao et al., 2022; Tan and Bellec, 2023; Sawaya et al., 2023), low-rank matrix estimation (Deshpande et al., 2017; Macris et al., 2020; Montanari and Venkataramanan, 2021), and various other models (Montanari et al., 2019; Loureiro et al., 2021; Yang and Hu, 2021; Mei and Montanari, 2022). These investigations focus primarily on the convergence limits of estimation and prediction errors. Theoretical analyses have shown that classical statistical estimation often fails to accurately estimate standard errors and may lack key desirable properties such as asymptotic unbiasedness and asymptotic normality.

In such analyses, the following theoretical tools have been employed: (i) the replica method (Mézard et al., 1987; Charbonneau et al., 2023), (ii) approximate message passing algorithms (Donoho et al., 2009; Bolthausen, 2014; Bayati and Montanari, 2011a; Feng et al., 2022), (iii) the leave-one-out technique (El Karoui et al., 2013; El Karoui, 2018), (iv) the convex Gaussian min-max theorem (Thrampoulidis et al., 2018), (v) second-order Poincaré inequalities (Chatterjee, 2009; Lei et al., 2018), and (vi) second-order Stein’s formulae (Bellec and Zhang, 2021; Bellec and Shen, 2022). Although these tools were initially proposed for analyzing linear models with Gaussian design, they have been extensively adapted to a diverse range of models. In this study, we apply observable adjustments based on second-order Stein’s formulae (Bellec, 2022) to directly estimate the asymptotic bias and variance of coefficient estimators. Furthermore, we provide a comprehensive proof of marginal asymptotic normality, extending the work of Zhao et al. (2022) to a wider array of estimators.

Define $[z]=\{1,\ldots,z\}$ for $z\in\mathbb{N}$. For a vector $\bm{b}=(b_{1},\dots,b_{p})\in\mathbb{R}^{p}$, we write $\|\bm{b}\|:=(\sum_{j=1}^{p}b_{j}^{2})^{1/2}$. For a collection of indices $\mathcal{S}\subset[p]$, we define a sub-vector $\bm{b}_{\mathcal{S}}:=(b_{j})_{j\in\mathcal{S}}$ as a slice of $\bm{\beta}$. For a matrix $\bm{A}\in\mathbb{R}^{p\times p}$, we define its minimum and maximum eigenvalues by $\lambda_{\rm min}(\bm{A})$ and $\lambda_{\rm max}(\bm{A})$, respectively. For a function $F:\mathbb{R}\to\mathbb{R}$, we say $F^{\prime}$ the derivative of $F$ and $F^{(m)}$ the $m$th-order derivative. For a function $f:\mathbb{R}\to\mathbb{R}$ and a vector $\bm{b}\in\mathbb{R}^{p}$, $f(\bm{b})=(f(b_{1}),f(b_{2}),\dots,f(b_{p}))^{\top}\in\mathbb{R}^{p}$ denotes a vector by elementwise operations.

We organize the remainder of the paper as follows: Section 2 presents our estimation procedure. Section 3 describes the asymptotic properties of the proposed estimator and develops a statistical inference method. Section 4 provides several experiments to validate our estimation theory. Section 5 outlines the proofs of our theoretical results. Section 6 discusses alternative designs for estimators. Finally, Section 7 concludes with a discussion of our findings. The Appendix contains additional discussions and the complete proofs.

In this step, we use the first subset $(\bm{X}^{(1)},\bm{y}^{(1)})$. We define the pilot estimator as the ridge estimator, $\tilde{\bm{\beta}}=(({\bm{X}^{(1)}})^{\top}\bm{X}^{(1)}+n_{1}\lambda\bm{I}_{p})^{-1}({\bm{X}^{(1)}})^{\top}\bm{y}^{(1)}$ where $\lambda>0$ is the regularization parameter. Further, we consider inferential parameters $(\mu_{1},\sigma_{1})$ of $\tilde{\bm{\beta}}$, which satisfy $\sqrt{p}\tau_{j}(\tilde{\beta}_{j}-\mu_{1}\beta_{j})/\sigma_{1}\overset{\mathrm{d}}{\to}\mathcal{N}(0,1)$ for $j\in[p]$ such that $\sqrt{p}\tau_{j}\beta_{j}=O(1)$. Using these parameters, we develop an estimator $W_{i}$ for the index $\bm{\beta}^{\top}{\bm{X}_{i}^{(1)}}$ as follows:

for each $i\in[n_{1}]$. Here, $\tilde{\mu}$ and $\tilde{\sigma}^{2}$ are estimators of $\mu_{1}$ and $\sigma_{1}$, defined as

where $\tilde{\gamma}:=\kappa_{1}/(\tilde{v}+\lambda)$ and $\tilde{v}=n_{1}^{-1}\mathrm{tr}(\bm{I}_{n}-\bm{X}^{(1)}({(\bm{X}^{(1)}})^{\top}\bm{X}^{(1)}+n_{1}\lambda\bm{I}_{p})^{-1}{(\bm{X}^{(1)}})^{\top})$. These estimators are obtained by the observable adjustment technique described in Bellec (2022).

This index estimator $W_{i}$ is approximately unbiased for the index $\bm{\beta}^{\top}{\bm{X}_{i}^{(1)}}$, yielding the following asymptotic result:

We will provide its rigorous statement in Proposition 1 in Section 5.1.

There are other options for the pilot estimator besides the ridge estimator $\tilde{\bm{\beta}}$. If $\kappa_{1}\leq 1$ holds, the least squares estimator can be an alternative. If $y_{i}$ is a binary or non-negative integer, the MLE of logistic or Poisson regression can be a natural candidate, respectively, although the ridge estimator $\tilde{\bm{\beta}}$ is valid regardless of the form that $y_{i}$ takes. In each case, the estimated inferential parameters $(\tilde{\gamma},\tilde{\mu},\tilde{\sigma}^{2})$ should be updated accordingly. Details are presented in Section 6.

We develop an estimator of the link function $g(\cdot)$ using $W_{i}$ in (4). If we could observe the true index $\bm{\beta}^{\top}{\bm{X}_{i}^{(1)}}$ with the unknown coefficient $\bm{\beta}$, it would be possible to estimate $g(x)=\operatorname{\mathbb{E}}[y_{1}\mid\bm{\beta}^{\top}\bm{X}_{1}=x]$ by applying standard nonparametric methods to the pairs of responses and true indices $(\bm{y}^{(1)},\bm{X}^{(1)}\bm{\beta})$. However, as the true index is unobservable, we must estimate $g(\cdot)$ using given pairs of responses and contaminated indices ${(y_{i}^{(1)},W_{i})}_{i=1}^{n_{1}}$, where $W_{i}\approx\bm{\beta}^{\top}\bm{X}_{i}^{(1)}+\mathcal{N}(0,\tilde{\sigma}^{2}/\tilde{\mu}^{2})$. The type of error $\mathcal{N}(0,\tilde{\sigma}^{2}/\tilde{\mu}^{2})$ involving the regressor leads to an attenuation bias in the estimation of $g(\cdot)$, known as the errors-in-variables problem. To address this issue, we utilize a deconvolution technique (Stefanski and Carroll, 1990) to remove the bias stemming from error-in-variables asymptotically. Further details of the deconvolution are provided in Section C.

We define an estimator of $g(\cdot)$. In preparation, we specify a kernel function $K:\mathbb{R}\rightarrow\mathbb{R}$, and define a deconvolution kernel $K_{n}:\mathbb{R}\to\mathbb{R}$ as follows:

where $h_{n}>0$ is a bandwidth, $\mathrm{i}=\sqrt{-1}$ is an imaginary unit, $\tilde{\varsigma}^{2}=\tilde{\sigma}^{2}/\tilde{\mu}^{2}$ is a ratio of the inferential parameters, and $\phi_{K}:\mathbb{R}\rightarrow\mathbb{R}$ and $\phi_{\tilde{\varsigma}}:\mathbb{R}\to\mathbb{R}$ are the Fourier transform of $K(\cdot)$ and the density function of $\mathcal{N}(0,\tilde{\varsigma}^{2})$, respectively. We then define our estimator of $g(\cdot)$ as

where $\mathcal{R}[\cdot]$ is a monotonization operator, specified later, which maps any measurable function to a monotonic function, and $\breve{g}(\cdot)$ is a Nadaraya-Watson estimator obtained by the deconvolution kernel. We will prove the consistency of this estimator in Theorem 1 in Section 3.

The monotonization operation $\mathcal{R}[\cdot]$ on $\breve{g}(\cdot)$ is justifiable because the true link function $g(\cdot)$ is assumed to be monotonic. One simple choice for $\mathcal{R}[\cdot]$, applicable to any measurable function $f:\mathbb{R}\to\mathbb{R}$, is

This definition holds for all $x\in\mathbb{R}$. Another effective alternative is the rearrangement operator (Chernozhukov et al., 2009). This operator monotonizes a measurable function $f:\mathbb{R}\to\mathbb{R}$ within a compact interval $[\underline{x},\overline{x}]\subset\mathbb{R}$:

This operator, which sorts the values of $f(\cdot)$ in increasing order, is robust against local fluctuations such as function bumps. Thus, it effectively addresses bumps in $\breve{g}(\cdot)$ arising from kernel-based methods.

We next propose our estimator of $\bm{\beta}$ using $\hat{g}(\cdot)$ obtained in (8). In this step, we consider the link estimator $\hat{g}(\cdot)$ from $\bm{X}^{(1)}$ as given, and estimate $\bm{\beta}$ using $\bm{X}^{(2)}$. To facilitate this, we introduce the surrogate loss function for $\bm{b}\in\mathbb{R}^{p}$, with $\bm{x}\in\mathbb{R}^{p}$, $y\in\mathbb{R}$, and any measurable function $\bar{g}:\mathbb{R}\to\mathbb{R}$:

where $\bar{G}:\mathbb{R}\to\mathbb{R}$ is a function such that $\bar{G}^{\prime}(t)=\bar{g}(t)$. This function can be viewed as a natural extension of the matching loss (Auer et al., 1995) used in generalized linear models. If $\bar{g}(\cdot)$ is strictly increasing, then the loss is strictly convex in $\bm{b}$. Moreover, the surrogate loss is justified by the characteristics of the true parameter as follows (Agarwal et al., 2014):

provided that $G(\cdot)$ is integrable. The surrogate loss aligns with the negative log-likelihood when $g(\cdot)$ is known and serves as a canonical link function, thereby making the surrogate loss minimizer a generalization of the MLEs in generalized linear models.

Using the second dataset $\bm{X}^{(2)}$ with any given function $\bar{g}(\cdot)$, we define our estimator of $\bm{\beta}$ as

where $J:\mathbb{R}^{p}\to\mathbb{R}$ is a convex regularization function. Finally, we substitute the link estimator $\hat{g}(\cdot)$ into (13) to obtain our estimator $\hat{\bm{\beta}}(\hat{g})$. The use of a nonzero regularization term, $J(\cdot)$, is beneficial in cases where the minimizer (13) is not unique or does not exist; see, for example, Candès and Sur (2020) for the logistic regression case.

We finally study estimators for the inferential parameters of our estimator $\hat{\bm{\beta}}(\hat{g})$, which are essential for statistical inference as discussed in Section 1.1. As established in (3), it is necessary to estimate the asymptotic bias $\mu_{*}(\hat{g})$ and variance $\sigma_{*}^{2}(\hat{g})$ that satisfy the following relationship:

conditional on $(\bm{X}^{(1)},\bm{y}^{(1)})$ and consequently on $\hat{g}(\cdot)$.

We develop estimators for these inferential parameters using observable adjustments as suggested by Bellec (2022), in accordance with the estimator (13). For any measurable function $\bar{g}:\mathbb{R}\to\mathbb{R}$, we define $\bm{D}=\mathrm{diag}(\bar{g}^{\prime}(\bm{X}^{(2)}\hat{\bm{\beta}}(\bar{g})))$ and $\hat{v}_{\lambda}=n_{2}^{-1}\mathrm{tr}(\bm{D}-\bm{D}\bm{X}^{(2)}(({\bm{X}^{(2)}})^{\top}\bm{D}\bm{X}^{(2)}+n_{2}\lambda\bm{I}_{p})^{-1}({\bm{X}^{(2)}})^{\top}\bm{D})$ for $\lambda\geq 0$. When incorporating $J(\bm{b})=\lambda\|\bm{b}\|^{2}/2$ into (13) with $\lambda>0$, we propose the following estimators:

In the case where $J(\cdot)\equiv\bm{0}$ holds, we define

A theoretical justification for the asymptotic normality with these estimators and their application in inference is provided in Section 3.3.

As a preparation, we present some assumptions.

Assumption 1 requires that the split subsamples, as described at the beginning of Section 2, have the same order of sample size.

It is common to assume Gaussianity of covariates in the proportionally high-dimensional regime, as mentioned in Section 1.2. The condition $\bm{\beta}^{\top}\bm{\Sigma}\bm{\beta}=1$ is necessary to identify the scale of $\bm{\beta}$, which ensures the uniqueness of the estimator in the single-index model with an unknown link function $g(\cdot)$. For example, without this condition, it would be impossible to distinguish between $g(\bm{X}_{1}^{\top}\bm{b})$ and $f(2\bm{X}_{1}^{\top}\bm{b})$, where $f(t)=g(t/2)$, for any $\bm{b}\in\mathbb{R}^{p}$. Furthermore, the assumption of upper and lower bounds on the eigenvalues of $\bm{\Sigma}$ implies that $\|\bm{\beta}\|=\Theta(1)$.

Assumption 3 restricts the class of link functions to those that are monotonic. This class has been extensively reviewed in the literature, with Balabdaoui et al. (2019) providing a comprehensive discussion. It encompasses a wide range of applications, including utility functions, growth curves, and dose-response models (Matzkin, 1991; Wan et al., 2017; Foster et al., 2013). Furthermore, under a monotonically increasing link function, the sign of $\bm{\beta}$ is identified, so that we can identify $\bm{\beta}$ only by the scale condition $\bm{\beta}^{\top}\bm{\Sigma}\bm{\beta}=1$.

The continuity of $m_{2}(x)$ is maintained in many commonly used models, particularly when $g(\cdot)$ is continuous. For instance, the Poisson regression model defines $m_{2}(x)=\exp(x)(1+\exp(x))$, and binary choice models specify $m_{2}(x)=g(x)$.

We demonstrate the uniform consistency of the link estimator $\hat{g}(\cdot)$ in (8) over closed intervals. We consider the $m\mathrm{th}$-order kernel $K(\cdot)$ that satisfies

for $\ell\in[m-1]$. We then obtain the following:

Here, according to Fan and Truong (1993), the logarithmic rate $O_{\mathrm{p}}(1/(\log n_{1})^{m/2})$ reaches a lower bound, indicating that this rate cannot be improved.

This section demonstrates the marginal asymptotic normality of our estimator $\hat{\bm{\beta}}(\hat{g})$ for $\bm{\beta}$, facilitated by the estimators of the inferential parameters, $\hat{\mu}(\hat{g})$ and $\hat{\sigma}(\hat{g})$. These results are directly applicable to hypothesis testing and the construction of confidence intervals for any finite subset of the $\beta_{j}$’s.

We validate the normal approximation of the index estimator $W_{i}$ as shown in (6). For cases where $n>p$, we set $(n,p)=(500,50)$ for the Cloglog model and $(n,p)=(500,200)$ for the other models. For cases where $n<p$, we set $(n,p)=(250,500)$ and apply the ridge regularized estimator to all models. We assign the maximum likelihood estimator (MLE) of logistic regression to the pilot estimator for (i) Cloglog, the MLE of Poisson regression for (ii) xSqrt, and the least squares estimator for both (iii) Cubic and (iv) Piecewise models. We calculate $\tilde{\mu}(\bm{W}-\bm{X}\bm{\beta})/\tilde{\sigma}$ using $1,000$ replications for each setup.

Figure 1 displays the results. In all settings, the difference between the index estimator and the index follows a Gaussian distribution, as expected.

Next, we evaluate the numerical performance of the link estimator $\hat{g}(\cdot)$, constructed from $(W_{1},\ldots,W_{n})$, using a fixed bandwidth for each $n$. Figure 2 (left panel) shows that the estimation error of $\hat{g}(\cdot)$ for (iv) Piecewise uniformly approaches zero as the sample size increases. The right four panels in Figure 2 display the squared losses of $\hat{g}(\cdot)$ evaluated over the interval $[-3,3]$, which all decrease as $n$ increases.

We examine the asymptotic normality of each coordinate of the estimator $\hat{\bm{\beta}}(\hat{g})$ for the true coefficients. As in Section 4.2, we construct the estimator using a fixed bandwidth and apply the rearrangement operator $\mathcal{R}_{\mathrm{r}}[\cdot]$ as defined in (10) over the interval $[-3,3]$ to obtain $\hat{g}(\cdot)$. We then compute $\hat{\bm{\beta}}(\hat{g})$ according to (13) using $J(\cdot)\equiv\bm{0}$ when $n>p$ and $J(\bm{b})=\|\bm{b}\|^{2}$ when $n\leq p$. Figure 3 shows the marginal normal approximation of the estimators under these conditions. All histograms closely resemble the standard normal density, corroborating the asymptotic normality stated in Theorems 2 and 3.

Finally, we compare the estimation efficiency of the proposed estimator with several pilot estimators. We use the effective asymptotic variance $\sigma^{2}_{*}/\mu_{*}^{2}$ as an efficiency measure, which is the inverse of the effective signal-to-noise ratio as described in Feng et al. (2022). We estimate this variance using the statistic ${\hat{\bm{\beta}}^{\top}\hat{\bm{\beta}}}/({\hat{\bm{\beta}}^{\top}\bm{\beta}})-1$ for an estimator $\hat{\bm{\beta}}$. This statistic is a reasonable approximation of the asymptotic variance of the debiased version of $\hat{\bm{\beta}}$ and converges almost surely to the effective asymptotic variance under certain conditions (see Section 5 for details).